Molecular insights into kinetic stabilization of amorphous solid dispersion of pharmaceuticals

Abstract

Poor aqueous solubility of crystalline active pharmaceutical ingredients (APIs) restricts their bioavailability. Amorphous solid dispersions with biocompatible polymer excipients offer a solution to overcome this problem, potentially enabling a broader use of many drug candidate molecules. This work addresses various aspects of the in silico design of a suitable combination of an API and a polymer to form such a binary solid dispersion. Molecular interactions in such bulk systems are tracked at full atomic resolution within molecular-dynamics (MD) simulations, enabling to identify API–polymer pairs that exhibit the most beneficial interactions. Importance of these interactions is manifold: increasing the mutual miscibility, kinetic stabilization of their amorphous dispersions and impedance of the spurious recrystallization of the API component. MD tools are used to investigate the structural and cohesive properties of pure compounds and mixtures, with a special emphasis on molecular interactions, microscopic structures and internal dynamics. This analysis is then accompanied by a macroscopic image of the energetic compatibility and vitrification tendency of the mixtures in terms of their excess enthalpies and glass transition temperatures. Density-functional theory (DFT) and non-covalent interaction (NCI) analysis fortify our computational conclusions and enable us to map the intensities of particular NCI among the individual target materials and relevant molecular sites therein. Three archetypal polymer excipients and four API molecules are included in this study. The results of our computational analysis of molecular interactions in bulk systems agree with the experimentally observed trends of solubility of the given API in polymers. Our calculations confirm PVP as the most potent acceptor of hydrogen bonding among the three considered polymer excipients, whereas ibuprofen molecules are predicted to be the most efficient hydrogen bond donors among our four target APIs. Our simulations also suggest that carbamazepine does not exhibit particularly strong interactions with the considered polymer excipients. Although current MD cannot offer quantitative accuracy of many of the discussed descriptors, current computational models focusing on NCI of APIs with polymer excipients contribute to understanding of the behavior of these materials at the molecular level, and thus also to the rational design of novel efficient drug formulations.

---

1. Introduction

Crystalline forms of active pharmaceutical ingredients are advantageous because of their better stability during long-term storage, reproducibility of the crystallization procedure and more reliable predictions of material properties at the molecular level under defined conditions.¹ Concurrently, a vast majority of newly proposed API candidate molecules have a predominant hydrophobic character, imparting a poor solubility of their crystals in water,² which adversely limits the bioavailability, dissolution, and distribution of such APIs through the treated organism. This fact impedes potential wider use of numerous API candidates as solid drugs in medical treatment. Improving the efficiency of oral administration of poorly soluble APIs has become one of the biggest challenges in the pharmaceutical industry.³ Strategies for rendering the otherwise promising API candidates more soluble have to be searched for to overcome these limitations. These rely on a suitable chemical functionalization, conversion of an API to its salts,⁴ formulation of cocrystals,⁵ or preparation of amorphous forms of APIs.^3,6 The latter option seems very beneficial as any amorphous phase inherently exhibits a larger solubility than respective crystals due to its thermodynamic metastable nature. However, kinetic aspects of stabilization of the amorphous API forms with respect to its crystallization over sufficient periods of time need to be also considered.^7,8 During processing, storage, and after contact with water or humidity, the thermodynamically metastable amorphous forms tend to recrystallize, which is a price to be paid for exploiting the higher solubility of a metastable API formulation.⁹ Formulation of amorphous solid dispersions of API in various matrices that stabilize the amorphous API form relies predominantly on biocompatible polymers which seems very promising in this context.¹⁰

Half a century ago, the earliest account of the so-called first-generation solid dispersion was provided by exploring eutectic mixtures that exhibit an enhanced rate of API release and bioavailability.¹¹ First-generation solid dispersions were built from crystalline carriers such as urea or sugars, forming crystalline solid dispersions. The second generation of solid dispersions was then based on replacing crystalline carriers with amorphous carriers such as polymers, forming an amorphous product in which APIs are dissolved. There is also a third generation of solid dispersions using a surfactant carrier or a combination of amorphous polymers and surfactants.⁶ The motivation of the latter strategy is to kinetically stabilize the amorphous solid phases of APIs by beneficial molecular interactions between the API and its excipient, which averts the rearrangement of API molecules into a crystal lattice even at supersaturated API concentrations.^7,8 This spurious process may even begin with an amorphous–amorphous separation, followed by a separate recrystallization of the API from its pure phase.¹²

Combinatorial chemistry techniques and high-throughput screening have led to an immense increase in the quantity of proposed poorly soluble candidate API molecules. Concurrently, there are vast numbers of known polymer materials that can be considered biocompatible and that exhibit the potential for acting as an excipient for an API. There are accurate and reliable calorimetric pathways to systematically characterize such complex multi-component systems.¹³ These purely experimental approaches become, however, very tedious when high-throughput screening is due. To rationalize the design of viable amorphous formulation for a particular API, in silico approaches are inevitable to avoid the laborious trial and error scheme throughout the experimental stages of drug development.¹⁴ Molecular simulations enable to focus on individual molecular interactions and mechanistic degrees of freedom in bulk materials at an all atom resolution.¹⁵ Such insights are invaluable for identifying the most important molecular interactions that a real API molecule can exhibit and preselecting a short list of potential excipients that are capable of forming beneficial interactions and spatial packing to stabilize the amorphous API dispersion.

The glass-transition temperature (T_g) of a bulk amorphous system is a key parameter governed by its inner dynamics and molecular interactions. Robust approaches for T_g predictions from MD simulations have been developed and validated recently.^16–19 These approaches rely on analysis of the simulated density as a function of temperature which can be extracted from molecular simulations with a fairly low uncertainty. Resulting predictions of the glass transition do not exhibit quantitative accuracy,¹⁹ but represent a valuable macroscopic descriptor of the internal dynamics of complex amorphous materials. As such, analysis of the simulated molecular interactions at the microscopic level and the glass transition temperatures at the macroscopic scale go hand in hand when designing suitable API–excipient pairs and assessing their kinetic stability.¹⁵

In this work, we select three archetypal biocompatible polymers, namely polyethylene glycol (PEG),²⁰ polylactic acid (PLA),^21–23 and polyvinyl pyrrolidone (PVP)^24–31 that all have (at least as building blocks of the true polymer) a real significance to be used as excipients in API formulations. Monomer units of all these materials are depicted in Fig. 1. To enable validation of the simulated results, we select common API molecules, the amorphous solid dispersions of which have been prepared and experimentally characterized recently.^13,23,30,31 Namely, we consider carbamazepine (CBZ),²⁹ ibuprofen (IBU),^21,23 indomethacin (IND),^22,25–30 and naproxen (NAP).^13,24 Molecular formulae of the target API are given in Fig. 1.

Fig. 1 Monomer units of target polymers (top row: PLA, PEG, and PVP) and API molecules (middle row: carbamazepine, ibuprofen, and bottom row: indomethacin, naproxen) considered in this work.

Among these materials, there are recent experimental data on equilibrium solubilities of APIs in the PLA based polymers, which range to low units of weight percent at 298.15 K (however, corresponding to dozens of molar percent due to significant differences in molecular weights of the components). Experimental literature data reveal that mixtures with composition around molar fraction x_API = 0.85 represent slightly oversaturated mixtures in terms of the API concentration.²³ Several attempts at preparation of even more oversaturated mixtures with x_API exceeding 0.90 have been exerted, leading in some cases to two-component systems that remain homogenous over limited time periods before the phase separation occurs.^23,30

Solubility of the current target API in PLA based polymers was found to increase in a row: IND < NAP < IBU,^21–23,31 assuming that PLA and its copolymer PLGA exhibit very similar behavior. Additionally, the solubility of IBU, IND, and NAP was experimentally found to vary significantly among the polymers studied, with PVP exhibiting a notably higher solubilizing capacity compared to PLA and PEG.^{13,20–22,24,27,31,32} These hierarchies indicate some variance in the molecular and thermodynamic compatibility of the components or of the extent of kinetic stabilization of such oversaturated mixtures with respect to recrystallization of the API.²³ A comparison of the rankings of thermodynamic and kinetic factors with the observed long-term stabilities points to the importance of looking at these aspects simultaneously. For instance, IBU was proved to be more soluble in PLA while IND experiences a large kinetic stability of its oversaturated dispersions in that polymer.²³ At this point, the path of atomistic simulations is adopted to provide a detailed interpretation of the observed phenomena.

Classical equilibrium MD simulations with a full resolution of atoms are performed in this work to investigate the microscopic structure, spatial packing, energetics and prevalence of individual molecular contacts and interactions and compare all these characteristics between pure amorphous materials and their binary mixtures. In addition, non-equilibrium simulations of rapid quenching of liquid mixtures to amorphous glassy phases are performed to assess the macroscopic kinetic stability of the amorphous mixtures and to interpret its trends among various target materials at the level of individual molecular interactions. Despite the inner dynamics of bulk amorphous materials can be captured more accurately with explicitly polarizable simulations,^19,33 we limit the scope of this work to non-polarizable simulations which impart a fair compromise between the accuracy of results and associated computational costs and a DFT analysis of non-covalent interactions within representative molecular clusters.

On purpose, we selected for our simulations systems with compositions of the two-component mixtures that slightly exceed the experimentally observed solubility of the target API in the given polymers. We aim at an atomistic interpretation of the experimentally observed solubility trends and diverse kinetic stabilities of various oversaturated mixtures that are prone to a phase separation from the thermodynamic point of view. Simulations of the local structure, inner dynamics and molecular interactions within the two-component mixtures are performed for this purpose. Note that the goals of this work do not include in silico predictions of the equilibrium solubility of the API.

2. Computational methods

2.1. General simulation setup

Polymer chains were created from dimer units by their replication on a one-dimensional lattice as described in our previous work.³⁴ Such a procedure leads to extended linear chains in all-trans conformation that need to be folded before packing in a simulation box. PLA chains were prepared from 200 monomer units of D-lactic acid so that the weight of a single PLA chain exceeds 14 400 g mol⁻¹. PEG chains were prepared from 239 monomer units of ethylene glycol, corresponding to a molecular weight over 10 500 g mol⁻¹. Both PLA and PEG backbone chains were terminated with hydroxyl moieties at both ends. PVP chains were created from 90 vinyl pyrrolidone units, corresponding to a molecular weight over 10 000 g mol⁻¹. Methyl groups were selected to terminate the PVP backbone. To simplify the computational setup, all polymer samples were simulated as monodisperse systems with the given chain lengths, since the current level of MD simulations was demonstrated not to distinguish³⁴ the differences of thermodynamic properties that occur in reality among various polymer samples differing in their molecular mass or polydispersity.¹³ Literature OPLS force field models for PLA,³⁵ PEG,^36–38 and PVP,^39,40 which had been united to the OPLS format⁴¹ and validated previously in another our study,³⁴ were used to describe the behavior of the polymer molecules.

