Narges
Karimi
a,
Maryam
Heydari Dokoohaki
a,
Amin Reza
Zolghadr
*a and
Axel
Klein
*ab
aDepartment of Chemistry, School of Science, Shiraz University, Shiraz 71946-84795, Iran. E-mail: nargec_karimi@yahoo.com; maryamhdi5844@gmail.com; arzolghadr@shirazu.ac.ir
bUniversity of Cologne, Faculty for Mathematics and Natural Sciences, Department of Chemistry and Biochemistry, Institute for Inorganic Chemistry and Materials Chemistry, Köln, Germany. E-mail: axel.klein@uni-koeln.de
First published on 30th June 2025
Molecular dynamics (MD) simulation and density functional theory (DFT) analyses were carried out on the solubilization of the heteroaromatic drugs allopurinol, losartan, and omeprazole in reline, which is a deep eutectic solvent (DES) composed of choline chloride/urea in 1:
2 molar ratio. The aim was to gain a deep understanding of the molecular properties of the different components in the solutions and the dynamic interactions between them. Especially, the impact of hydrogen bonding and π-stacking interactions between rings planes in drug⋯drug interactions leading to unwanted aggregations was carefully evaluated using combined angular/radial distribution functions and also natural bond orbital theory. Addition of the drugs to the DES led to decreased conductivity and diffusion coefficients (ranging from 1.52 to 1.07 × 10−12 m2 s−1). The effects of the small allopurinol on the DES structure are more pronounced compared with the larger losartan and omeprazole molecules. This is in line with the formation of aggregates, which are markedly smaller for allopurinol (∼2.6 molecules in average) than for omeprazole (∼4.3) and losartan (∼5.5). DFT-calculated dimers of the drug molecules are stabilized by π-stacking with energies ranging from −10 (allopurinol) to −32 kcal mol−1 (losartan) and show interplanar distances ranging from 0.36 to 0.47 nm. Our observations might pave the way for designing optimized DES for drug delivery.
Aromatic and heteroaromatic moieties are frequently found in drugs and an estimate says that oral drugs have an average of 1.6 aromatic rings in their structures.4 Detailed analyses have emphasized that rising the number of aromatic rings in their molecular structure of drugs lowers their overall developability/applicability, by decreasing their solubility in water and raising their lipophilicity. However, an important factor of aromatic and heteroaromatic structures in drug discovery is the observed enhancement of compound potency by increasing the ligand⋯receptor binding energy through reducing the entropy term.3
The biopharmaceutical classification system (BCS)5,6 sorts drugs along their solubility and their permeability (Fig. 1) and for drug development both aspects must be considered. The use of co-solvents to increase the solubility of drugs with inherently low water solubility has been previously employed including ionic liquids (IL)7 and deep eutectic solvents (DES).8–11 Various use of DESs for therapeutic applications has been reported in the last two decades.10–13 including the solubilization of heteroaromatic drugs such as the nonsteroidal anti-inflammatory drugs celecoxib or piroxicam by the use of choline-based DESs.11–14 Very recently, the mechanism of solubility improvement of some drugs belonging to BCS class II was studied in natural DESs based on solid organic acids, proline, choline chloride, and sugars, which revealed the formation of drug-DES supramolecular assemblies through hydrogen bonding (H bonding).15 Therapeutic DESs (TDES) and volatile DESs (VODES) have recently been explored as novel approaches for enhancing the solubility of active pharmaceutical ingredients (APIs) and controlling APIs crystal growth in a scalable way at room temperature, with the additional possibility of control over polymorphism.2,16
Molecular dynamics (MD) simulations have evolved as a powerful tool to describe the complex behavior of the components in such DES.12,17–26 An early MD study was devoted to the structure and dynamics of a 1:
2 eutectic mixture of choline chloride (ChCl) and ethylene glycol.23 A first-principles MD study on reline and its mixture with water was followed in 201825 and recently refined.19 Very similar to these studies, an MD study on dimethylsulfoxide (DMSO) as an aprotic solvent and reline as DES showed that DMSO has only minor effects on the intermolecular interactions between the reline components in comparison to water.24 Furthermore, while water molecules preferentially interacted with chloride ions, the main interactions of DMSO were shown to be with the urea molecules.24 In a series of ab initio molecular dynamics (AIMD) studies, the solvation shell structures of ammonia, aromatic and aliphatic compounds, and also CO2 and SO2 in DES were reported recently.20–22
The dissolution of heteroaromatic drugs in DES was the topic of a theoretical study on the dissolution of lidocaine (Scheme 1) in two selected DES.26 Very recently, the solubility of celecoxib (Scheme 1) in two ChCl-based DES containing ethylene glycol and 1,2-propanediol hydrogen bond (H bond) donors (HBD) was investigated using MD simulations showing that attractive interactions between the drug and the DES were prevalent over water⋯celecoxib interactions.12
![]() | ||
Scheme 1 Schematic structures of choline chloride (ChCl), lidocaine, celecoxib, losartan, omeprazole, and allopurinol. |
We selected losartan, omeprazole, and allopurinol (Scheme 1), all containing heteroaromatic groups, from the list of the 100 most important drugs.27 Losartan is used as a medicine to treat high blood pressure and prevent heart damage and belongs to class III drugs in the BCS classification (Fig. 1).28 Omeprazole is a stomach acid inhibitor and is found in Class II of the BCS classification.6,15,29 Omeprazole has low stability in acidic solutions, poor water solubility, and undergoes first pass metabolism markedly reducing its concentration before reaching the systemic circulation. All of these factors result in low bioavailability.29,30 Allopurinol is effective in treating gout by reducing the amount of uric acid in the blood and belongs to class IV of BCS classification.31 Due to its low permeability and limited solubility, allopurinol has limited therapeutic effects when used through oral or rectal administration.32 A previous experimental study has shown that the aqueous solubility of allopurinol in various ChCl-based deep eutectic solvents (DESs) combined with urea, glycerol, citric acid, and glycol increases as the weight fraction of DESs increases.33 In contrast to a good number of experimental studies,8–10,13–15,17,33,34 insight on the molecular level mechanism of drug-DES interactions and the role of DESs in the solubility of various drugs through molecular dynamics (MD) studies remained scarce.12,18,19,26
We herein report on an MD study on the solvation properties of losartan, omeprazole, and allopurinol in the reline DES (ChCl/urea 1:
2 mixture), where we observed some degree of self-association through π-stacking and further intermolecular interactions of the heteroaromatic structures.34–37 However, the aggregates are transient and limited in size, rather than persistent or large. Density functional theory (DFT) computational methods were also used in this study to gain insight into the hydrophobic π-stacking interactions. To our knowledge, this is the first study to provide a molecular-level picture of how solvation and aggregation mechanisms vary across different drug scaffolds in a DES medium. By integrating MD simulations with DFT, our work offers new insight into the subtle balance between solvation, self-aggregation, and drug⋯solvent interactions in DES, which is crucial for guiding rational drug formulation and delivery system design.
![]() | ||
Fig. 2 DFT-optimized structures of the DES components (a and b) and the heteroaromatic drugs (c–e) with atom labels. |
The MD topology for each component was generated using the automated topology builder (ATB). MD modeling was performed with the GROMACS 5.1.2 program40 using the GROMOS96 53A6 force field, which was chosen for its reliability in capturing H bonding and solvation behavior in polar and H bond-rich environments,41–45 employed periodic boundary conditions in xyz directions with a short-range cut-off of 1.2 nm. The PME (Particle Mesh Ewald) procedure was utilized to model the electrostatic interactions of atom pairs further apart than 1.5 nm.46
A ChCl/urea 1:
2 molar ratio was modeled using 1000 ion pairs of ChCl and 2000 urea molecules, randomly distributed in a low-density box. This system size was chosen to provide good statistical sampling while maintaining the correct eutectic ratio. The system was energy-minimized and equilibrated to allow the DES structure to form. After equilibration, the simulation box was enlarged to accommodate drug molecules. 12 drug molecules were then randomly introduced into the DES simulation box to construct a representative model system as a compromise between capturing statistically relevant intermolecular interactions and maintaining computational treatability. A second equilibration was performed to allow the drugs to interact with the solvent before starting the production simulations. As previous studies have shown, only a small number of molecules must be treated within the simulation box to achieve relevant data.41–45,47–51
To equilibrate each system five steps were executed: (1) Primary energy minimization was carried out to eliminate close contacts. (2) An equilibration NVT (constant number of particles, volume, temperature) run was executed beyond 5 ns in 1 × 10−3 ps time steps, while the system was linked to a velocity-rescaling thermostat using 0.1 ps coupling constants.52 (3) Consequent configurations were utilized to launch the NPT (constant number of particles, pressure, temperature) runs for 5 ns using the Parrinello–Rahman barostat.53 (4) After the primary equipoise, the systems were captured as the initial configurations for the production simulations of 10 ns, with 2 × 10−3 ps time steps at 450 K, ensuring that the mobility of molecules is high enough for achieving a reliable equilibration. (5) Simulations were carried out for 20 ns at 298 K. Considering the simulation annealing steps in which the temperature of each system was decreased at intervals of 50 K from 400 K to 298 K, the simulation was performed for 500 ps under constant NPT conditions at each temperature. The overlap of RDF graphs considering various atom⋯atom interactions at different simulation times proves that the systems are well equilibrated (see Fig. S1–S3, ESI†).42,48–51
To address potential biases from temperature scaling in the NVT/NPT ensembles, we performed additional simulations in the NVE (constant N, V, and E) ensemble for 20 ns. The NVE ensemble allows for energy conservation without external thermostats, providing a more direct evaluation of the system's dynamical properties.
Radial distribution function (RDF) analysis was carried out by means of the gmx rdf tool in Gromacs40 to provide insight into the probability distribution of each DES or drug molecule around the other molecule.
Combined radial/angular distribution function (CDF) plots are a potent analysis for the study of H bonding and CDF analyses were performed by combining radial and angular distribution functions by trajectory analyzer and visualizer (TRAVIS) package.54 With this method, both inter- and intramolecular interactions have been previously studied.55
The formation of H bonds between the DES components as well as in DES⋯drug and drug⋯drug interactions was calculated using the gmx hbond criterium through the GROMACS package.40 The lifetime of H bonds was evaluated through time autocorrelation functions (TCF), using the TRAVIS software.54 The size of aggregations of the drug molecules was estimated using gmx clustsize through the GROMACS package.40 A histogram of aggregate (cluster) sizes was obtained from the simulation trajectory files with a threshold of lower than 0.35 nm indicative for aggregates.
To examine the dynamics of the drugs in the DES, the center of mass (COM) mean-square displacement (MSD) dataset was obtained during the MD trajectories using gmx msd through the GROMACS program.40 Also, the diffusion coefficients Di were computed. In general, the average distance squared when an atom or molecule moves away from beginning point throughout the time interval τ (eqn (1)).
![]() | (1) |
![]() | (2) |
![]() | (3) |
The molar electrical conductivity of the solutions was determined using the Nernst–Einstein (eqn (4)).
![]() | (4) |
Viscosity calculations were performed for each system using the Einstein–Stokes equation (eqn (5)).
![]() | (5) |
Interaction energy = EDimer − 2EMonomer | (6) |
![]() | (7) |
In general, the simulated snapshots as pictures of individual configurations or even averaged structures do not comprise any numerical details of full simulation length. To explore the structure of the solutions in detail and to examine how the drugs are dissolved in DES on a molecular level, we used radial distribution function (RDF), spatial distribution function (SDF), and combined radial/angular distribution function (CDF) analyses.
These interactions align closely with findings from prior investigations into pure reline. In these studies, the essential role of reline in achieving the eutectic properties is attributed to the specific H bonding interactions between chloride anions and the hydroxyl H atom of the choline cations, as well as with the H atoms of urea.12,18,19,62,63
The H3 atoms of allopurinol show a very high RDF with the Cl− anions slightly below 0.2 nm, while further, less intense interactions of H3 and H2 with O(Urea) and O(Chol) atoms cluster range between 0.18 (H3⋯O(Urea)) and 0.25 nm (H2⋯O(Urea) (Fig. 5c). While these drug⋯DES interactions are short-ranged, the drug⋯drug interactions O1, N1, and N4 around H2 and H3 are slightly longer but probably dominant as can be seen from their far higher RDFs (Fig. 5d and Fig. S2, ESI†). The strongest peaks are found for H2⋯N3 and H3⋯N2 at distances above 0.35 nm. A third, less intense RDF is observed for H2⋯O1 at 0.40 nm. Some intense RDFs are found for C1⋯C1, C1⋯C5 and C2⋯C3 at 0.47 nm (Fig. S5, ESI†), which lies in the typical range for intermolecular π-stacking interactions of either π–π or CH–π type.62,63 These RDFs exhibit lower intensity when compared with the prevalent drug⋯drug H bond interactions.
Additionally, there are peaks of low probability for H(Chol)⋯O1 and H(Chol)⋯O3 at 0.18 nm. The interaction of omeprazole with urea is less important for the solvation, but the H(Urea)⋯O2, ⋯O3, and ⋯N3 interactions are all found around 0.18 nm and represent an extended H bond network (Fig. 7b). The acidic N–H1 maintains similar weak, but short interactions with Cl, N(Urea), and O(Urea) (Fig. S7, ESI†). As for allopurinol and losartan, the drug⋯drug H bonding network in the DES/omeprazol system seems to be dominant over the DES⋯drug H bonds as can be interpreted from the short H1⋯O2 and H1⋯N2 interactions that comes along with a longer but more intense H1⋯N3 interaction (Fig. 7c). π-stacking interactions can be inferred from intense the C3⋯C6 RDF peak and further smaller peaks in the range 0.35 to 0.5 nm (Fig. 7d). The highest RDF approaches a level of intensity similar to that of the most potent drug⋯drug H bond.
The sharp rise in the first peak intensity in the H(chol)⋯Cl RDF plot in the presence of drug molecules (Fig. 4, 5a, 6b, and 7a) can be attributed to the reorganization of the solvent structure caused by drug⋯solvent interactions. The drugs seem to induce stronger H bonding or van der Waals interactions with the choline cation and chloride anions, leading to a denser local environment. This results in a higher first peak in the RDF, suggesting more frequent or stronger interactions between H(chol) and Cl. We found that the first peak height in the H(chol)⋯Cl RDF was lowest for losartan, slightly larger for allopurinol, and markedly larger for omeprazol. Thus, the size of the drug molecules is not decisive.
The intense and strong peaks in the RDFs for the H bonding between the H atoms of drug molecules and chloride align with previous studies, which revealed that H bonding between chloride and various H atoms in β-cyclodextrin, MeOH, EtOH and water occurs at 0.2 nm.62–66
More details of the H bonding were revealed by CDF analyses and in the following, an H bond is defined when the D–H⋯A angle is between 160 and 180° and the A⋯H distance is lower than 0.3 nm. Pronounced H bonds are found in the DES/allopurinol system between allopurinol H2, H3, and O1 and the Cl, O(Urea), H(Urea), and H(Chol) DES components (Fig. 8 and Fig. S8, ESI†). The H3⋯Cl distance is found around 0.19 nm with a linear N3–H3⋯Cl (180°) H bond (Fig. 8a). The strongest drug⋯urea interaction is found for H3⋯O(Urea) with an N3–H3⋯O(Urea) angle of 178° (Fig. 8b). A little bit less pronounced drug⋯urea and drug⋯choline H bond interactions are found for H(Urea)⋯O1 (Fig. 8c) and H(Chol)⋯O1 (Fig. 8d) with distances less than 0.2 nm and angles of 175°. This is in complete agreement with the conclusions drawn from the RDF plots (see above).
![]() | ||
Fig. 8 CDF plots for DES/allopurinol system (a–d). Red relates to low probability and dark blue relates to high normalized probability. |
The same is true for the DES/losartan system (Fig. S9, ESI†). Here, the most pronounced interactions occur between the DES H atoms H(Chol)/H(Urea) and the N1/O1/N4 atoms of losartan. By far the highest probabilities were found for the O1 atom with H(Chol) and N(Urea) with the O(Chol)–H(Chol)⋯O1 angle at 178° and the N(Urea)–H(Urea)⋯O1 angle at 180°. In turn, for the acidic H1 and H6 atoms of losartan, the strongest H bonds are found with the chloride ion of the DES and O of the choline cation (H1⋯Cl and H6⋯Cl at 0.18 nm with angles ranging from 165 to 180°, and 0.17 nm with an 180° angle for H6⋯O(Chol)). Overall, these H bonds are similar in character to those of allopurinol. However, there is a tendency of higher probability for the corresponding N–H⋯Cl interaction and an additional very strong N–H⋯O(Chol) H bond for the losartan systems, which is unmatched in the DES/allopurinol system.
CDF plots of the DES/omeprazole system show strong O(Chol)–H(Chol)⋯O2 interactions with angles close to 180° and sub-0.2 nm distances, confirming the important role of the SO group for H bonding (Fig. S10, ESI†). Less intense interactions were detected for N1–H1⋯O(Urea) and O(Chol)–H(Chol)⋯O3, both showing 178° angles and distances below 0.2 nm.
Previous studies on solvation of rather hydrophilic (ammonia)20 and hydrophobic (methane)67 solutes showed, that both are mainly surrounded by the Chol hydroxyl groups, with NH3 H bonds having specific geometric criteria (≤0.25 nm and H(Chol)-O(Chol)-N(NH3) angle of ≤30°). Most OH groups in Chol show a tangential orientation at 60° to 80° and are about 0.4 nm from the hydrophobic solute. In arginine-based DESs, anesthetic drug solubility studies revealed strong HBA (hydrogen bond acceptor) and HBD (hydrogen bond donor) interactions at 0.19 nm and 160°.68
Insight into the three-dimensional space came from SDF analysis (Fig. 9 and Fig. S11, S12, ESI†). In the DES/allopurinol system, the spatial density of the choline cation is predominantly located near the H3 and N3 atoms of allopurinol. The chloride anions are exclusively grouped close to the allopurinol H3 atom. Urea approaches the drug molecule clearly at its O1 atom. Allopurinol molecules are grouped in a clamp-like fashion around the drug molecule with prevailing presence in well-defined regions above and below the double-ring plane. So, for allopurinol, all intermolecular interactions (H bond and π-stacking) are rather localized at specific parts of the molecule.
The SDF analysis of DES/losartan (Fig. S11, ESI†) is less simple. Choline, Cl−, and urea show several components of spatial density around the losartan molecule. For the chloride ion the major component is located very close to the OH group of losartan, in line with the RDF (Fig. S6, ESI†) and CDF (Fig. S9a, ESI†) plots. Urea SDF densities are vastly delocalized, which is in line with the plethora of H bonds of the carbonyl H-accepting and the NH2 H-donating units with various OH, NH and even CH groups of the drug, that can be assumed from the RDF and CDF plots. The SDF density of the drug found in close proximity to the alkyl chain and in the space far from the losartan molecule, can be interpreted as dimerization, leading to special closeness of the alkyl chains.
The SDF pictures of the DES/omeprazole system (Fig. S12, ESI†) show several density components for the choline cation, the chloride, and urea. Interestingly, the urea density clamps the drug molecule between the SO and the MeOPh groups, while the density of omeprazole is essentially located in the binding pocket of the concavely arranged omeprazole molecule, which also points to marked dimerization.
These observations match an earlier study on NH3, showing urea and chloride near the NH3 H atoms and Chol near the N atom.20 Chloride ions in alcohol–reline mixtures primarily interact with the alcohol hydroxyl group H atoms.65 This is fully in line with the isosurface distribution of chloride in the current study. Arginine-based DESs around anesthetic drugs show localized HBA and HBD molecules with minimal neighboring drug distribution.68 In ChCl:lactic acid DES, Chol cations surround the lidocaine phenyl, NH, and CO groups, (Scheme 1) with less prominent chloride ions, and no significant drug self-clustering.26 Lidocaine, due to its amphiphilic nature, could potentially replace one of the DES components and form gel-like structures under certain conditions, as suggested by a previous study.69 While this was not the focus of the present study, the possibility of gel formation in DESs incorporating lidocaine represents an interesting future research direction. This could have implications for improving drug release properties and stability in DES-based drug delivery systems.
The formation of ChCl⋯urea, ChCl⋯drug, and urea⋯drug H bonds was calculated using the gmx hbond criterium through the GROMACS package (Fig. S13 and Table S1, ESI†).40 In the presence of drug molecules, the number of H bonds between solvent components is only slightly reduced from 812 (per time frame) to 808 for omeprazole, to 795 for allopurinol and to 773 for losartan, while essentially urea⋯drug (27 to 38) and to a smaller extent ChCl⋯drug (2 to 16) H bonds are formed (Table S1, ESI†). Lifetimes of these H bonds as determined through time autocorrelation functions (Fig. S14, ESI†) range from 3.8 to 5.5 ms and slightly decrease along the series: choline > urea > chloride (Table S2, ESI†). Omeprazole forms significantly more H bonds with both urea and ChCl, compared to allopurinol and losartan. Especially the number of H bonds with ChCl is far higher (16 compared with 2). This is probably due to the unique and very polar sulfinyl SO group, adding to the acidic NH function. The latter is also found in allopurinol and losartan.
To further study possible π-stacking interactions, the CDFs between different aromatic and heteroaromatic rings are explored by monitoring both the distance (d) and the angle (α) (Fig. 10 and Fig. S15, S16, ESI†). As in previous studies,70–73 we defined π–π stacking along two criteria: (i) d < 0.5 nm and (ii) 0° < α < 40° or 150° < α < 180°.
![]() | ||
Fig. 10 CDFs (angular radial distribution functions) of omeprazole ring planes in six combinations. Here, colors from bright red to dark blue indicate the probability distribution of stacking angles. CoR and RN are the center of ring and ring normal, respectively, for the labeling of the rings, see Fig. 2. |
For allopurinol, the highest CDFs lie in the range of 0.75 to 1 nm distance with 40° and ∼140° angles (Fig. S15, ESI†). This is in line with the intense atom⋯atom RDF, observed for C2⋯C5 at 0.75 nm (Fig. S5, ESI†). Only a very weak RN1 (RNi = ring normal of ring i)/CoR1⋯CoR1 (CoRi = center of ring i) interaction is observed at around <0.5 nm with angles around 50 and 130°. For losartan, very similarly to allopurinol, most of the interactions lie markedly above 0.5 nm, only the RN1 (ring normal, ring 1)/CoR1⋯CoR1 (center of ring 1), CoR1⋯CoR3 (center of ring 3), and CoR1⋯CoR4 (center of ring 4) interactions lies at around 0.47 nm with 0° < α < 40° or 150° < α < 180° (Fig. S16, ESI†). The interaction of the RN2/CoR2⋯CoR4 occurs at 0.47 nm with 0° < α < 30° or 160° < α < 180° in line with the intense RDF C11⋯C20 of the tetrazole ring (R4) and biphenyl (R2) at 0.47 nm (Fig. 6d). For the biphenyl moiety, a strong interaction of the RN2/CoR2⋯CoR2 is observed. For omeprazole, well defined maxima slightly above 0.4 nm with angles of around 50 and 140° for the RN2/CoR2⋯CoR2 interaction fulfill roughly our π–π stacking criteria (Fig. 10). This is fully in line with the intense C6⋯C6 RDF peak (Fig. 7d). The same is true for the intense RN1/CoR1⋯CoR1 interaction observed at distances <45 nm and angles varying from 10 to 170°, which fully match the shorter RDF peak of C9⋯C11 at around 0.42 nm. Thus, it seems that omeprazole and losartan exhibit pronounced π–π stacking interactions in the DES solution, while allopurinol shows only weak π–π stacking.
The average size of aggregations of the drug molecules was estimated using gmx clustsize through the GROMACS program (Fig. S17, ESI†).40 Production MD simulations were extended to 100 ns to ensure sufficient sampling of interactions, particularly those involving transient aggregation events. The average drug cluster sizes were 2.585 for allopurinol (136.112 g mol−1), 5.511 for losartan (422.91 g mol−1), and 4.255 for omeprazole (345.42 g mol−1). It should be noted that these values reflect aggregation behavior under the specific simulation conditions used in this study (12 drug molecules per box), and would be expected to vary with changes in drug concentration.
Overall, the aggregation numbers remain small throughout the simulations, suggesting that although transient self-association occurs, the extent and stability of the aggregation is limited. This highlights the efficiency of reline DES in solvating these drugs and preventing the formation of large or insoluble aggregates. The smaller aggregate size of allopurinol suggests that it remains more dispersed in solution, which can facilitate mass transfer, whereas larger aggregates of losartan and omeprazole may reduce overall mobility within the DES.
These findings are consistent with previous studies showing that the solvation of doxorubicin varies across solvents, affecting its aggregation size,42 and that tetracycline forms larger clusters in salt solution compared to EtOH due to hydrophobic interactions.74 Additionally, 1-ethyl-3-methyl imidazolium acetate IL effectively disrupts Dopamine aggregation, resulting in smaller aggregates and improved dissolution.75
Experimental studies have consistently demonstrated that DESs can support the formation of nanostructures through self-assembly of both surfactants and drug molecules confirmed by Tyndall scattering and dynamic light scattering (DLS). For instance, the self-aggregation behavior of various surfactants in reline-based DES has revealed a broad range of aggregate sizes (1–200 nm), with some micellar systems comprising as many as ∼60–75 surfactant molecules.76–80 In the context of pharmaceuticals, experimental studies have shown that depending on the nature of the drug and the DES system, drug molecules may either remain molecularly dissolved or form colloidal aggregates.80 Further, it was reported that several poorly water-soluble drugs form colloidal dispersions (10–200 nm) in thymol-based DES media.81 In contrast, other drugs did not exhibit aggregation and remained truly dissolved.81 Our MD simulations, though limited to small cluster sizes considering a model system, offer insight into the molecular-level interactions that may initiate such behavior, and capture the initial stages of such aggregation phenomena.
System | D chol+ | D Cl− | D urea | D drug | t chol+ | t Cl− | Λ |
---|---|---|---|---|---|---|---|
a Diffusion coefficients obtained from the MSD-time curves. Molar electrical conductivity Λ determined using the Nernst–Einstein equation (eqn (4)). t = transference (or transport) numbers outlined according to the anionic and cationic diffusion coefficient values ![]() |
|||||||
DES | 1.52 ± 0.15 | 0.76 ± 0.02 | 0.93 ± 0.11 | — | 0.67 | 0.33 | 8.57 ± 0.23 |
DES/allopurinol | 1.14 ± 0.05 | 0.61 ± 0.06 | 0.80 ± 0.09 | 0.48 ± 0.06 | 0.65 | 0.35 | 6.58 ± 0.14 |
DES/losartan | 1.12 ± 0.13 | 0.58 ± 0.04 | 0.76 ± 0.12 | 0.28 ± 0.03 | 0.66 | 0.34 | 6.39 ± 0.09 |
DES/omeprazole | 1.07 ± 0.08 | 0.55 ± 0.07 | 0.74 ± 0.04 | 0.21 ± 0.01 | 0.66 | 0.34 | 6.09 ± 0.12 |
Previously reported experimental values for the diffusion coefficients of the studied drugs across different media are as follows: allopurinol in neutral aqueous medium (930.7 × 10−12 m2 s−1),82 omeprazole in MeOH (188.0 × 10−12 m2 s−1)83 and an aqueous buffer supporting electrolytes (231.0 × 10−12 m2 s−1),84 and losartan in aqueous solution (260.0 × 10−12 m2 s−1), DMSO (160.0 × 10−12 m2 s−1) and in a micellar environment (43.0 × 10−12 m2 s−1).85 The lower diffusion coefficients observed in our simulations are expected, given the significantly higher viscosity and complex H bonding network of the reline DES compared to typical organic solvents. Additionally, strong drug⋯solvent and drug⋯drug interactions in DES, such as H bonding and π-stacking, can further hinder molecular mobility. Nevertheless, the observed trends remain qualitatively consistent and useful for comparing relative mobility among the drugs in the DES environment.
To assess the reliability of our results, we compared the dynamical properties obtained from NPT (Table 1) and NVE (Table S3, ESI†) simulations. The NVE simulations further confirmed the accuracy of the observed dynamic properties, demonstrating that temperature fluctuations do not significantly alter the diffusion coefficients or other transport parameters. This reinforces the reliability of our results and the validity of the NPT ensemble choice in capturing the key behaviors of the studied systems. The differences in the diffusion coefficients and other dynamic properties were minimal, suggesting that our original findings are robust and unaffected by the choice of ensemble.
As shown in Table 1, the values of cationic transport numbers (ti) for the choline cations are larger than 0.5 both in the presence or absence of the drug molecules in the DES. This means that the cation plays a major role in carrying electric current in these systems.
The calculated molar electrical conductivity of the solutions showed that addition of drug molecules to the DES decreases the mobility of ionic species, which negatively impacts the ionic conductivity of the mixtures. The dissolution of drug molecules lowered the electrical conductivity of the solutions along the series: pure DES > DES/allopurinol > DES/losartan > DES/omeprazole. Allopurinol has a relatively minor effect on the ionic mobility of the DES components due to its smaller size. Losartan and omeprazole show larger aggregation, which possibly leads to a stronger disruption of the ion movement in the solution, and reducing the free ions available for conductivity.86 The decrease in conductivity upon drug addition, further induces alterations in the DES structure, suggesting potential effects on the mass transfer of the solutes.9
Reline is known for its high viscosity, which can significantly impact the drug solubility and diffusion dynamics. Hence, viscosity calculations were performed for each system using the Einstein–Stokes equation (eqn (5)). The calculated viscosity for the pure DES of 1037 mPa s at 298.15 K, lies within the range of experimental values (750–1750 mPa s), which validates the computational model used in this study.87–90
Upon the addition of the drugs, the calculated viscosity increased to 1277 mPa s for allopurinol, 1336 mPa s for losartan, and 1395 mPa s for omeprazole. The larger the aggregates, the more resistance there is to the molecular mobility, leading to higher viscosity values. The dissolution of drug molecules also likely affects the solvent structure of the DES, leading to a more ordered or structured arrangement of the DES components around the drug molecules. This structuring effect can also increase viscosity, as the solvent components become less mobile. The increased viscosity upon the addition of the drugs, particularly the larger drugs, implies a more “viscous” environment that could slow down the mass transfer process.91–93
For the bicyclic allopurinol structure, which comprises a pyrazole ring fused to a hydroxy-substituted pyrimidine ring, we observed an intermolecular rotation (∼120°) in the stacked dimer with an energy of about −10 kcal mol−1 compared with the monomer. The interplanar distance was calculated to be around 0.36 nm. The shortest CoR1⋯CoR2 interaction is found at around 0.45 nm in the DES, based on the RDF and CDF plots (Fig. S5 and S15, ESI†). In the dimer model for losartan, the molecules were stacked through the biphenyl rings resulting in a markedly higher interaction energy of about −32 kcal mol−1 compared with that of allopurinol. The interplanar distances were calculated to lie between 0.36 and 0.47 nm, in good agreement with RDF and CDF interacting distances of 0.47 nm in the DES (Fig. 6 and Fig. S13, ESI†). For the omeprazole dimer, both the substituted pyridine and the substituted benzimidazole rings are stacking with each other, the energy gain being about −26 kcal mol−1. The interplanar distances ranging between 0.34 and 0.42 nm, closely match with the primary π–π stacking interactions depicted in Fig. 7d and 10. It is important to note, that the dimer interaction energies were computed in the gas phase, and thus do not account for solvent effects present in the reline DES environment. As a result, the reported values likely represent upper-bound estimates of interaction strength. These gas-phase calculations provide qualitative insight into the relative tendencies of the drugs to self-associate, which were further supported by MD simulations performed in the presence of solvent.
The formation of dimers is thus favorable in all three systems with an increasing energy gain along the series: allopurinol < omeprazole < losartan. Among the three drugs studied, losartan also showed the highest average aggregation number in MD simulations, which can be attributed to its larger molecular size, extended aromatic system, and hydrophobic character. The biphenyl system in losartan provides a large, conjugated, and rigid π-surface that favors π–π stacking, a key driving force for self-aggregation in non-aqueous environments like DES. This structural feature, combined with the presence of additional polar and H bonding groups, promotes stable dimer and small aggregate formation. Notably, an experimental study94 has also reported the formation of dimeric degradation products of losartan under acidic or neutral conditions, in which nucleophilic substitution reactions occur between functional groups of losartan molecules. While these degradation pathways differ from our observed non-covalent aggregates, they suggest that losartan has a predisposition for molecular self-association, especially under certain solvent conditions.
A previous DFT-AIM study found stacking interaction energies between drug fragments such as indole, benzofuran, and benzothiophene and the natural nucleobases guanine, adenine, thymine, and cytosine falling below 12 kcal mol−1, particularly those involving one or two rings.95 N-based organic fused heterocycle stacked dimeric systems like allopurinol have been reported to exhibit effective π–π stacking interactions in the gas phase, with interaction energies typically below 10 kcal mol−1.59 Calculations on imidazole stacking interactions, have shown inter-ring distances ranging from 0.38 to 0.43 nm,96 which aligns with the π–π stacking distances between the benzimidazole rings in omeprazole calculated by us.
To assess the physico-chemical properties of the dimer models, we compared the HOMO–LUMO energies (Table S4, ESI†). The HOMO–LUMO gap, representing the overall reactivity, decreased in the order: allopurinol < omeprazole < losartan. This is in line with the highest energy-gain through dimerization of losartan within this series. The value of global hardness for the losartan dimer (1.487 eV) compared with allopurinol (1.893 eV) and omeprazole (1.598 eV) agrees also with the idea that the negative dimerization energies go along with reduced reactivity. The losartan dimer structure showed the highest softness within the series and less resistance against changing in its electronic configuration.
The most important donor–acceptor interactions within the drug dimers obtained from NBO analysis are reported in Tables 2 to 4 as well as in Tables S5, and S6, ESI.† π–π stacking interactions are found for all the dimers, whereas CH–π stacking interactions are observed for losartan and omeprazole (Tables S5 and S6, ESI†). As expected, the overall strength of the π–interactions were considerably larger for losartan and omeprazole compared with allopurinol, in good agreement with the RDF and CDF results (Fig. 5, 6, 7, 8, 10 and Fig. S5 to S10, ESI†) from the MD study.
Dimer structure | ||
---|---|---|
π–π stacking interactions | ||
Donor | Acceptor | E (2) (kcal mol−1) |
a BD(2) and BD*(2) are bonding and anti-bonding natural orbital of π and π* type, respectively. a and b refer to the two molecules as shown. | ||
aBD(2) N1⋯C1 | bBD*(2) C5⋯N4 | 0.36 |
bBD(2) C5⋯N4 | aBD*(2) C2⋯O1 | 0.11 |
bBD(2) C3⋯C4 | aBD*(2) C3⋯C4 | 0.14 |
bBD(2) N1⋯C1 | aBD*(2) C5⋯N4 | 0.36 |
aBD(2) C5⋯N4 | bBD*(2) C2⋯O1 | 0.11 |
Dimer structure | ||
---|---|---|
π–π stacking interactions | ||
Donor | Acceptor | E (2) (kcal mol−1) |
a BD(2) and BD*(2) are bonding and anti-bonding natural orbital of π and π* type, respectively. a and b refer to the two molecules as shown. | ||
aBD(2) C12⋯C13 | bBD*(2) C11⋯C12 | 2.29 |
aBD(2) C12⋯C13 | bBD*(2) C9⋯C10 | 3.94 |
aBD(2) C12⋯C13 | bBD*(2) C17⋯C18 | 1.00 |
aBD(2) C12⋯C13 | bBD*(2) N4⋯N5 | 2.66 |
aBD(2) C12⋯C13 | bBD*(2) N3⋯C20 | 1.50 |
aBD(2) C8⋯C9 | bBD*(2) C11⋯C12 | 0.48 |
aBD(2) C10⋯C11 | bBD*(2) C9⋯C10 | 0.13 |
bBD(2) C14⋯C19 | aBD*(2) C12⋯C13 | 0.80 |
bBD(2) C14⋯C19 | aBD*(2) C17⋯C18 | 0.89 |
aBD(2) C15⋯C16 | bBD*(2) C17⋯C18 | 0.41 |
bBD(2) C14⋯C19 | aBD*(2) C15⋯C16 | 1.06 |
Dimer structure | ||
---|---|---|
π–π stacking interactions | ||
Donor | Acceptor | E (2) (kcal mol−1) |
a BD(2) and BD*(2) are bonding and anti-bonding natural orbital of π and π* type, respectively. a and b refer to the two molecules as shown. | ||
aBD(2) C11⋯C12 | bBD*(2) C9⋯C10 | 0.19 |
aBD(2) C6⋯N3 | bBD*(2) C5⋯C7 | 13.77 |
aBD(2) C6⋯N3 | bBD*(2) C3⋯C4 | 43.62 |
aBD(2) C6⋯N3 | bBD*(2) C2⋯C8 | 9.05 |
aBD(2) C5⋯C7 | bBD*(2) C3⋯C4 | 4.42 |
aBD(2) C5⋯C7 | bBD*(2) C2⋯C8 | 0.71 |
aBD(2) C3⋯C4 | bBD*(2) C5⋯C7 | 2.37 |
aBD(2) C3⋯C4 | bBD*(2) C2⋯C8 | 1.99 |
aBD(2) C3⋯C4 | bBD*(2) C3⋯C4 | 28.61 |
aBD(2) C2⋯C8 | bBD*(2) C5⋯C7 | 0.24 |
aBD(2) C2⋯C8 | bBD*(2) C3⋯C4 | 1.19 |
aBD(2) C2⋯C8 | bBD*(2) C2⋯C8 | 0.28 |
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5cp01312g |
This journal is © the Owner Societies 2025 |