For the sake of compatibility with polymer treatment, OPLS force fields were adopted to describe the behavior of API molecules as well. Namely, literature parameter sets were used for CBZ,¹⁵ IBU,⁴² IND,⁴³ and NAP.⁴⁴ If needed, these force-field models have been also unified to the OPLS format and subsequently tested in our earlier work.¹⁵ Note that we selected the OPLS force field⁴¹ as a simple non-polarizable, yet versatile all-atom model that has been previously tested to provide a fairly accurate description of both crystal and liquid phases of many molecular materials^41,45–47 and their mixtures.^48–50 In the context of the all-atom OPLS model for pure API and polymer materials that are targeted in this work, our previous studies evaluated the mean errors of the MD-simulated bulk liquid densities of the API to range to 4–6% at 400 K,¹⁵ and of the bulk monodisperse polymers with molar mass around 10 kg mol⁻¹ to be 7–15%.¹⁹ Fusion enthalpies of the pure API were captured within 10–50% of the experiment.¹⁵ Typical errors in simulated glass-transition temperatures amount to 60–80 K for the pure API¹⁵ and to 40–80 K for pure polymers³⁴ when the given non-polarizable model is used. Such performance of the model, however, predetermines that quantitative predictions of subtle thermodynamic properties cannot be targeted by the current simulations. Instead, this work focuses on a qualitative interpretation of the previously observed phenomena.

LAMMPS software package⁵¹ (version 5 May 2020) was used for all MD simulations. At first, a single isolated linear polymer chain in its all-trans conformation was simulated at 300 K for 2 ns. This procedure was important to bend the linear extended chain conformations to a more reasonable (but rather random) globular conformation to be packed later in simulation boxes mimicking bulk mixtures.³⁴ Subsequently, simulation boxes for amorphous binary mixtures of APIs with polymers were prepared using the Packmol code,⁵² which places a user-defined number of molecules at random non-overlapping positions within a cubic box. Molecular topology and input files for the simulations in LAMMPS were then generated using the fftool script.⁵³ Compositions of the simulation boxes used in this work are summarized in Table 1.

Table 1 Composition of simulation boxes (x and w being the molar and mass fractions, respectively) mimicking amorphous dispersions of API and PLA considered in this work

PLA–API systems
N API	100	200	300
N PLA	17	17	17
x API	0.85	0.92	0.95
w CBZ	0.09	0.16	0.22
w IBU	0.08	0.14	0.20
w IND	0.13	0.23	0.30
w NAP	0.09	0.16	0.22
POLY-IBU systems
POLY	PEG	PVP
N POLY	18	20
N IBU	100	100
x IBU	0.85	0.83
w IBU	0.10	0.09

---

A single pre-folded globular conformation, described above, was always packed in a simulation box to facilitate building up the simulation boxes mimicking the bulk mixtures. In our previous work,³⁴ we found that equilibration of bulk polymer at a sufficiently high temperature (e.g. 500 K for PLA and PEG and 800 K for PVP) leads to a complete erasing of the conformational memory of individual chains in our simulations. Properties such as the density, total energy or even the mean square end-to-end distance of a polymer chain proved to be insensitive to the imposed initial conformation when this high-temperature equilibration phase was included in the simulation protocol.

Equilibrium simulations of bulk systems were always simulated with isotropic periodic boundary conditions within the isothermal–isobaric NpT ensemble, adopting the Nosé–Hoover thermostat and barostat⁵⁴ and the velocity Verlet integrator as implemented in LAMMPS.⁵¹ A 12 Å cut-off distance was imposed for the short-range interactions, being accompanied by integral tail corrections for interactions beyond this threshold. A long-range PPPM solver was employed to sum the Coulombic interactions of atomic point charges.⁵⁵ The SHAKE algorithm⁵⁶ was applied to constrain the lengths of covalent bonds terminating in hydrogen atoms, which enabled us to perform the simulations with a constant time integration step 1 fs.

Following our search for an efficient computational treatment of polymer samples,³⁴ equilibration simulation runs of bulk systems were initiated from randomly packed simulation boxes and were equilibrated at the temperature of 500 K and pressure of 1 bar over 2 ns. Subsequently, these pre-equilibrated systems were adjusted to additional two target temperatures, 300 K and 800 K, by an additional 2 ns equilibration period at both temperatures. Finally, production runs spanning 10 ns at both 300 K and 500 K were performed to sample the equilibrium properties of the pure compounds and their binary mixtures. Enthalpic, structural (radial distribution functions, RDFs) and dynamic (mean square displacement, MSD) characteristics of the simulated ensembles were sampled every 1 ps. Block-averaging scheme was employed to estimate the statistical sampling uncertainty of the observed properties, such as equilibrium density and enthalpy.⁵⁷ Static structure factors corresponding to small-angle neutron scattering (SANS) were calculated for all the considered mixtures. For this purpose, the simulated trajectories (10 ns period) were post-processed by the diffraction utility,⁵⁸ using coherent scattering lengths for elements from ref. 59. Such computational setup was previously validated for bulk nano-structured materials.^19,60

MD results for pure amorphous APIs¹⁵ and pure polymer systems³⁴ were taken from our recent studies validating the computational setup of atomistic MD simulations of structure and glass transitions in these materials. The computational setup adopted to treat pure materials in the MD simulations in these studies was identical to that adopted in this work.

2.2. Localization of glass transition temperatures

To determine the glass transition temperature of the mixtures, rapid non-equilibrium quenching of the bulk was simulated.³⁴ A cooling rate of 30 K ns⁻¹ was imposed to cool the system from 800 K down to 200 K. To minimize the hazards associated with random freezing of the simulated amorphous system in an anomalous glassy state and to gain an insight on the actual computational uncertainty of this quenching procedure, it was repeated five times for each system.¹⁷ Initial configurations of the simulation boxes for these five simulation replicas were selected from mutually distant time steps from the previous trajectory at 800 K, randomizing also the initial atomic velocities for each of these cooling runs to maximize the independent character of these replicas.³⁴

Density as a function of the instantaneous temperature was sampled every 1 ps during these non-equilibrium quenching runs. Further data processing was then limited to the temperature interval from 600 K to 200 K. Hyperbola-fit method was applied to determine the glass transition temperature of all simulated bulk mixtures.⁶¹ We adopted the following equation to fit the obtained density–temperature data:

where a, b, and c are adjustable parameters, ρ₀ is a reference density and T₀ is the actual glass transition temperature. This procedure has been demonstrated to minimize human intervention in the analysis of the simulated density trend shift and can be efficiently automated.^16,61 It can be used also in situations where the glass transition is not accompanied by an abrupt trend shift of the density with respect to temperature. Note that eqn (1) provides a smooth curve when the parameter c amounts to finite numbers (with T₀ corresponding to the hyperbola center), whereas it becomes a broken line for c = −∞ (with the break at T₀). Both settings (finite or infinite c values) typically produce very similar T_g results when the density trend shift is clearly visible.⁶¹ The latter setting proves to be more useful for cases with a more indistinct density trend shift due to glass transition.³⁴ For most bulk simulated systems and proper simulation protocols, it holds that both T₀ values are very similar so that one can assume that it holds T₀ = T_g.

2.3. Identifying important non-covalent interactions

To localize and quantify the non-covalent intermolecular interactions (NCIs) between the API and its polymer carrier, additional analyses were performed within the framework of quantum theory or atoms in molecules.⁶² We focused on the critical points of the reduced density gradient (RDG) versus the sign of the second-density Hessian eigenvalue (λ₂). Note that the RDG critical points correspond to positions around the interacting molecules, where the RDG reaches or at least sharply approaches the zero value, and concurrently these sites match regions where most intense non-covalent interactions among individual atoms occur.⁶³

For this purpose, a cluster of a single API molecule and an oligomer containing four to six monomer units (so that the chain covers the key carboxyl or amide moiety in the API molecule) were prepared. A simulated annealing protocol, cooling slowly this cluster from 400 K to 0.01 K over 1000 ns, was performed at the level of non-polarizable MD to find a reasonable stable cluster conformation.⁶⁴ Geometry of such an annealed cluster was then reoptimized within the DFT framework using the B3LYP-D3(BJ)/6-311+G(d,p) level of theory^65,66 in Gaussian, version 16 (RevB.01).⁶⁷ The MultiWFN package⁶⁸ was used to perform all manipulation of the quantum-chemically computed electron structure properties.

To analyze optimum molecular contacts due to the hydrogen bonding of API molecules and fragments of the polymer chains, clusters containing a single API molecule and a capped chain monomer or a tetramer were prepared. Geometries of these clusters were obtained from a classical-MD simulated annealing procedure followed by a density-functional theory B3LYP-D3/6-311+G(d,p) optimization. Finally, scaled zeroth order symmetry-adapted perturbation theory sSAPT0⁶⁹ along with the calendar jun-cc-pVDZ basis set⁷⁰ were used for calculations^69,70 of their pair interaction energies. There calculations were performed in PSI4⁷¹ to dissect individual mechanistic components of the interaction and to estimate the hydrogen-bonding strength by subtracting the dispersion interaction from the total SAPT pair interaction energy. Note that the SAPT approach was selected for calculations of the pair interactions as it is inherently free of the basis-set superposition error.⁷²

3. Results and discussion

This section is organized as follows. First, excess molar volumes and enthalpies of the API–polymer mixtures are discussed to present an initial macroscopic picture of the character of the mixtures. Interactions and internal dynamics are then analyzed in more detail at the microscopic level. Finally, simulated glass transition temperatures are presented to unite the microscopic and macroscopic characterization of the mixtures and their kinetic stability.

3.1. Excess properties of mixtures

Excess quantities express deviations of thermodynamic properties of real mixtures from ideal mixtures, i.e. the extent of asymmetry of heteromolecular interactions with respect to the situation in pure compounds. In the current case of amorphous mixtures of molecules with dramatic differences in molecular size and capability to interact via hydrogen bonding, appreciable non-zero excess quantities can be expected.

First, excess molar volumes V^E of API–PLA mixtures are negative for all considered API at 300 K as depicted in Fig. 2 and listed in Table S1 (ESI†). This indicates a good steric compatibility of API molecules with PLA chains. The spatial packing of their mixtures is then slightly more efficient than that of pure compounds. Particular V^E values in Fig. 2 and Fig. S1 (ESI†) correspond to less than 1% of the actual molar volumes of the mixtures with PLA and up to 2% of those with PVP in the absolute value at 300 K. Naturally, the magnitude of this effect increases as the API molar fraction decreases from unity. The most important negative deviation from ideality is exhibited by ibuprofen mixed with PLA, whereas the analogous value for naproxen is lower by a factor of two.

Fig. 2 Excess molar volumes V^E (in cm³ mol⁻¹) of simulated mixtures of target API and selected polymers at various compositions and temperatures 300 K (left box, full symbols) and 500 K (right box, empty symbols).

At 500 K, V^E values are higher for all considered mixtures than at 300 K. The only API–PLA mixture that retains at least slightly negative V^E at 500 K is IND–PLA. Otherwise, simulated V^E values are effectively zero as their magnitudes in the vicinity of the zero level are comparable with their statistical sampling uncertainties. V^E is convincingly positive at 500 K for the NAP–PLA mixture. Such results indicate that the higher temperature significantly disrupts the beneficial packing that the API molecules can exhibit within the PLA matrix. From the perspective of the spatial packing of the bulk, carbamazepine and ibuprofene form near-to-ideal mixtures with PLA at 500 K. Notably, IBU–PLA mixtures, exhibiting the most negative V^E at 300 K among the considered API, exhibit also the largest increase of their V^E upon heating to 500 K (by almost 10 cm³ mol⁻¹ at x_API = 0.85), which may be an indication for limited thermal stability of the respective beneficial ibuprofen–PLA interactions. Predicted V^E values decrease less significantly for the other systems (6 to 8 cm³ mol⁻¹ at x_API = 0.85).

Despite the beneficial packing features of the discussed API–PLA mixtures under ambient conditions, their excess enthalpies H^E are rarely negative, as depicted in Fig. 3. It is only ibuprofen that appears to form significant beneficial interactions with PLA above the level of an ideal mixture. Whereas the computed H^E values are nearly zero for CBZ–PLA mixtures, analogous values are appreciably positive when PLA is mixed with naproxen and indomethacin, indicating certain deficiencies of these mixtures from the perspective of molecular interactions. Interestingly, the predicted H^E values are higher at 500 K, exhibiting positive values for all the systems. This may be an impact of the loss of the spatial packing advantage that manifests at 300 K and that can partly compensate for the disruption of strongly attractive API – API interactions that occur in bulk preferentially at high API concentrations. The presence of polymers in a bulk mixture attenuates the API – API interactions, which are not fully compensated by creating new intermolecular interactions between API and PLA in the mixtures, there is minimum to no enthalpic driving force for mixing of the given APIs with PLA. Should such a mixing leading to homogenous two-component phases occur, it would have to rather rely on entropic factors.

Fig. 3 Excess molar enthalpies H^E (in kJ mol⁻¹) of simulated mixtures of target API and selected polymers at various compositions and temperatures 300 K (left box, full symbols) and 500 K (right box, empty symbols).

The trends of calculated H^E, which are appreciably higher for mixtures with lower API fractions, support this hypothesis further. Notably, trends of H^E simulated among individual APIs qualitatively match the ranking of their experimental solubilities in PLA based polymers: IND < NAP < IBU. This is an encouraging finding, calling for additional structural analyses at the microscopic scale.

An explanation for the small magnitudes of H^E for CBZ–PLA mixtures at ambient temperature is rooted in compensation of two effects: (i) an endothermic impact of the partial disruption of the homo-molecular CBZ–CBZ hydrogen bonds due to the polymer presence in the mixture which will be analyzed below; and (ii) an exothermic effect due to the improvement of the spatial packing efficiency of the bulk mixture leading to more intense dispersion cohesive interactions which is supported by the negative excess volume shown in Fig. 2 for the CBZ–PLA mixture at 300 K. Larger (positive) H^E values for CBZ–PLA mixtures at 500 K are then an imprint of the loss of the packing efficiency benefits which propagate to the endothermic mixing.

Excess properties of mixtures of IBU with individual target polymorphs are depicted in Fig. S1 (ESI†). Both excess molar volumes and enthalpies of systems containing PVP are significantly more negative than for the remaining two polymers, which indicates that strong attractive IBU–PVP interactions are established upon mixing which further stabilizes that mixture. On the other hand, the IBU–PEG system exhibits excess properties very close to zero or even positive already at 300 K, indicating only weak hetero-molecular interactions in that system. PLA remains in the middle of the current polymer hierarchy with respect to the excess properties. Such a trend again agrees with the experimental hierarchy of IBU solubilities in individual polymers which was established above.

3.2. Equilibrium structure and interactions

Recently, super-saturated mixtures of some of the target API and polymers have been experimentally prepared and demonstrated to remain homogenous for various macroscopic time periods, ranging from minutes (e.g. for IBU–PLA) to months (e.g. IND–PLA).²³ To validate the structural homogeneity of our simulated samples at the microscopic level, we characterized our samples from the perspective of their nano-structuration via simulations of their SANS diffraction patterns. Fig. S2 (ESI†) depicts static structure factors S(q) mimicking SANS outputs. The only one convincing broad Bragg peak near the reciprocal coordinate q = 1 Å⁻¹ value in all such plots can be attributed to the adjacency of individual molecules in the amorphous bulk phase. The position of this signal corresponds then to an efficient radius of the API molecules or blobs in the polymer chain which amounts to around 2π/q ≈ 6 Å. Naturally, S(q) signals are smooth at 500 K as all the samples are fluid at these conditions, whereas S(q) signals are relatively noisy at 300 K, corresponding to the loss of fluidity of the glassy systems. Since there are no significant pre-peaks on the S(q) signals that would correspond to a regular formation of larger domains, one can rule out immediate nano-segregation of the API molecules from the polymer carrier over the simulated ten-nanosecond time horizon. Simulated samples can thus be assumed to be fairly homogenous over the time scale of the generated MS trajectories.

Comparison of simulated SANS profiles for mixtures of IBU with individual polymers, depicted in Fig. S3 (ESI†), reveals a qualitative similarity of the mixtures regardless of the present polymer. The position of the main broad peak of those SANS profiles varies among the present polymers. It is shifted to the largest q values for the IBU–PEG mixture which can be attributed to the simplest chain structure of PEG, and thus to the smallest size of its chain molecules. There is a very weak shoulder near 0.8 Å⁻¹ for the IBU–PVP mixture, possibly indicating some weak medium-range ordering of that amorphous mixture with a characteristic length of roughly 7.9 Å. Interestingly, experimental X-ray diffraction on a similar IBU–PVP amorphous mixture revealed a broad shoulder of the diffraction starting at roughly 0.7 Å⁻¹.⁷³ Our simulations predict that there are no such S(q) pre-peak features for the IBU–PEG or IBU–PLA mixtures that would exceed the computational noise.

To investigate how the ordering of the amorphous mixtures differs from what can be observed for pure amorphous API, Fig. S4 (ESI†) depicts simulated SANS profiles for pure amorphous API at 300 K. Pure amorphous IBU exhibits two distinct pre-peaks at 0.35 Å⁻¹ and 0.51 Å⁻¹ (corresponding to characteristic lengths of 18.0 Å and 12.3 Å, respectively) that can be interpreted as imprints of an alternation of hydrogen-bonded layers where the polar and non-polar segments of IBU molecules tend to segregate in the bulk. Note that the linear dimension of the closest carboxyl-dimer of IBU molecules (depicted in Fig. S5, ESI†) amounts to roughly 18 Å which reasonably agrees with the position of the pre-peak in the simulated SANS profiles. Note that a significant S(q) shoulder has been also observed experimentally for amorphous IBU under ambient conditions,⁷⁴ confirming our simulated results. Since such significant pre-peaks are missing in the simulated SANS profiles of mixtures of IBU with any of the polymers in Fig. S2 (ESI†), one can conclude that formulation of these API–polymer dispersions is a successful strategy to disrupt the most important cohesive features that occur in neat API phases.

There is also experimental evidence of similar structural features of amorphous IND, the diffraction pattern of which exhibits a shoulder at 0.82 Å⁻¹ (corresponding to a characteristic length 7.8 Å),^75,76 which can be also attributed to some medium-range ordering of the hydrogen-bonded IND dimers in the bulk. Fig. S5 (ESI†) depicts the scale of the closest carboxyl-dimer of IND molecules that agrees with this characteristic length. Our simulations also indicate a very weak shoulder in that S(q) region. Simulated S(q) profiles for pure amorphous NAP are qualitatively similar to those of IND. Finally, there is also an experimental report of a diffraction experiment on an amorphous pure CBZ phase, that lacks any pre-peaks,⁷⁷ which is in agreement with our simulations predicting a smooth decay of the S(q) in the low-q region.

RDF represents a suitable structural quantity for interpretation of the observed behavior of excess quantities of the mixtures and to identify important close atomic contacts and related non-covalent interactions that have the highest impact on material cohesion. Since all considered API molecules are capable of the formation of strong hydrogen bonds in both their condensed phases, and all the target polymers contain oxygen atoms that behave as acceptors of hydrogen bonds, this type of non-covalent interaction naturally needs to be analyzed.

At first, homo-molecular API–API interactions in the mixtures will be discussed and their occurrence will be compared to what occurs in neat amorphous API.

Fig. 4 depicts homo-molecular RDF simulated for pure API and their mixtures with PLA. Large amplitudes and narrow shapes of the first RDF peak indicate that strong hydrogen bonding dominates the cohesive interaction in the pure amorphous API phases. Current RDF data suggest that CBZ–CBZ hydrogen bonding is the weakest in bulk amorphous CBZ at 300 K among the target API, whereas IND–IND interactions are the strongest. This can be glimpsed also from the positions of the first RDF peaks, which are located at 0.5 Å shorter molecular separations for IND than for CBZ. The second RDF peak that is exhibited by CBZ systems at around 3.5 Å corresponds only to the stoichiometry of the NH₂ moiety, not to formation of an immediate second solvation shell.

Fig. 4 Radial distribution functions g(r) for homo-molecular contacts of the API carboxyl/amide hydrogen atom (H_O/H_N) with the API carbonyl oxygen atom (O_C) in pure amorphous API and API–PLA mixtures simulated at 300 K. Signals from the mixtures are renormalized so that the unity value always corresponds to the average contact counts in neat bulk API. Inset molecular fragments serve as tentative illustrations of the possible close contacts in bulk API.

For the mixtures, all RDF signals are normalized so that the g = 1 value always corresponds to the average particle density in neat amorphous API. Thanks to this normalization, one can directly observe how the dilution of the API by dispersing it in the polymer excipient affects the hydrogen bonding structure. For CBZ, IBU and NAP, there is a significant decrease in the RDF peak amplitudes upon API dilution with PLA (by a factor from 4 to 8 for x_API = 0.85), whereas this amplitude drop is less pronounced for the IND–PLA mixtures (by a factor of 1.5). Such data indicate that molecules of the former three APIs can be considered as fairly dispersed in the PLA matrix and that their homo-molecular interactions that impart the strong cohesion of their pure crystals are largely suppressed in the mixture. On the other hand, IND molecules remain significantly bound among themselves even in the mixture with PLA. This strong hydrogen bonding among IND molecules still leaves some potential for their aggregation in the bulk mixtures, possibly leading to structural inhomogeneity on longer time scales. Results presented in Fig. 4 itself cannot distinguish between formation of dimers or catemer chains among API molecules in the bulk. Therefore, an additional analysis of this phenomenon will be included below.

Analogous homo-molecular API–API RDF functions extracted from our MD simulations at 500 K are depicted in Fig. S6 (ESI†). Amplitudes of all RDF peaks naturally diminish as the system becomes more chaotic. Ratios of the first-peak amplitudes at 500 K and 300 K range from 0.28 (IBU) over 0.45 (CBZ) and 0.50 (NAP) to 0.93 (IND) at x_API = 0.85.

Such observation of the RDF can be naturally supported by evaluation of the integral coordination numbers over the first coordination shell in the simulated bulk systems, this time normalized always with respect to the actual average density of each simulated system. Fig. 5 confirms the dramatic losses of homo-molecular API–API hydrogen bonding upon mixing with PLA for all considered APIs except IND. Notably, relevant carboxyl moieties are saturated with hydrogen bonding from two-thirds (for CBZ and NAP) over roughly 80% (for IND) up to nearly 90% (for IBU) at 300 K in API glasses. Data in Fig. 5 depict that this homo-molecular saturation drops to near 10% at x_API = 0.85 and 300 K for all API except IND, which retains almost 60% of its carboxyls hydrogen-bonded with other IND molecules at these conditions. Interestingly, PLA mixtures with IND and NAP exhibit an almost linear relationship of the API–API coordination number with the API concentration in the considered range at 300 K, but the coordination number decreases more steeply for the CBZ and IBU containing mixtures.

Fig. 5 Coordination numbers for intermolecular contacts of the API carboxyl/amide hydrogen atom (H_O/H_N) with the API carbonyl oxygen atom (O_C) in pure amorphous API and API–PLA mixtures simulated at 300 K (full columns) and 500 K (empty columns).

Fig. 5 also compares how these coordination numbers change upon heating all bulk phases to 500 K. This significant heating has no visible impact on all IND systems, confirming the resistance of the homo-molecular IND–IND hydrogen bonding to the temperature-induced decay, which is a sign of its massive strength.^46,50 A slight decrease in the API–API hydrogen bonding can be glimpsed for NAP systems at 500 K, whereas this heating leads to a drop to 60% of the intensity at 300 K for CBZ and IBU and similarly for their mixtures with PLA.

To assess the character of hydrogen-bonded homo-molecular clusters of API molecules in the bulk, one can consider the dominant cohesive features that occur in pure crystalline API as a reference. These represent the closest carboxyl dimers, bound by two short hydrogen bonds, in pure ibuprofen⁷⁸ and indomethacin⁷⁹ (with the H_O⋯H_O contact distance around 2.2 Å as depicted in Fig. S7, ESI†) whereas enantiopure naproxen molecules are packed in hydrogen-bonded catemer chains in their crystals (with the H_O⋯H_O contact distance above 3.1 Å).⁸⁰ Carbamazepine molecules then form amide dimers in their pure crystals.⁸¹

To gain an insight on these packing motifs in our target amorphous systems, we analyzed RDF that correspond to mutual contacts of two hydroxyl or amide hydrogen atoms in adjacent API molecules at 300 K, which are depicted in Fig. S8 (ESI†). These results indicate that both carboxyl packing motifs are present in all considered systems, but their occurrence varies. The closest carboxyl dimer with H_O⋯H_O contact distances below 2.5 Å is dominant only in systems with high NAP concentrations and its relative abundance wanes upon dilution with PLA. Both the closest dimer and the chained dimer are similarly populated in IBU and CBZ systems, with the dimer populations being significantly lower for CBZ in general. Interestingly, IND is predicted to form predominantly the catemers with H_O⋯H_O contact distances covering a wide range from roughly 3.0 to 4.5 Å which indicates also a variety of conformations of these chained IND clusters. Coordination numbers corresponding to these H_O⋯H_O contacts reveal that up to tetramers are formed in pure amorphous IBU as there are 2.7 H_O atoms within 5 Å from a reference H_O atom on average. Pure IND and NAP are predicted to form rather trimers (H_O⋯H_O coordination number 2.1). Bulk CBZ is modelled to contain dimers and trimers of the amide moieties as the respective H_N⋯H_N coordination number within 5 Å is 3.0, which corresponds to 1.5 coordinated amide group for each reference H_N atom. These coordination numbers drop down to 0.3–0.4 when RDF for analogous H_O⋯H_O contacts in mixtures are analyzed. The only exception is the INC–PLA mixtures, where significant IND–IND clustering persists even in the mixtures (coordination number 1.4).

To conclude this discussion of the simulation results on homo-molecular API–API hydrogen bonding, it appears to be fairly resistant to both temperature- and dilution-induced decay in IND-containing amorphous systems. NAP systems are relatively resistant to temperature-induced decay but are prone to dilution-induced decay, whereas CBZ and IBU-containing systems are prone to both effects. Analysis of the IBU–IBU coordination numbers related to their homo-molecular hydrogen bonding among the mixtures of IBU with individual polymers revealed a trend increasing in a: PVP (0.006) < PEG (0.020) < PLA (0.107) at 300 K. This high capability of PVP to disperse IBU in its bulk and to minimize the homo-molecular IBU contacts seems to be consistent with the relatively largest solubility of IBU therein.

Next, it is also important to analyze the hetero-molecular interactions API–PLA that form in bulk mixtures. Within a PLA molecule, the most active acceptors of hydrogen bonds are presumably its carbonyl oxygen atoms. Fig. 6 depicts RDF corresponding to the O/N–H_API⋯O_PLA hydrogen bonds in the given mixtures at 300 K. Clearly, large amplitudes of the first RDF peak that are significantly larger than one indicate important interactions. These occur in PLA mixtures with NAP and IBU, as amplitudes of the respective RDF peaks exceed values of three or four, respectively, already at very low PLA molar fractions. Interactions of PLA with IND augment higher PLA concentrations, whereas CBZ–PLA hydrogen bonding remains relatively unimportant regardless of the PLA concentrations.

Fig. 6 Radial distribution functions for contacts of the API carboxyl/amide hydrogen atom (H_O/H_N) with the PLA carbonyl oxygen atom (O_C) in API–PLA mixtures simulated at 300 K. Depicted molecular fragments illustrate the particular close contacts (with API being the bottom fragment).

For any of the considered systems, there is no such specific O/N–H_API⋯O_PLA contacts beyond the first coordination shell, as the respective RDF signal only gradually monotonously converges to one after it reaches its minimum around 2.7 Å. Note that the doublet character of the RDF for the CBZ systems corresponds to the stoichiometry of the NH₂ group where the mutual distance of the considered hydrogen atoms is always constrained.

It is interesting to compare the ratios of the first hetero-molecular RDF peak amplitudes at 500 K and at 300 K, which gives an account of the thermal decay of hydrogen bonding in the material.⁵⁰ These ratios range here from 0.42 (for IND–PLA) over 0.51 (NAP–PLA) and 0.55 (IBU–PLA) to 0.63 (CBZ–PLA), which are all significantly lower values than the analogous ratios observed for API–API interactions. It is, however, only ibuprofen for which the given RDF peak exceeds the value of two at 500 K, indicating that the specific H^API_O/N⋯O^PLA_E contacts prevail to a relevant extent in the IBU–PLA mixture even at 500 K. Respective RDFs for H^API_O/N⋯O^PLA_E contacts at 500 K are shown in Fig. S9 (ESI†). The vulnerability of such API–PLA hydrogen bonds to temperature-induced decay suggests that the interaction energies associated with these bonds are lower in absolute value than those of the API–API hydrogen bonds. Results on RDFs for H^API_O/N⋯O^PLA_E contacts with the ester oxygen atom from the polymer indicate that these interactions are even less important as depicted in Fig. S10 (ESI†). This fact naturally restrains the affinity of all hydrogen-bond donors to PEG containing only the ether oxygen atoms. Still, Fig. S11 (ESI†) shows that the simulated structural impact of this IBU-POLY hydrogen bonding is similar for both PEG and PLA.

Equimolar mixtures of IBU with the three target polymers exhibit very similar coordination numbers of key oxygen atoms from the polymer chain that surround a reference IBU hydroxyl hydrogen atom. These values vary between 1.00 and 1.01 indicating that the spatial packing of the bulk mixture and the high numbers of oxygen atoms within the polymer chains always enforce a H_O⋯O_C/E coordination in the bulk. The sole information on the coordination numbers without information about the energetics of those interactions does not suffice to rank the mutual affinity of an API with individual excipients.

To support the discussion on how the hydrogen bonds between the API molecules and polymer chains wane upon heating, Fig. 7 compares coordination numbers that denote how many PLA carbonyl oxygen atoms can be found on average around a hydrogen atom (amide from CBZ or hydroxyl elsewhere) in the first solvation shell. For mixtures containing IBU and NAP, these hetero-molecular contacts are predicted to be significantly more frequent in the first coordination shell than the homo-molecular API–API hydrogen bonds in the bulk. Especially for IBU, the coordination numbers exceeding the value of 1.0 indicate that more than a single carbonyl oxygen atom from PLA can be in the vicinity of IBU hydroxyl moiety. Both homo-molecular and hetero-molecular hydrogen bonds occur to a similar extent in systems containing CBZ, whereas the hetero-molecular bonds are significantly less common in systems with IND despite the general high concentration of the PLA carbonyl oxygen atoms in all systems. For comparison, analogous coordination numbers concerning the hydrogen bonds of API molecules to the ester oxygen atoms of PLA are depicted in Fig. S12 (ESI†).

Fig. 7 Coordination numbers for intermolecular contacts of the API carboxyl/amide hydrogen atom (H_O/H_N) with the PLA carbonyl oxygen atom (O_C) in API–PLA mixtures simulated at 300 K (full columns) and 500 K (empty columns).

Importantly, trends of the calculated API–PLA coordination numbers depicted in Fig. 7 match the hierarchy of experimental solubility of the target API in PLA based excipients. The most soluble API (IBU) forms the highest number of hydrogen bonds to PLA, whereas the lowest soluble IND forms the fewest hydrogen bonds with that polymer.

For the sake of completeness, we also identified the most frequent contacts between two polymer moieties, which proved to occur between the carbonyl oxygen (when present) and hydrogen atom bound to the α-carbon atom next to the carboxyl group. Addition of API molecules into the PLA bulk does not alter the shape of the PLA–PLA RDF and the position of its first peak as shown in Fig. S13 (ESI†). Variation of the RDF amplitudes corresponds then only to changes in the API concentration. Fig. S14 (ESI†) depicts that RDF corresponding to these contacts are sharp for PLA and PVP systems where the carbonyl group is present, but relatively uniform for PEG (where only the ether oxygen could be considered).

Apart from the contact distances, it is relevant to analyze also the spatial distribution of coordinating molecules in the first shell around a reference API molecule. Fig. 8 depicts spatial distribution functions (SDF) simulated at 300 K for pure amorphous IBU and its mixtures with PLA at x_API = 0.85. Notably, isosurfaces corresponding to the local particle density that is 30 times larger than the bulk average are depicted in Fig. 8. This analysis shows that the IBU–IBU hydrogen bonding is more localized in space around the reference carboxyl moiety in pure IBU than in the IBU–PLA mixture as the isosurfaces for the former case cover smaller areas. The closest carboxyl dimers are not that common in pure IBU as the isosurface for the O_C distribution spans the region typically rather than for IBU catemer clusters. Lowering the SDF isovalue from 30 in pure IBU would enlarge the O_C isosurfaces which would then cover also the regions typical for the closest carboxyl dimer, which is not completely ruled out in bulk IBU by the above discussed shapes of the RDF for H_O⋯H_O contacts that are depicted in Fig. S8 (ESI†).

Fig. 8 Spatial distribution functions around a fixed carboxyl frame of IBU molecules simulated at 300 K. Image legend: (a) pure amorphous IBU – dark red, light red, and white surfaces depict distribution of O_C, O_H, and H_O atoms from neighboring IBU molecules, respectively, all isovalues equal to 30; (b) IBU–PLA mixture with x_API = 0.85 – same coloring and atom types as in (a); (c) the same IBU–PLA mixture – dark red and light red surfaces depict the distribution of O_C and O_E atoms from neighboring PLA molecules, respectively, all isovalues equal to 3; (d) the same IBU–PLA mixture – dark red and light red surfaces depict distribution of O_C atoms from neighboring PLA and IBU molecules, respectively, all isovalues equal to 10.

Notably, imprints of both the IBU–IBU chains and closest dimers can be found in the shapes of the IBU⋯IBU contacts in IBU–PLA mixture, which is generally less localized. This is a logical consequence of the dilution by the polymer. SDF shapes for the contacts of IBU with the carbonyl O_C atoms in PLA indicate that the same spatial regions around a reference IBU carboxyl are occupied by both neighboring IBU and PLA molecules which thus inevitably compete for the hydrogen bonding with other IBU molecules. These SDFs also confirm that the ester O_E atoms from PLA are significantly less active in approaching the IBU carboxyl moieties. Finally, Fig. 8 shows that the heteromolecular IBU⋯PLA hydrogen bonds in the mixture are significantly more spatially restricted, whereas the homomolecular IBU⋯IBU hydrogen bonds in the same system can exhibit more spatial orientations from the reference IBU carboxyl. Small API molecules can fit more easily into interstitial spaces and can possibly form important interactions with a larger directional flexibility. Very similar observations can be made from SDF simulated for the remaining target API which are depicted in Fig. S15–S17 (ESI†). The only exception is that those SDFs admit formation of the closest carboxyl or amide dimers in pure amorphous API even at high isosurface values for those materials.

The cohesion of a bulk material is not affected only by the number of individual interactions, but also by their strength. First, lengths of the respective non-covalent O⋯H or N⋯H contacts extracted from MD simulations of the bulk can be analyzed as a tentative indicator of the strength of any site-specific interactions. A balance between the attractive components of a hydrogen bond with the exchange repulsion leads to an optimum length of the hydrogen bond. One can thus assume that the stronger the attractive components, the closer the O⋯H atoms can approach a hydrogen bond. Fig. 9 shows that the most probable observed non-covalent O⋯H contact distances span an interval from 1.5 to 2.1 Å. These contact distances are appreciably shorter (thus assumedly more intense due to stronger attraction) for homo-molecular hydrogen bonds of CBZ, IND and NAP molecules than the lengths of hydrogen bonds of these molecules to the PLA chain. Contact distances of API–API and API–PLA interaction types are comparable in systems containing IBU.

Fig. 9 The most probable contact distances due to hydrogen bonding API⋯API and API⋯PLA in bulk amorphous mixtures at 300 K and x_API = 0.85 obtained from MD simulations. Left – homo-molecular API⋯API contacts; right – hetero-molecular API⋯PLA contacts.

Furthermore, one can also compare the radii of the first solvation shells, which correspond to the distance coordinate of the first minimum on the RDF signals. While these radii range from 2.3 to 3.0 Å for the homo-molecular hydrogen bonds at 300 K, the radii for the hetero-molecular interactions amount to 2.4–3.3 Å for the H^API_O/N⋯O^PLA_C contacts or even to 4.1–4.6 Å for the H^API_O/N⋯O^PLA_E contacts. The latter can be then hardly taken for significant hydrogen bonds due to those too large interatomic separations.

To deepen the analysis of the key pair interaction for the cohesion of the target dispersions, we compared the most probable contact distances due to hydrogen bonding in bulk amorphous mixtures at 300 K and x_API = 0.85 extracted from MD simulations with their analogues obtained by quantum-chemical optimizations of isolated clusters of a single API molecule with a chain fragment representing a methyl-capped monomer or a tetramer.

The results of this quantum-chemical analysis are listed in Table 2. Since MD provides results for a densely packed bulk phase at a finite temperature, thermal effects and steric constraints are expected to propagate to somewhat longer contact distances than those obtained from a static DFT calculation for an isolated B3LYP-optimized cluster of an API with a monomer or oligomer without any vibrational effects. This trend is indeed visible for CBZ–PLA, NAP–PLA and IBU–PVP interactions (concerning carbonyl oxygen atoms in polymers). Since the current MD predicts shorter contacts due to IND–PLA interactions in the bulk than what the DFT-optimized clusters yield, one can conclude that the current force field overestimates these attractive IND–PLA interactions. Note that SAPT typically predicts an optimum contact of a hydrogen bond at 1.8 to 2.0 Å.⁸²

Table 2 Calculated contact distances (d, in Å) corresponding to hydrogen bonds between the target API molecules (H_O atom) with individual polymers (O_C and O_E atoms), corresponding interaction energies (in kJ mol⁻¹) and signed electron density (ρ_EL, in atomic units) coordinate of the reduced density gradient that is important for the classification of non-covalent interactions

System	d ^RDFMD ^a	d ^monoDFT ^b	d ^oligoDFT ^c	ε SAPT	ε DL	sign(λ2) ρEL^f
a Position of the first RDF peak due to HO/N⋯O contacts based on MD simulations of bulk amorphous mixtures at 300 K. b Atomic distance of the HO/N⋯O contact in an isolated cluster of API and chain monomer optimized using B3LYP-D3(BJ)/6-311+G(d,p). c Distance of HO/N⋯O atoms in an isolated cluster of API and chain oligomer optimized using B3LYP-D3(BJ)/6-311+G(d,p). d Pair interaction energy within an isolated cluster of API and chain monomer calculated using sSAPT0/jun-cc-pVDZ. e Dispersion-less interaction energy within an isolated cluster of API and chain monomer calculated using sSAPT0/jun-cc-pVDZ as: εDL= εSAPT – εDISP according to the definition in ref. 82. f Output of the NCI analysis^63 – electron density in atomic units multiplied with the sign of the second eigenvalue of the electron-density Hessian matrix, which corresponds to a critical point of the reduced-density gradient due to the API–polymer hydrogen bond.
Polymer site	HO⋯OC contact
PLA–CBZ	2.09	1.96	2.38	–70.9	–17.7	−0.019
PLA–IBU	1.77	1.77	1.76	–64.5	–20.3	−0.035
PLA–IND	1.50	1.78	1.88	–51.2	–16.6	−0.027
PLA–NAP	1.81	1.78	1.81	–51.5	–12.2	−0.032
PVP–IBU	2.09	1.96	2.38	–69.3	–23.6	−0.050
Polymer site	HO⋯OE contact
PLA–CBZ	5.31	2.11	2.06	–27.0	0.0	−0.020
PLA–IBU	3.90	1.85	1.97	–39.3	–9.8	−0.021
PLA–IND	3.74	1.88	2.12	–41.2	–3.5	−0.020
PLA–NAP	3.92	1.80	1.89	–34.9	–16.9	−0.028
PEG–IBU	5.31	2.11	2.06	–48.3	–16.0	−0.029

---

The predicted too close contact distances of hydrogen bonded IND molecules are probably a computational artifact of the underlying force field that overestimates the electrostatic attractive and/or underestimates the repulsive components of this atomic contact. Among the force-field models for the target API, atomic charges of the hydroxyl hydrogen atoms are relatively similar (ranging from 0.42 to 0.46 q_e) but the amino hydrogen atoms in CBZ are outlying with an atomic charge of 0.31 q_e. The most active acceptor of hydrogen bonds in the API molecules, being their carbonyl oxygen atom from their carboxyl/amide moiety, possesses an atomic charge ranging from −0.44 to −0.60 q_e (IND molecule contains the most negative oxygen). This particular parametrization contributes to the trends of hydrogen hierarchy modeled for individual API systems.

DFT-based geometry optimizations for API clusters with longer oligomer chains were supposed to partly mimic the steric constraints that may impede the hydrogen-bonding contact distance in bulk. These distances were predominantly longer (by up to 0.40 Å) than those in the clusters with the chain monomers. In addition, DFT calculations admit that the hydrogen bonding to the ether/ester oxygen atoms in polymer chains should be relatively strong, yet somewhat weaker (longer by 0.1–0.2 Å) than those related to the carbonyl oxygen atoms in the optimum isolated clusters.

Despite this relatively small DFT-predicted difference in the hydrogen-bond lengths, MD predicts that the H^API_O/N⋯O^PLA_E contact distances are larger by a factor of 2–3 in the bulk, indicating the importance of the competition of individual molecular interactions in densely packed materials. Unlike in a small isolated cluster, only the most beneficial interactions (being H^API_O/N⋯O^PLA_C) prevail in a densely packed bulk environment. Still, given this difference in the assessment of the H^API_O/N⋯O^PLA_E contacts by MD and DFT, one can assume that the current MD model overstates the API⋯O^POLY_C hydrogen bonds at the expense of their API⋯O^POLY_E competitors.

SAPT calculations offer another metric of the hydrogen-bonding strength that is based on analysis of the dispersion-less interaction energies as listed in Table 2. These sums of electrostatic, exchange and induction contributions suggest that the IBU–PLA hydrogen bonds are the strongest among individual API–PLA systems, and the IBU–PVP hydrogen bonds are the strongest among the considered IBU–polymer systems.

In complex molecular clusters or in bulk environments, however, where too many atoms get in close contacts, calculations of total pair interactions, or even their purely electrostatic or dispersion-less components can be misleading to assess the strength of a hydrogen bond.⁸² On the basis of the quantum-chemical analysis of non-covalent interactions (NCIs), molecular sites corresponding to the most intense interactions between the API and polymer molecules, such as hydrogen bonding, can be identified and their strength compared among various conformers even when there are numerous other degrees of freedom that can blur comparison of total interaction energies of large molecules.⁶³ Fig. S18–S21 (ESI†) depict the so-called RDG plots for clusters of PLA oligomers with individual target API. Critical RDG points (sharp spikes approaching zero RDG values) are mapped in NCI plots to isosurfaces in particular intermolecular regions, also comparing the intensities of hydrogen bonds accepted by the carbonyl or ester oxygen atoms in PLA chains. Signed electron-density coordinates of the relevant RDG critical points are listed in Table 2. It is obvious that the hydrogen bonds from the API hydroxyl groups to the PLA carbonyl oxygen atom are appreciably stronger in such isolated optimum clusters as the RDG critical points lie at more negative signed electron densities (far more to the left in the RDG plots in Fig. S18–S21, ESI†) than to what occurs for the hydrogen bonds to the ester oxygen atom. On the other hand, both hydrogen-bond oxygen acceptors in PLA are predicted to behave very similarly when the amide moiety of CBZ donates the hydrogen bond, which is, however, predicted to be weak. Comparing these positions of RDG critical points for individual API, the hydrogen bonding strength with PLA increases in a row: CBZ < IND < NAP < IBU, which again agrees with the experimental hierarchy of API solubilities in PLA based polymers.

These findings from NCI are in agreement with the shapes of simulated RDF and corresponding coordination numbers, indicating that the interactions in bulk liquid and an isolated cluster are very similar. Furthermore, one can locate vast isosurfaces in the NCI plots that correspond to attractive dispersion interactions between the aromatic cores of the target API molecules with the planar carbonyl moieties of the PLA chain. Current calculations thus suggest that polymer carriers containing some segments capable of π–π stacking in their structure can attain more beneficial spatial packing and dispersion interactions with the dispersed API if that contains some aromatic moieties.

Fig. S19, S22 and S23 (ESI†) then compare interactions of IBU molecules with individual considered polymers. The hydrogen-bonding strength of IBU to the polymers increases in the row: PEG < PLA < PVP. Since PEG contains no carbonyl groups, it can accept only weaker hydrogen bonds with its ether oxygen atoms. Carbonyl moiety from PVP is then predicted to accept very strong hydrogen bonds from IBU with the corresponding RDG critical point located below −0.050 a.u., being an imprint of the presence of the electron-rich side-chain pyrrolidone moieties in its structure.

Such findings about the relative strength of the hydrogen bonding between API and PLA molecules are in agreement with the observed positive excess enthalpies of the mixing for these systems, which can be interpreted as if the former strong API–API interactions are partially replaced by somewhat weaker API–PLA interactions which is unfavorable in terms of enthalpy.

Notably, experimentally determined the solubility of the target API within PLA based polymers at 298.15 K matches this trend of the modelled hydrogen-bonding strength, as well as the trend of solubility of IBU in various polymers that were mentioned in the Introduction section.

3.3. Internal dynamics of mixtures

The ability of molecules to diffuse at the microscopic level can affect the macroscopic phase behavior of materials. In the current context, more mobile molecules may impart a facile route towards the spurious recrystallization of amorphous solid dispersions. Molecular dynamics enables us to analyze the time-dependence of mean-square displacements (MSD) of molecules and to derive their diffusivities in various environments via the Einstein theory. First, there are significant differences in self-diffusivities of the API molecules in their neat amorphous state at 500 K when all systems are convincingly liquid. Fig. S24 (ESI†) depicts that the Einstein diffusive regime with a steady linear increase of the MSD over time is established.⁸³ These differences among the diffusivities remain visible in the simulated data even at 300 K, although the behavior of such MSD data corresponds rather to the glassy state or to a sluggish liquid due to the very low MSD values reached after 10 ns and frequent difficulties with reaching the steady linear diffusive regime as depicted in Fig. S25 (ESI†).

Converting the MSD into diffusivities of individual API, Fig. 10 reveals that IND molecules are the least mobile in their neat bulk phase, which agrees with the above explored strongest hydrogen bonding between its molecules in MD simulations which, together with the largest molecular size among the target API, possibly imposes hindrance to diffusion of hydrogen bonded IND associates through the bulk. Interestingly, IBU molecules are more mobile by nearly an order of magnitude than the remaining target API, which correlates well with the lowest experimental melting and glass-transition temperatures among the considered API.⁸⁴

Fig. 10 Diffusivities of API molecules as functions of time in various bulk systems simulated at 500 K.

Fig. 10 also compares how the diffusivity of the API molecules varies in the mixtures with PLA within the considered concentration range. The most significant simulated feature is the massively attenuated diffusivity of IBU molecules in their mixtures with PLA. The respective MSD is lower in the mixture by a factor of up to six (regardless of API concertation in this range) than what was found for the pure API at both 300 K and 500 K. To explain this behavior, one should consider first that IBU molecules were shown to be very mobile in their neat phase, and concurrently, the diffusivity of neat PLA was simulated to be two orders of magnitude lower as depicted in Fig. 10. Furthermore, the existence of site-specific PLA–IBU interactions that bind the small API molecules to the long chain, which does not usually exhibit any significant translation movement as a whole molecule, can contribute to an explanation of this sluggish behavior of IBU molecules in its mixtures with PLA. The observed decrease of IBU diffusivity upon lowering its concentration in the mixture seems to be in agreement with this interpretation.

In this context, Fig. 10 also confirms that D values for PLA remain very low in the mixtures with all considered API, with a weak decreasing tendency as API concentration rises. Fig. S26–S29 (ESI†) then compare the MSD of the centers of mass of PLA chains with MSD evaluated only for the chain termini. Especially at 300 K, the end groups are somewhat more mobile in the bulk due to the conformational flexibility of polymer chains, but still, the chain termini remain significantly less mobile in the bulk mixtures than the API molecules themselves.

Both neat CBZ and NAP are predicted to exhibit diffusivities that are roughly five-fold lower than IBU, but still two orders of magnitude larger than those of neat PLA, as shown in Fig. 10. Diffusivity of both APIs is then predicted to drop down only very little in the mixtures as their concentration decreases. While CBZ was shown not to interact with PLA any importantly, the observed slight decrease of its diffusivity can be attributed to a mechanistic hindrance imposed by the presence of the large polymer molecules.

In contrast, diffusion of IND molecules is predicted to be appreciably augmented by the presence of PLA at elevated temperatures, but IND diffusivity is rather indifferent to the PLA's presence at ambient temperatures. To explain this behavior, one should consider that current simulations yield similar (very low) diffusivity values for both neat PLA and IND. Elevated temperature promotes large conformational variability of the PLA chain which augments at least the mobility of individual chain segments if not the whole large polymer molecules. Fig. S26–S29 (ESI†) show that PLA termini reach similar diffusivities to the API molecules in their mixtures at 300 K. One can thus assume that the IND molecules that are dispersed around PLA chains simply follow the conformational chaos of those chains, which in turn speeds up their diffusion in the bulk. Notably, IND diffusivity is here predicted to increase as the PLA concentration rises, which supports the hypothesis that the polymer dynamics affect the IND molecules.

Since the conformational flexibility of the polymers can significantly contribute to the transport behavior in the bulk, we analyzed whether the simulated polymer chains can visit various conformation states in the bulk mixture with IBU. We selected the instantaneous end-to-end distance L of individual chains as an easily accessible metric to be derived from the MD simulations. We also assume that L is correlated with the effective radius of the polymer molecule. Conformations of the initially pre-folded globules constituting the bulk were let to evolve independently in the simulation, as much as the environment and the conformational preferences of the polymer enabled.

For our production trajectories of the mixtures at 300 K and at 500 K, we constructed histograms of the instantaneous L values of the three polymer chains that are depicted in Fig. S32 (ESI†). At both temperatures, the observed distribution of the end-to-end distances is relatively broad, indicating that individual chains have high chances of visiting conformations that are largely distinct from the initial pre-folded chain structures, including both more unfolded or collapsed states (the latter not for PVP). Obviously, the systems are more fluid at the higher temperature, resulting also in a better sampling of the phase space, which in turn propagates to smoother histograms of the instantaneous end-to-end distances of individual chains over the production MD simulations. Importantly, the end-to-end distances are longer at the higher temperature, which indicates a larger extent of chain unfolding at 500 K. Notably, a principal constraint limiting the chain unfolding to some extent is the finite size of the simulation box, which is closely connected with the numbers of atoms included in the box. We observed this finite-size effect in the given simulations in our recent paper,³⁴ having found, however, that densities or enthalpies of the bulk are insensitive to the conformation-constraints imposed due to simulation-box size limited to the current 60 Å or larger size.

Comparing the diffusivities of IBU in its mixtures with individual polymers, it increases in a row: PVP ≪ PLA < PEG (denoting the other component) at 500 K as depicted in Fig. S27, S30 and S31 (ESI†). This ranking is governed by the very slow diffusivity of the PVP in the current simulations, which anchor also the small IBU molecules in the bulk mixture to a large extent. Under ambient conditions, the MSD of IBU across individual polymers does not differ that largely, but the PVP system is still the least diffusive one.

Finally, the simulated diffusivities of API molecules in the bulk can contribute to the interpretation of the extent of the kinetic stabilization of super-saturated API–POLY mixtures. The trend of rates of the observed phase separation of super-saturated mixtures that increase in a row: IND < NAP < IBU²³ matches the hierarchy of the simulated diffusivities of API molecules in the bulk systems. The larger mobility of small API molecules in the mixture naturally facilitates any phase transitions and structural rearrangements of super-saturated systems.

The current hierarchy of diffusivities at 500 K does not, nevertheless, explain the varying glass forming tendencies of the target API. Both the API with the highest and lowest simulated diffusivities belong to very potent glass formers (class III according to the established classification⁸⁵). On the other hand, it is very difficult to prepare an amorphous form of NAP (class I⁸⁵), the diffusivity of which is predicted to be between those of IBU and IND.

3.4. Glass transition temperatures

Computational analysis of T_g represents a link between the microscopic interpretation and macroscopic behavior of the target systems. Illustrative data sets that were also used for determination of T_g with the hyperbola fits are displayed in Fig. 9 for each API–PLA system.

Note that the hyperbola fit of the volumetric data from non-equilibrium MD was employed to minimize any arbitrariness due to human intervention within analyses of the density vs. temperature data.¹⁶ For systems with a clearly visible ρ(T) trend shift, which corresponds to a sudden glass transition, the determination of T_g is more reliable and has a lower uncertainty.¹⁸ If the transition between the liquid and the glassy state appears to be rather gradual, computational determination of T_g becomes more sensitive to any noise and the associated uncertainty of the resulting value is higher. It becomes important to repeat the simulated cooling runs multiple times especially for such stiff systems to rationalize the computational uncertainty of the predicted T_g.³⁴ While IBU systems represent the first case with a sudden trend shift and low computational uncertainties, IND systems are representatives of the latter scenario exhibiting a rather gradual ρ(T) trend shift. As a metric of how the trend shift obvious is, one can perform the hyperbola fit with its c parameter constrained to c = −∞ (the result is a broken line) and with a fully flexible c (resulting in a smooth curve). Proximity of the T_g results from these two fits of the same raw density data correlates strongly with the abruptness of the trend shift, as can be seen in Fig. 11. Uncertainty of the calculated T_g is then naturally lower when the trend shift is well localizable.

Fig. 11 Illustration of the computational localization of the glass transition temperature T_g of API–PLA systems with x_API = 0.85 based on non-equilibrium cooling of the melt and hyperbola fit of the simulated density vs. temperature data. Two vertical dotted lines correspond to different values of the c parameter in the hyperbola equation: T^∞_g value corresponds to c = −∞; T^H_g value corresponds to adjustable c.

Uncertainties of individual T_g data points corresponding to a single cooling run were estimated from synthetic perturbation of the simulated ρ(T) data set. Within a single trajectory, the root mean squared residuals (RMSR) value for the hyperbola fit of the raw simulated density data was evaluated at first. Subsequently, each of the simulated density entries at an instantaneous time step during the cooling was augmented by the RMSR value multiplied by a factor, which was randomly selected for each entry from −1, 0 or 1 value. This procedure was repeated 40 times independently. Finally, the standard deviation of T_g values determined from these perturbed ρ(T) data set was taken as the uncertainty of T_g based on a single ith cooling run, denoted σ_i. Since the cooling was repeated in five independent replicas, another uncertainty term, σ_R, related to the reproducibility of the cooling was accounted for. This was estimated from the standard deviation of T_g among the individual cooling runs. Finally, the uncertainty of the resulting T_g value combined both terms according to the equation:

Results on T_g determined from individual cooling runs are summarized in Table S2 (ESI†). T_g values averaged over the individual independent runs for each system are depicted in Fig. 12. Simulation results on T_g for neat PLA were taken from our earlier work also using all-atom MD simulations of non-equilibrium cooling and the hyperbola fit.³⁴ Note that the simulated T_g values for pure API were taken from another our work, which used the same force field as in this work, but T_g was determined there only using the equilibrium graphical intersection method,¹⁵ the results of which are usually burdened by somewhat larger uncertainties.³⁴

Fig. 12 Calculated glass transition temperatures T_g of API–PLA mixtures with composition x_API = 0.85 and their standard uncertainties with results obtained earlier for pure API¹⁵ and pure PLA.³⁴ Lines connecting the explicitly calculated data point serve only to guide the eye.

As long as T_g calculated for the neat polymer anchors data sets for all considered API at the same temperature for x_API = 0, T_g data for pure API play a crucial role in affecting the resulting trends of T_g with respect to the mixture composition. Fig. 12 depicts that this dependence is likely monotonous, meaning that T_g for the mixture ranges in between the values for pure compounds, but is appreciably non-linear in most cases. Such results on the T_g trends are in qualitative agreement with the experimental findings²³ for similar mixtures and also with simulation results performed earlier for mixtures of the target API with low-molecular excipients.¹⁵

Obviously, it generally holds that the higher the diffusivity of a material is, the faster molecular motion occurs, which leads in turn to a lower T_g value for that material. In this context, current simulations indicate that IBU and NAP, both exhibiting appreciably lower T_g than neat PLA, experience a large T_g increase upon mixing with PLA. This behavior is consistent with the reported important API–PLA interactions and the decrease of MSD of both APIs upon mixing with the polymer carrier. For mixtures with a visible decrease of diffusivity of the API upon mixing with the polymer carrier, one should anticipate a higher T_g value of the mixture than what the neat API exhibits. Pure CBZ was predicted to exhibit a similar T_g as PLA does, and concurrently, the hetero-molecular interactions in this system were shown to be rather weak so there is no reason for an outlying behavior of T_g of the simulated CBZ–PLA mixtures that exhibit a relatively lesser non-ideality.

Statistical uncertainties associated with T_g for the IND containing systems are somewhat larger due to the gradual character of the glass transition within the simulations. As T_g for pure IND is predicted to be higher than for neat PLA, mixing of these compounds yields a lower T_g than what pure IND exhibits. Hypothetically, should a multi-component system exhibit higher T_g than all of its components in their pure state, specific strong hetero-molecular interactions would be the cause. However, such a behavior was not observed among the considered systems although certain interactions between the API and the polymer carrier were identified in the simulations.

Recent computational studies have predominantly reported overestimated T_g results (by up to higher tens of Kelvin) for various molecular materials.^{15,17–19,34} This overestimating trend can be visible also in this work when actual T_g computed for individual pure compounds are compared with experimental data, as listed in Table S3 (ESI†). One of the reasons may be the too fast cooling rate imposed in the simulations which is orders-of-magnitude larger than that applied in real experiments.⁴³ On the other hand, variation of the cooling rate imposed on simulated polymer samples did not alter the obtained T_g results in a systematic way.³⁴ Real computational uncertainties of our T_g (related to the difference between simulation and experiment) for the mixtures should be expected to be somewhat higher than their statistical sampling uncertainties (related to data processing at a given level of theory).

4. Conclusions

Molecular-dynamics simulations with all-atom resolution and quantum-chemical calculations were performed to investigate the properties of amorphous dispersions of four archetypal active pharmaceutical ingredients with three biocompatible polymers, including their six binary mixtures. This effort provided an insight on the inner structure and dynamics at the microscopic scale and of the energetics of the mixing and the vitrification tendency at the macroscopic level. Identification of individual molecular interactions between the active pharmaceutical ingredient and the polymer chain enabled us to interpret experimentally phenomena that had been observed earlier.

Computational conclusions about the potential of individual APIs to donate hydrogen bonds increasing in a row: IND < NAP < IBU match the real solubility trends of these APIs in PLA based excipients. CBZ is modeled to interact with PLA even weaker than IND does, but we are not aware of any experimental evidence to verify this computational conclusion. Computationally ranked potential of polymers to accept hydrogen bonding, increasing in a row: PEG < PLA < PVP also agrees with the appreciably larger experimental solubility of IBU in PVP.

The larger miscibility of ibuprofen with derivatives of polylactic acid can be attributed to the strong hydrogen bonding of the API to the polymer where the ibuprofen molecule acts as the donor and PLA as the acceptor.

Trends of the simulated diffusivities of individual API molecules within the mixtures agree with the experimental observations of the kinetic stability of super-saturated mixtures, which decreases in a row: IND > NAP > IBU. The high kinetic stability of the dispersion of indomethacin with derivatives of polylactic acid can be attributed to the very slow dynamics of the indomethacin molecules themselves in the bulk which are orders of magnitude than those for the remaining API molecules. Even though the polymer appears to be more diffusive at elevated temperatures than IND, the latter retains very low diffusivities even in the mixtures.

Polymer chains impose appreciable constraints on the diffusion of ibuprofen molecules in the mixture, which imparts some kinetic stabilization with respect to recrystallization of ibuprofen dispersed in the polymer when compared to pure amorphous ibuprofen. Nevertheless, a comparison with the experimental evidence reveals that this mobility hindrance is not sufficient to prevent recrystallization from massively oversaturated dispersions.

Presented analyses illustrate how molecular simulations can contribute to interpretation of kinetic and thermodynamic phenomena that are crucial to understanding upon design of novel pharmaceutic formulations. Current computational results offer a consistent description of structural, energetic and dynamic aspects that are compared for pure compounds and their mixtures. Incorporating explicit polarizability models into such simulations in the future is expected to further improve the quality of the predicted internal dynamics of the bulk systems and lower the associated computational uncertainties.

Author contributions

Vladislav Aulich – investigation, formal analysis, writing – original draft. Jan Ludík – investigation, formal analysis, visualization. Michal Fulem – conceptualization, funding acquisition, project administration, writing – review & editing. Ctirad Červinka – conceptualization, methodology, resources, supervision, data curation, writing – review & editing.

Data availability

The data supporting this article have been included as part of the ESI.†

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

The authors acknowledge financial support from the Czech Science Foundation under the project 22-07164S. Computational resources were provided by the e-INFRA CZ project (ID:90254), supported by the Ministry of Education, Youth and Sports of the Czech Republic.

References

1. V. Caron, L. Tajber, O. I. Corrigan and A. M. Healy, Mol. Pharmaceutics, 2011, 8, 532–542.
2. J. Brouwers, M. E. Brewster and P. Augustijns, J. Pharm. Sci., 2009, 98, 2549–2572.
3. C. Leuner and J. Dressman, Eur. J. Pharm. Biopharm., 2000, 50, 47–60.
4. L.-F. Huang and W.-Q. Tong, Adv. Drug Delivery Rev., 2004, 56, 321–334.
5. E. Batisai, ChemistryOpen, 2021, 10, 1260–1268.
6. T. Vasconcelos, B. Sarmento and P. Costa, Drug Discovery Today, 2007, 12, 1068–1075.
7. A. Mathers, M. Pechar, F. Hassouna and M. Fulem, Int. J. Pharm., 2022, 623, 121855.
8. A. Newman and G. Zografi, Mol. Pharmaceutics, 2022, 19, 378–391.
9. A. Prudic, Y. Ji and G. Sadowski, Mol. Pharmaceutics, 2014, 11, 2294–2304.
10. P. Srinarong, H. de Waard, H. W. Frijlink and W. L. J. Hinrichs, Expert Opin. Drug Delivery, 2011, 8, 1121–1140.
11. K. Sekiguchi, N. Obi and Y. Ueda, Chem. Pharm. Bull., 1964, 12, 134–144.
12. C. Luebbert, F. Huxoll and G. Sadowski, Molecules, 2017, 22, 296.
13. A. Mathers, M. Pechar, F. Hassouna and M. Fulem, Int. J. Pharm., 2023, 648, 123604.
14. J. Pavliš, A. Mathers, M. Fulem and M. Klajmon, Mol. Pharmaceutics, 2023, 20, 3960–3974.
15. C. Červinka and M. Fulem, Pharmaceutics, 2021, 13, 1253.
16. L. Alzate-Vargas, M. E. Fortunato, B. Haley, C. Li, C. M. Colina and A. Strachan, Modell. Simul. Mater. Sci. Eng., 2018, 26, 065007.
17. F. J. Carmona Esteva, Y. Zhang, Y. J. Colón and E. J. Maginn, J. Phys. Chem. B, 2023, 127, 4623–4632.
18. F. J. Carmona Esteva, Y. Zhang, E. J. Maginn and Y. J. Colón, J. Chem. Phys., 2024, 161, 014108.
19. M. Klajmon and C. Červinka, J. Phys. Chem. B, 2022, 126, 2005–2013.
20. R. Schneider, L. Taspinar, Y. Ji and G. Sadowski, Int. J. Pharm., 2020, 582, 119317.
21. C. Luebbert, F. Huxoll and G. Sadowski, Molecules, 2017, 22, 296.
22. A. Prudic, A.-K. Lesniak, Y. Ji and G. Sadowski, Eur. J. Pharm. Biopharm., 2015, 93, 88–94.
23. A. Iemtsev, F. Hassouna, M. Klajmon, A. Mathers and M. Fulem, Eur. J. Pharm. Biopharm., 2022, 179, 232–245.
24. K. Lehmkemper, S. O. Kyeremateng, O. Heinzerling, M. Degenhardt and G. Sadowski, Mol. Pharmaceutics, 2017, 14, 157–171.
25. P. J. Marsac, T. Li and L. S. Taylor, Pharm. Res., 2009, 26, 139–151.
26. A. Mahieu, J.-F. Willart, E. Dudognon, F. Danède and M. Descamps, Mol. Pharmaceutics, 2013, 10, 560–566.
27. Y. Sun, J. Tao, G. G. Z. Zhang and L. Yu, J. Pharm. Sci., 2010, 99, 4023–4031.
28. M. M. Knopp, N. E. Olesen, P. Holm, P. Langguth, R. Holm and T. Rades, J. Pharm. Sci., 2015, 104, 2905–2912.
29. M. B. Rask, M. M. Knopp, N. E. Olesen, R. Holm and T. Rades, Int. J. Pharm., 2018, 540, 98–105.
30. A. Iemtsev, A. Zemánková, F. Hassouna, A. Mathers, M. Klajmon, M. Slámová, L. Malinová and M. Fulem, Int. J. Pharm., 2022, 613, 121424.
31. A. Mathers and M. Fulem, Int. J. Pharm., 2024, 664, 124613.
32. A. Iemtsev, M. Klajmon, F. Hassouna and M. Fulem, Pharmaceutics, 2023, 15, 645.
33. K. Goloviznina, Z. Gong and A. A. H. Padua, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2022, 12, e1572.
34. M. Klajmon, V. Aulich, J. Ludík and C. Červinka, Ind. Eng. Chem. Res., 2023, 62, 21437–21448.
35. J. H. McAliley and D. A. Bruce, J. Chem. Theory Comput., 2011, 7, 3756–3767.
36. C.-E. Fang, Y.-C. Tsai, C. Scheurer and C.-C. Chiu, Polymers, 2021, 13, 1131.
37. N. S. V. Barbosa, Y. Zhang, E. R. A. Lima, F. W. Tavares and E. J. Maginn, J. Mol. Model., 2017, 23, 194.
38. A. N. Leonard, R. W. Pastor and J. B. Klauda, J. Phys. Chem. B, 2018, 122, 6744–6754.
39. Y. Zhou, W. A. Saidi and K. A. Fichthorn, J. Phys. Chem. C, 2014, 118, 3366–3374.
40. A. Kyrychenko, O. M. Korsun, I. I. Gubin, S. M. Kovalenko and O. N. Kalugin, J. Phys. Chem. C, 2015, 119, 7888–7899.
41. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
42. M. Greiner, E. Elts, J. Schneider, K. Reuter and H. Briesen, J. Cryst. Growth, 2014, 405, 122–130.
43. T.-X. Xiang and B. D. Anderson, Mol. Pharmaceutics, 2013, 10, 102–114.
44. F. Römer and T. Kraska, Mol. Simul., 2012, 38, 152–160.
45. J. N. Canongia Lopes and A. A. H. Pádua, Theor. Chem. Acc., 2012, 131, 1129.
46. V. Štejfa, H. Lyshchuk, K. Babková, M. Krupička, J. Ludík, M. Fulem, C. Červinka and K. Růžička, J. Mol. Liq., 2024, 394, 123804.
47. C. E. S. Bernardes and A. Joseph, J. Phys. Chem. A, 2015, 119(12), 3023–3034.
48. A. J. M. Sweere and J. G. E. M. Fraaije, J. Chem. Theory Comput., 2017, 13, 1911–1923.
49. S. H. Jamali, T. v Westen, O. A. Moultos and T. J. H. Vlugt, J. Chem. Theory Comput., 2018, 14, 6690–6700.
50. V. Pokorný, V. Štejfa, K. Růžička and C. Červinka, Phys. Chem. Chem. Phys., 2021, 23, 26874–26886.
51. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
52. L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, J. Comput. Chem., 2009, 30, 2157–2164.
53. A. A. H. Pádua, K. Goloviznina and Z. Gong, fftool: XML force field files, ENS Lyon, 1.2.1 edn, 2023.
54. M. E. Tuckerman, J. Alejandre, R. López-Rendón, A. L. Jochim and G. J. Martyna, J. Phys. A: Math. Gen., 2006, 39, 5629.
55. R. W. Hockney and J. W. Eastwood, Computer simulation using particles, Taylor&Francis, New York, 1988.
56. J.-P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, J. Comput. Phys., 1977, 23, 327–341.
57. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, 2nd edn, Oxford, 2017.
58. A. Dequidt, Diffraction Code, Universite Clermont Auvergne, Aubière, France, 2016.
59. V. F. Sears, Neutron News, 1992, 3, 26–37.
60. X. Mao, P. Brown, C. Červinka, G. Hazell, H. Li, Y. Ren, D. Chen, R. Atkin, J. Eastoe, I. Grillo, A. A. H. Padua, M. F. Costa Gomes and T. A. Hatton, Nat. Mater., 2019, 18, 1350–1357.
61. P. N. Patrone, A. Dienstfrey, A. R. Browning, S. Tucker and S. Christensen, Polymer, 2016, 87, 246–259.
62. R. F. W. Bader, Chem. Rev., 1991, 91, 893–928.
63. E. R. Johnson, S. Keinan, P. Mori-Sánchez, J. Contreras-García, A. J. Cohen and W. Yang, J. Am. Chem. Soc., 2010, 132, 6498–6506.
64. C. Červinka, M. Klajmon and V. Štejfa, J. Chem. Theory Comput., 2019, 15, 5563–5578.
65. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
66. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
67. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman and G. B. Scalmani, V. Petersson, G. A. Nakatsuji, H. Li, X. Caricato, M. Marenich, A. V. Bloino, J. Janesko, B. G. Gomperts, R. Mennucci, B. Hratchian, H. P. Ortiz, J. V. Izmaylov, A. F. Sonnenberg, J. L. Williams-Young, D. Ding, F. Lipparini, F. Egidi, F. Goings, J. Peng, B. Petrone, A. Henderson, T. Ranasinghe, D. Zakrzewski, V. G. Gao, J. Rega, N. Zheng, G. Liang, W. Hada, M. Ehara, M. Toyota, K. Fukuda, R. Hasegawa, J. Ishida, M. Nakajima, T. Honda, Y. Kitao, O. Nakai, H. Vreven, T. Throssell, K. Montgomery, J. A. Peralta, J. E. Ogliaro, F. Bearpark, M. J. Heyd, J. J. Brothers, E. N. Kudin, K. N. Staroverov, V. N. Keith, T. A. Kobayashi, R. Normand, J. Raghavachari, K. Rendell, A. P. Burant, J. C. Iyengar, S. S. Tomasi, J. Cossi, M. Millam, J. M. Klene, M. Adamo, C. Cammi, R. Ochterski, J. W. Martin, R. L. Morokuma, K. Farkas, O. Foresman, J. B. Fox, Gaussian 16, Revision B.01, Gaussian, Inc., Wallingford, 2016.
68. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
69. T. M. Parker, L. A. Burns, R. M. Parrish, A. G. Ryno and C. D. Sherrill, J. Chem. Phys., 2014, 140, 094106.
70. E. Papajak and D. G. Truhlar, J. Chem. Theory Comput., 2011, 7, 10–18.
71. D. G. A. Smith, L. A. Burns, A. C. Simmonett, R. M. Parrish, M. C. Schieber, R. Galvelis, P. Kraus, H. Kruse, R. Di Remigio, A. Alenaizan, A. M. James, S. Lehtola, J. P. Misiewicz, M. Scheurer, R. A. Shaw, J. B. Schriber, Y. Xie, Z. L. Glick, D. A. Sirianni, J. S. O’Brien, J. M. Waldrop, A. Kumar, E. G. Hohenstein, B. P. Pritchard, B. R. Brooks, H. F. Schaefer, III, A. Y. Sokolov, K. Patkowski, A. E. DePrince, III, U. Bozkaya, R. A. King, F. A. Evangelista, J. M. Turney, T. D. Crawford and C. D. Sherrill, J. Chem. Phys., 2020, 152, 184108.
72. J. Garcia, R. Podeszwa and K. Szalewicz, J. Chem. Phys., 2020, 152, 184109.
73. T. K. Giri, H. R. Badwaik, A. Alexander and D. K. Tripathi, Int. J. Appl. Biol. Pharm. Technol., 2010, 1, 793.
74. A. Minecka, E. Kamińska, K. Jurkiewicz, D. Heczko, B. Hachuła, W. Pisarski, K. Kamiński and M. Paluch, Phys. Chem. Chem. Phys., 2020, 22, 295–305.
75. S. Bates, G. Zografi, D. Engers, K. Morris, K. Crowley and A. Newman, Pharm. Res., 2006, 23, 2333–2349.
76. P. Karmwar, K. Graeser, K. C. Gordon, C. J. Strachan and T. Rades, Int. J. Pharm., 2011, 417, 94–100.
77. S. J. L. Billinge, T. Dykhne, P. Juhás, E. Božin, R. Taylor, A. J. Florence and K. Shankland, CrystEngComm, 2010, 12, 1366–1368.
78. K. Ostrowska, M. Kropidłowska and A. Katrusiak, Cryst. Growth Des., 2015, 15, 1512–1517.
79. P. J. Cox and P. L. Manson, Acta Crystallogr., Sect. E: Struct. Rep. Online, 2003, 59, o986–o988.
80. G.-M. Tang, J.-H. Wang, C. Zhao, Y.-T. Wang, Y.-Z. Cui, F.-Y. Cheng and S. W. Ng, CrystEngComm, 2015, 17, 7258–7261.
81. K. S. Eccles, S. P. Stokes, C. A. Daly, N. M. Barry, S. P. McSweeney, D. J. O'Neill, D. M. Kelly, W. B. Jennings, O. M. Ni Dhubhghaill, H. A. Moynihan, A. R. Maguire and S. E. Lawrence, J. Appl. Crystallogr., 2011, 44, 213–215.
82. J. Hoja, A. F. Sax and K. Szalewicz, Chem. – Eur. J., 2014, 20, 2292–2300.
83. C. Cadena, Q. Zhao, R. Q. Snurr and E. J. Maginn, J. Phys. Chem. B, 2006, 110, 2821–2832.
84. V. Štejfa, V. Pokorný, A. Mathers, K. Růžička and M. Fulem, J. Chem. Thermodyn., 2021, 163, 106585.
85. J. A. Baird, B. Van Eerdenbrugh and L. S. Taylor, J. Pharm. Sci., 2010, 99, 3787–3806.

---

Footnote

† Electronic supplementary information (ESI) available: Numerical lists of calculated excess properties. Additional structural aspects. Detailed results of NCI analyses. Simulated MSD trends. Numerical lists of calculated T_g. See DOI: https://doi.org/10.1039/d4cp03557g

---