Skip to main content


Front. Pharmacol., 20 April 2021
Sec. Translational Pharmacology
Volume 12 - 2021 |

Molecular Dynamics Studies of Therapeutic Liquid Mixtures and Their Binding to Mycobacteria

  • 1LAQV, REQUIMTE, Chemistry Department of NOVA School of Science and Technology, Caparica, Portugal
  • 2Red Glead Discovery AB, Lund, Sweden

Tuberculosis is an highly contagious disease still considered by the WHO as one of most infectious diseases worldwide. The therapeutic approach, used to prevent and treat tuberculosis targets the Mycobacterium tuberculosis complex, comprises a combination of drugs administrated for long periods of time, which, in many cases, could cause several adverse effects and, consequently, low compliance of the patient to the treatment and drug-resistance. Therefore, therapeutic liquid mixtures formulated with anti-tuberculosis drugs and/or adjuvants in tuberculosis therapy are an interesting approach to prevent toxic effects and resistance to anti-tuberculosis drugs. The herein formulated therapeutic liquid mixtures, including ethambutol, arginine, citric acid and water under different molar ratios, were studied through a molecular dynamics approach to understand how ethambutol and arginine could be stabilized by the presence of citric acid and/or water in the mixture. To gain insights on how the uptake of these mixtures into the mycobacteria cell may occur and how a mycobacterial ABC transporter could contribute to this transport, multiple simultaneous ligand docking was performed. Interactions between citric acid and ethambutol involving the carboxyl and hydroxyl groups of citric acid with the amines of ethambutol were identified as the most critical ones. Water molecules present in the mixture provides the necessary network of hydrogen bonds that stabilize the mixture. Molecular docking additionally provided an interesting hypothesis on how the different mixture components may favor binding of ethambutol to an ABC importer. The data presented in this work helps to better understand these mixtures as well as to provide cues on the mechanisms that allow them to cross the mycobacterial cell membrane.


Tuberculosis persists as one of the major causes of morbidity and mortality caused by an infection, worldwide. It is an infectious disease triggered by Mycobacterium tuberculosis complex that is pathogenic in humans and animals and belongs to a group of bacteria named Mycobacteria. This group of bacteria are frequently resistant to several antibiotics, mostly due to slow absorption of the drugs through the highly hydrophobic mycobacterial cell envelope (Viveiros et al., 2003). The mycobacterial cell wall is composed by peptidoglycans, arabinogalactans, and mycolic acids (Torfs et al., 2019), which leads to the formation of an intrinsic barrier that hampers drugs penetration and can contribute to drug resistance to anti-tuberculosis drugs (Silva et al., 2013). In spite of a descending tendency observed on the evolution of tuberculosis over the last years, in 2019 the WHO reported 1.4 million deaths by this infectious disease (Fitzgerald et al., 2015; WHO, 2020). The treatment for tuberculosis includes a combination therapy of multiple antibiotics administrated for several months. One of the major problems in tuberculosis treatment is the emergence of multidrug resistance (MDR) that can occur during the therapy due to chromosomal mutations mechanisms (e.g., single nucleotide polymorphisms), drug-drug interactions, malabsorption of drugs and prolonged administration of multiple drugs, which in turn leads to a low compliance of the patient to the treatment (Forbes et al., 2018; Hameed et al., 2018).

Since this disease continues to present some gaps in the effective treatment, particularly when MDR occurs, new strategies are most welcome. In addition, optimization of the existing ones to reduce toxic and MDR associated events related to tuberculosis therapy may also play an important role in the treatment of the disease (WHO, 2015; Rendon et al., 2017; Sotgiu et al., 2017). One possible strategy to overcome these gaps is the development of shorter regimens of drug administration which could be achieved using alternative and sustainable solvents as therapeutic liquid mixtures (Aroso et al., 2016; Duarte et al., 2017). These mixtures are considered low transition temperature mixtures (LTTMs) and are formed mixing two or more components that at certain molar ratios, often form a clear viscous liquid at room temperature with a lower melting point than those of the initial compounds (Dai et al., 2013; Francisco et al., 2013; Aroso et al., 2016; Abbott et al., 2017). Different theories have been reported to explain which interactions could be involved to promote their stability as liquid forms and their properties. Accordingly, Abbott et al. had proposed hydrogen bonds as important contributors for the decrease of entropic differences of the phase transitions of the components and, concomitantly, for the reduction of their melting point (Abbott et al., 2017). The different atoms present in the molecules and the number of HBDs and HBAs are, hence, among the main factors that could determine the capacity of the different components of the mixture to interact, together with van der Waals and/or electrostatic forces (Abbott et al., 2017; Santana et al., 2019). In this sense, the charge delocalization could modulate the physicochemical properties of these eutectic liquid mixtures, when compared with initial components (Francisco et al., 2013). Furthermore, being water part of the system and acting as a component of the mixture, its presence usually facilitates the preparation and manipulation of the mixture by decreasing their viscosity and, quite often, increasing the stability of the mixture. However, adding more than 50% of water to the system could lead to a disruption of the hydrogen bonds and dilution of the components, creating a solution (Dai et al., 2015; Craveiro et al., 2016). In previous studies, the possibility of pH itself affect the solubility of other components in the mixture was explored and it was observed that the slight variation of pH did not affect the solubility of these eutectic mixtures and due to that, the pH was not considered in molecular dynamics (MD) studies (Santos et al., 2018; Roda et al., 2021).

In the pharmaceutical field, these mixtures have been used over the decades and several studies reported that the combination of active pharmaceutical ingredients (API) with compounds that could interact by van der Walls and electrostatic forces and/or hydrogen bonds, in determined molar ratios, could provide an enhancement of the some properties of the API such as solubility, permeability, and absorption (Evers et al., 1985; Stott et al., 1998; Aroso et al., 2015). One of the examples was reported by Stott and coworkers, they prepared eutectic mixtures combining ibuprofen with different terpenes and formed a therapeutic liquid mixture for enhancement of skin permeation (Stott et al., 1998). The therapeutic liquid mixtures comprise thousands of different combinations, and although their preparation is still based on an approach of trial and error, gaining insights about the interaction between the components of the mixture provides more information towards the successful development of these mixtures. The knowledge about how therapeutic liquid mixtures work and how they could be applied in pharmaceutical and biomedical field remains limited. However, it is a field that has been continuously explored, in an attempt to improve the characteristics of active pharmaceutical ingredients and enhance their pharmacokinetics and pharmacodynamics. To overcome these questions, MD represents an important tool that could boost the understanding of these mixtures, particularly how the components interact with each other to form a stable mixture (Kumari et al., 2018; Gutiérrez et al., 2019), by providing insights on the behavior of atoms and molecules and allowing the prediction of the interactions of a particular system during a pre-determined period of time. In the case of eutectic systems, this type of analysis could give insights about the chemical and physical properties of the systems, such as dynamics of the molecules, radial distribution functions (which describes how density varies as a function of distance from a reference particle), spatial distribution functions of the atoms, self-diffusion coefficients, intermolecular interaction, polar solvation energies, percentage of volume expansion and type of hydrogen-bonding interactions. To this matter, MD simulations had already proven its value by providing valuable knowledge regarding the solubility of active pharmaceutical ingredients (Gutiérrez et al., 2018, 2019) or through the study of the internal structure of deep eutectic solvents at its molecular level (Fetisov et al., 2018; Kumari et al., 2018).

In the current study, a MD approach was used to understand how an anti-tuberculosis drug, ethambutol (EMB), and a natural compound as arginine (ARG), described to act as an adjuvant in tuberculosis therapy (Schön et al., 2011; Farazi et al., 2015; Santos et al., 2018), could be stabilized by the presence of citric acid (CIT) and/or water in the mixture. Furthermore, a molecular docking study was also performed to unveil a possible CIT-promoted entry pathway for EMB through a recently published mycobacterial ABC transporter (Rempel et al., 2020).

Materials and Methods

All molecules were parameterized in OPLS/AA force field (Jorgensen and Tirado-Rives, 1988) using LigParGen server (Jorgensen and Tirado-Rives, 2005; Dodda et al., 2017a; Dodda et al., 2017b) and are depicted in Figure 1. Four systems were considered for evaluation: citric acid monohydrate:arginine:water 1:1:7 (S1), citric acid monohydrate:ethambutol:water, 2:1:10 (S2) and 1:1:5 (S3), and citric acid monohydrate:ethambutol:arginine:water 2:1:1:7 (S4) (Table 1).


FIGURE 1. Therapeutic liquid mixtures components studied in this work. Atom names used through the manuscript are also indicated for clarity.


TABLE 1. Therapeutic liquid mixtures used for the molecular dynamics studies.

To avoid any bias due to system assembly, two smaller systems (S1 and S2) were initially built by creating a solvated box, replacing the number of water molecules by the species of interest and further removing the water molecules in excess to obtain the desired ratio. Both systems were energy minimized, with the system’s temperature being equilibrated using an NVT ensemble during 1-ns followed by a 1000-ns (1-µs) long MD simulation in an NpT ensemble. After this pre-equilibration procedure, larger boxes were generated by duplicating the initial box size in all xyz axis. The other systems (S3 and S4) were created by removing CIT (S3) or by adding ARG (S4), with the number of water molecules then adjusted to the desired ratio. Then, all systems were again simulated for 100 ps (NVT ensemble for temperature equilibration), followed by an 1-µs NpT ensemble run. Additional systems as ethambutol in water (EMB-H2O), citric acid in water (CIT-H2O) or arginine in water (ARG-H2O) were also built for comparison purposes and simulated in similar conditions as described above.

All MD simulations were performed using GROMACS 2020.3 (Pronk et al., 2013; Abraham et al., 2015). Molecule insertion and larger systems were obtained using internal GROMACS routines. While energy minimization steps used the steepest gradient method, all NVT ensemble runs for temperature equilibration used the velocity-rescale (V-rescale) thermostat (Bussi et al., 2007) and all NpT runs the Nosé-Hoover thermostat (310 K) and Parrinello-Rahman barostat (1 bar) (Nosé and Klein, 1983; Hoover, 1985; Nosé, 2002). In NVT ensemble runs, initial velocities were assigned from a Maxwell-Boltzmann distribution and periodic boundary conditions (PBC) were used in all simulations. To allow the calculation of non-bonded interactions in GPUs, a Verlet cut-off scheme was employed (Páll and Hess, 2013). All bond lengths were constrained using LINCS (Hess et al., 1997; Hess, 2008) (SETTLE for TIP5P waters) (Miyamoto and Kollman, 1992). All analysis were performed using the last 500-ns of the NpT ensemble MD simulation runs. Radial distribution functions (RDF) and spatial distribution functions (SDF) were calculated using GROMACS gmx rdf module and Volmap plugin in VMD (Humphrey et al., 1996), respectively. Images were rendered in VMD.

Multiple ligand simultaneous docking (MLSD) (Li and Li, 2010) was performed to evaluate the interaction of the different system’s components with a recently published mycobacterial ABC transporter, found to mediate hydrophilic compounds uptake (PDB ID: 6TQE) (Rempel et al., 2020). Standard molecular docking calculations were performed using a genetic algorithm (GA), in which the population size (ga_pop_size), the maximum number of energy evaluations (ga_num_eval) and the maximum number of generations (ga_num_generations) were set to 150, 5.000.000 and 25.000, respectively. The docking box, with a size of 48 × 68 × 62, grid spacing of 0.375 Å and centered at coordinates 119.15 × 118.874 × 112.419 (xyz coordinates) was defined to include a lower chamber separated from the cytosol by an arginine gate comprising a single Arg287 residues from both protomers. A total of 100 runs (ga_runs) were performed. Single docking using Autodock 4.2.6 was performed for comparison purposes, with ga_pop_size, ga_num_eval and ga_num_generation set to 150, 2.500.000 and 25.000, respectively. All other docking parameters were kept at default. Docking results clustering was performed using AutoDockTools (Morris et al., 2009) summarize_results4. py with a root-mean square tolerance of 10.

Results and Discussion

Molecular Dynamics Studies

Molecular dynamics is an unique tool that allows probing the internal structure of these therapeutic liquid mixtures and to gain insights on the interactions and configurations of each component at its molecular level (Kumari et al., 2018; Cao et al., 2019; Gutiérrez et al., 2019; Sapir and Harries, 2020). Herein, we employed MD simulations to study the molecular structure of each mixture, in order to gain insights on which chemical groups are involved in the formation of the eutectic mixture and to identify suitable chemical modifications that can further improve the properties of ethambutol, enabling them to be used as therapeutic delivery systems (Santos et al., 2018).

Hydrogen Bonding

An initial set of MD simulations were performed comprising each of the liquid mixture component in water (Figure 2), namely EMB/H2O, CIT/H2O, and ARG/H2O. While no aggregation could be observed in either CIT/H2O or ARG/H2O systems, the system EMB/H2O remained inhomogeneous, with EMB molecules aggregated mostly due to intermolecular hydrogen-bonds (HBs) between the amino (pKa1 = 6.35 and pKa2 = 9.35) (Beggs and Andrews, 1974) and hydroxyl groups. Such intermolecular HBs account for nearly 12% of all calculated HBs, with the remaining 88% corresponding to the interactions between the second hydroxyl group and the surrounding water molecules (Figure 3). In the ARG/H2O systems a similar trend is observed (8% ARG-ARG vs. 92% ARG-H2O) however, in the CIT/H2O system almost all calculated HBs are between citric acid and water (99 vs. 1% CIT-CIT interactions).


FIGURE 2. Final configuration (after 1 µs MD simulation time) for all studied systems.


FIGURE 3. Comparison on the average number of hydrogen-bonds between all liquid mixtures components in the studied systems.

When considering a liquid mixture as CIT/ARG/H2O (S1, Figure 3), most of all calculated HBs were observed between ARG/H2O and CIT/H2O pairs (46 and 37%, respectively). Quite interestingly, while intermolecular ARG-ARG hydrogen-bonding increases (10%), all novel ARG-CIT interactions between both components were accomplished by directly replacing water molecules of the first solvation shell by ARG. It seems that, even in a low percentage (6%), the newly established HBs are sufficient to improve ARG solubility and to drive S1 into a more homogeneous state when compared with ARG/H2O alone. For EMB-containing systems S2 and S3, again we observe that most of EMB/H2O HBs were replaced by the more favorable CIT/H2O interactions (71 and 51% of all calculated HBs for S2 and S3, respectively). In both systems, the number of HBs between EMB/CIT was low (2 and 6% for systems S2 and S3, respectively) but sufficient to improve the system’s homogeneity, as it can be seen by the absence of EMB aggregates in both S2 and S3 when compared with EMB/H2O mixture (Figure 3). Finally, in S4 the introduction of ARG replaces EMB/H2O interactions in a larger extent than of CIT/H2O, reaching similar values than those calculated for the S1 system regarding ARG/ARG, CIT/CIT, and ARG/CIT interactions (Figure 3). Therefore, the results seems to be in agreement with previous studies in which the increased solvation of CIT improves the solubility of both ARG and EMB (Santos et al., 2018). In the following section, we describe which chemical groups are responsible for the intermolecular interactions and which maintain their interactions with the surrounding environment.

Radial Distribution Functions of Therapeutic Liquid Mixtures Components

In order to probe the liquid mixtures structure of the different systems and observe the affinity between molecules that compose the mixtures, we calculated RDFs for both center-of mass (c-RDF, Figure 4) and atomic (a-RDF, Figures 5, 6) for all components of the mixture. Regarding c-RDF of systems S2 and S3, while changes in the components ratios only had minor alterations, the introduction of ARG (S4) was revealed by 1) the appearance of three peaks in the RDF profile at r = 0.38 nm (enabling intermolecular HB with EMB), r = 0.56 nm (HB with CIT) and r = 0.82 nm (also with CIT) and 2) by inducing a slight decrease in the water peak (when compared with S2) while reinforcing both citrates peaks at r = 0.54 nm and r = 0.85 nm. Therefore, the data suggests that by adding ARG, the structure of the liquid mixture is reinforced due to its ability to establish intermolecular interactions with both EMB and CIT molecules.


FIGURE 4. Center-of-mass radial distribution functions (c-RDFs) for citric acid (CIT), arginine (ARG), and water (H2O) around ethambutol (EMB) for systems S2-S4. Dashed grey line corresponds to the 0.35 nm cut-off for hydrogen bond calculations.


FIGURE 5. Atomic radial distribution functions (a-RDFs) between citric acid (CIT) and water (SOL) and ethambutol (EMB) for systems S2 and S3, in reference to ethambutol hydrogen-bond acceptor and donor groups. Dashed grey line corresponds to the 0.35 nm cut-off for hydrogen bond calculations.


FIGURE 6. Atomic radial distribution functions for system S4, between ethambutol amine (A), hydrogen-bond acceptor; (B), hydrogen-bond donor and hydroxyl (C), hydrogen-bond acceptor moieties; (D), hydrogen-bond donor groups). Dashed grey line corresponds to the 0.35 nm cut-off for hydrogen bond calculations.

Next, we used a-RDF to better understand which specific interactions could be responsible for the interactions between the various liquid mixture components. From Figure 5 is possible to verify that, concerning direct interactions between CIT and EMB, the interactions were established between the hydrogen from carboxylic acid group at position C-2 (HBD) and EMB nitrogen atoms (HBA), meaning that R2H-N ··· H-O-COCIT in S2 system (CIT/EMB/H2O, 2:1:10; Figure 5) are two times more frequent than with any of the terminal carboxylic acid groups. Thus, the latter will remain available to mediate interactions with the surrounding aqueous environment. However, as titration studies show that EMB has two dissociation constants (Beggs and Andrews, 1974), we also expected that only one amine group undergoes hydrogen-bonding with CIT in the former. This does not change in S3 (CIT/EMB/H2O, 1:1:5), nonetheless the change in the liquid mixture components ratio promoted additional interactions, via hydrogen-bonding, between the carbonyl moiety of carboxylic acid at C-2 and the hydrogen atoms bound to the nitrogen atoms (R2N-H ··· O=C-OH), between the amine moiety and water molecules (R2N-H ··· OH2) and between the hydroxyl groups of CIT and EMB (R2N-H ··· OH-R). Regarding the hydroxyl groups of EMB, in both systems only interactions with water molecules were observed, either as hydrogen-bond donor or acceptor.

Regarding S4 (CIT/EMB/ARG/H2O, 2:1:1:7), the addition of ARG to form a new system fundamentally changed the interaction modes of EMB with the surrounding environment but not with CIT (Figure 6). While the interaction of CIT and EMB was maintained and even reinforced, as it can be seen from the sharper peaks at 0.19 nm in Figure 6A, additional interactions between the amine group of EMB (HBA) and the carboxylic acid (HBD) moiety of ARG (R2H-N ··· H-O-COARG) were observed. Other pairs, although less frequent than the former, could also be observed in which the EMB amine group acted as a HBD (Figure 6B). Concerning EMB hydroxyl groups, they act preferentially as HBDs (Figures 6C vs. 6D). Quite remarkably, ARG appears to be able to disrupt the HB network of water molecules observed in all former systems, as it can be seen by the appearance of a sharp and very intense peak slightly shifted from that of water (r = 0.17 vs. 0.19 nm), assigned to the hydrogen bonding between the hydroxyl and the guanidinium nitrogen N08, and two additional peaks at r = 0.32 nm (guanidinium nitrogen N07) and r = 0.36 nm (α-amino group of ARG).

The results obtained with MD simulations allowed us to understand that 1) only one carboxylic group of CIT is essential to promote interactions with the amine groups of EMB, and that 2) ARG, although also interacting with the same amine groups, prefers to establish intermolecular HBs with EMB hydroxyl groups. Therefore, it is conceivable that both CIT and ARG provide a greater balance between hydrophobicity and hydrophilicity when interacting with EMB, while the remaining functional groups continue immersed in the aqueous environment structuring the therapeutic liquid mixture.

Spatial Distribution Functions

The interaction of the different liquid mixture components can also be observed from the corresponding SDF. Herein, SDF was used to determine specific interactions between the molecules, using a three-dimensional density distribution of CIT, ARG, and water molecules in a local coordinate system linked with ethambutol. From Figure 7A, it is possible to see that CIT mainly occupies a toroidal section around the amine groups of EMB, which is in agreement with results depicted in the RDF plots. Hence, increased water densities are only found in the vicinity of the hydroxyl groups of EMB, interacting with both hydroxyl groups of EMB and with the terminal carboxylic acid groups of CIT, thus assisting in the solubilization of EMB. In Figures 7A,B similar density next to the amino groups of EMB could be assigned to CIT molecules (green arrow), with ARG being found in the close vicinity of EMB hydroxyl groups (yellow arrows). Water, although also interacting with the hydroxyl groups of EMB, is more frequently found in the first hydration shell surrounding ARG or CIT. Thus, data suggests that the presence of ARG confers additional stability to the internal liquid mixture structure by engaging intermolecular interactions with several EMB and CIT molecules, as it can be seen by the identification of several clusters comprising 2–4 EMB molecules during the MD trajectory (data not shown, exemplified in Figure 7B). However, as CIT and ARG are highly functionalized molecules, the remaining groups are kept immersed in the aqueous environment and assist in the effective solubilization of the liquid mixture components.


FIGURE 7. Spatial distribution functions (SDFs) of the liquid mixtures components around ethambutol for (A) S2 and S3, and (B) S4 systems. Specific interactions between citric acid-ethambutol and arginine-ethambutol are identified in the latter by green and yellow arrows, respectively.

Molecular docking

The results presented in the above section provided valuable insights on the internal structure of the therapeutic liquid mixtures and which specific interactions between EMB, CIT, and ARG contribute to the formation of the eutectic system. However, another critical aspect concerning the therapeutic value of this liquid mixture resides on understanding how such systems may promote the incorporation of EMB into the mycobacterial cell.

It is known that the lipidic content in the mycobacterial cell wall, comprising both outer and inner membranes, is highly enriched in diacyl phosphatidylmannosides as Ac2PIM2 or Ac2PIM6, providing an efficient permeability barrier due to its low fluidity and slowing the passive influx of drugs (Sartain et al., 2011; Bansal-Mutalik and Nikaido, 2014). This is a particular question in which concerns EMB, that is expected to be positively charged at physiological pH. To that matter the recent publication of a mycobacterial ABC transporter (Rv1819c, PDB ID: 6TQE), described to be involved in the uptake of hydrophilic compounds through the inner mycobacterial membrane (Rempel et al., 2020), may provide an entry pathway for EMB which is worth to be explored.

In this study, it was explored both single and multiple simultaneous ligands docking approaches to understand 1) if EMB can be internalized via Rv1819c, and 2) if any of the therapeutic liquid mixture components may enhance the interaction of EMB with the transporter. Due to the large internal cavity (∼7.700 Å³), we narrowed our search space to its lower section, defining a docking box ranging from a proposed intracellular gate formed by two R287 residues up to a 17-aminoacid loop in both transmembrane helix three that protrudes into the cavity (Figure 8A).


FIGURE 8. Molecular docking of therapeutic liquid mixtures components. (A) Final configuration for the 100 docking poses from multiple ligand simultaneous docking (MLSD) procedure, with gate residues (Arg287) depicted as volume spheres; (B) top-ranked docking poses for citric acid (green), arginine (yellow) and ethambutol (red) obtained with standard molecular docking; (C), centers-of-mass for the top 20 poses from MLSD procedure for citric acid (green) and ethambutol (red), with gate residues (Arg287) depicted as licorice; and (D) top-ranked docking poses for citric acid (green) and ethambutol (red) from the MLSD procedure, with hydrogen-bonds depicted in cyan.

When EMB, ARG, and CIT are docked separately (Figure 8B), all clusters possess favorable binding energies (Table 2), but while CIT is found to be mostly located in the close vicinity of the R287, both EMB and ARG largest clusters are found in symmetrical locations, interacting with Q190/G201 in one protomer (EMB) or Q190/D193/N204 in the opposite protomer (ARG).


TABLE 2. Clustering results from single molecule molecular docking.

Interestingly, this changes when EMB and CIT are simultaneously docked. By visualizing the center-of-mass for both EMB and CIT obtained from the first 20 docking poses (Figure 8C), the results seem to suggest that most of the CIT molecules are located in the vicinity of the intracellular gate, interacting with R287 (cluster #1, 40 poses, ΔGbind−6.59 kCal.mol⁻1, LE −0.2441) or slightly above (cluster #2, 30 poses, ΔGbind−3.73 kCal.mol⁻1, LE −0.1381), but using its negatively charged carboxylate groups to strengthen the interactions between EMB and Rv1819c. A more detailed analysis of the interactions between the EMB/CIT top-ranked poses and the transporter corroborates this assumption. By establishing additional CIT/EMB intermolecular HBs (Figure 8D), the binding energies of the protein-ligand complex show a difference of only +1.64 kCal.mol⁻1 when compared with EMB and CIT binding energies obtained with the standard docking approach. Moreover, additional interactions with N279 and E280 may also contribute to decrease the binding energy of the protein-ligand complex. When considering the estimated dissociation constants (Kd) this corresponds to a ∼150-fold decrease, from 2.26 mM (EMB) to only 14.85 µM (EMB/CIT complex), which means that internalization of EMB via Rv1819c would be favored in the presence of CIT.

We evaluated additional combinations through MSLD, for comparison purposes. For instance, replacing CIT by ARG slightly improves its binding affinities (cluster #1, 50 poses, ΔGbind −6.79 kCal.mol⁻1, LE −0.2612; cluster #2, 34 poses, ΔGbind−6.04 kCal.mol⁻1, LE −0.2323), but again no molecules were found in the vicinity of R287 and no intermolecular interactions were found within the first 20 poses (except for pose 9). In another combination, including a second CIT, the binding affinity sharply decreases (cluster #1, 48 poses, ΔGbind−2.94 kCal.mol⁻1, LE −0.0735) and no significant intermolecular HBs were observed. Nonetheless, in the presence of an additional ARG molecule lower binding energies are restored and even improved (cluster #1, 21 poses, ΔGbind−8.26 kCal.mol⁻1, LE −0.1588; cluster #2, 14 poses, ΔGbind−7.48 kCal.mol⁻1, LE −0.1438; cluster #3, 32 poses, ΔGbind−5.51 kCal.mol⁻1, LE −0.1060). Again, intermolecular HBs between CIT/EMB/ARG, either in the vicinity of R287 (CIT/ARG) or next to the protruding loops (CIT/EMB) (clusters #1 and #2) or, alternatively, with a single CIT/EMB pair similar to that previously observed (cluster #3). Last but more important, Kd further increases to the nanomolar range (884.65 nM). A ∼15-fold increase was observed when compared with the one reported for the CIT/EMB complex.

Hence, our data suggests that the presence of CIT and/or ARG may also have an important role for EMB uptake by mycobacteria, apart from the prevention of EMB aggregation when present in the liquid mixture. Herein, the exploitation of a possible entry pathway through the ABC importer Rv1819c using single and multiple simultaneous ligand docking suggests that CIT may have a critical role in increasing the binding affinity of EMB towards Rv1819c by mediating specific interactions with its intracellular gate, formed by an ARG pair from both protomers (R287). Although we present a valid route for the uptake of EMB by mycobacteria, other entry routes, e.g. via direct interaction with the lipidic membrane or by means of other undisclosed transporters, cannot be discarded.


Molecular dynamics was applied in this work to understand the internal structure and the molecular interactions that occur with the different components of the therapeutic liquid mixtures combining ARG, EMB, CIT and H2O. It was observed that EMB in water tends to aggregate, but when present in the liquid mixture with CIT and water it is stabilized. The radial distributions functions and spatial distribution functions calculated for different systems demonstrated that different ratios of CIT in the mixtures S2 and S3 could promote additional interactions by hydrogen-bonding in the system S3. In the system S4, with one more component (ARG), the data suggested that the structure of the liquid mixtures was more stable, since stronger intermolecular interactions could be formed due to the ability of ARG to disrupt the water-mediated hydrogen-bond network while promoting the hydrogen-bonding with the other components of the mixture. The MD data of these therapeutic liquid mixtures corroborates that while the carboxyl groups of CIT establish intermolecular interactions with amine groups of EMB, when ARG is present in the liquid mixtures it prefers to interact with hydroxyl groups of EMB and contributing to the solubility of EMB and ARG in water.

Furthermore, from molecular docking, different approaches were used to observe if the internalization of EMB via the mycobacterial ABC transporter (Rv1819c) is propitious by the liquid mixtures, and if the presence of CIT and ARG in the systems increases the affinity to the transporter. According to our findings, the different ratios of CIT in the mixtures did not had a direct impact in the EMB affinity to the binding site of the transporter, but nonetheless an important role for CIT as a mediator for EMB interaction with Rv1819c could also be hypothesized.

This computational approach applied to these therapeutic liquid mixtures represents an important step to understand molecular interactions that can be established between the different components of these low melting transition liquid mixtures, and why they could promote some modifications in properties of active compounds that are present in the mixtures. This will certainly be important to facilitate the entry of the anti-tuberculosis drug EMB in the mycobacterial cell.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author Contributions

All authors have participated and contributed substantially to this manuscript. Design and supervise the study, FS, RF, and AD; performed and analyze the computational studies, HM, and RF; writing the manuscript, FS and RF; review and editing the manuscript, HM, FS, AP, RF, and AD; funding acquisition, AD.


This project has received funding from the European Union’s Horizon 2020 (European Research Council) under grant agreement No ERC-2016-CoG 725034 and was supported by the Associate Laboratory for Green Chemistry-LAQV which is financed by national funds from FCT/MCTES (UIDB/50006/2020).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor declared a shared affiliation with the authors at time of review.


MDR, multidrug resistance; LTTMs, low transition temperature mixtures; HBDs, hydrogen bond donors; HBAs, hydrogen bond acceptors; MD, molecular dynamics; API, active pharmaceutical ingredient; EMB, ethambutol; ARG, arginine.CIT citric acid; PBC, periodic boundary conditions; RDF, radial distribution functions; SDF, spatial distribution functions; MLSD, multiple ligand simultaneous docking; GA, genetic algorithm; HBs, hydrogen bonds.


Abbott, A. P., Ahmed, E. I., Prasad, K., Qader, I. B., and Ryder, K. S. (2017). Liquid Pharmaceuticals Formulation by Eutectic Formation. Fluid Phase Equilibria 448, 2–8. doi:10.1016/j.fluid.2017.05.009

CrossRef Full Text | Google Scholar

Abraham, M. J., Murtola, T., Schulz, R., Páll, S., Smith, J. C., Hess, B., et al. (2015). GROMACS: High Performance Molecular Simulations through Multi-Level Parallelism From Laptops to Supercomputers. SoftwareX 1-2 (2), 19–25. doi:10.1016/j.softx.2015.06.001

CrossRef Full Text | Google Scholar

Aroso, I. M., Craveiro, R., Rocha, Â., Dionísio, M., Barreiros, S., Reis, R. L., et al. (2015). Design of Controlled Release Systems for THEDES-Therapeutic Deep Eutectic Solvents, Using Supercritical Fluid Technology. Int. J. Pharmaceutics 492, 73–79. doi:10.1016/j.ijpharm.2015.06.038

CrossRef Full Text | Google Scholar

Aroso, I. M., Silva, J. C., Mano, F., Ferreira, A. S. D., Dionísio, M., Sá-Nogueira, I., et al. (2016). Dissolution Enhancement of Active Pharmaceutical Ingredients by Therapeutic Deep Eutectic Systems. Eur. J. Pharmaceutics Biopharmaceutics 98, 57–66. doi:10.1016/j.ejpb.2015.11.002

CrossRef Full Text | Google Scholar

Bansal-Mutalik, R., and Nikaido, H. (2014). Mycobacterial Outer Membrane Is a Lipid Bilayer and the Inner Membrane Is Unusually Rich in Diacyl Phosphatidylinositol Dimannosides. Proc. Natl. Acad. Sci. U.S.A. 111, 4958–4963. doi:10.1073/pnas.1403078111

PubMed Abstract | CrossRef Full Text | Google Scholar

Beggs, W. H., and Andrews, F. A. (1974). Chemical Characterization of Ethambutol Binding to Mycobacterium Smegmatis. Antimicrob. Agents Chemother. 5, 234–239. doi:10.1128/AAC.5.3.234

PubMed Abstract | CrossRef Full Text | Google Scholar

Bussi, G., Donadio, D., and Parrinello, M. (2007). Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 126, 014101–014107. doi:10.1063/1.2408420

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, C., Nian, B., Li, Y., Wu, S., and Liu, Y. (2019). Multiple Hydrogen-Bonding Interactions Enhance the Solubility of Starch in Natural Deep Eutectic Solvents: Molecule and Macroscopic Scale Insights. J. Agric. Food Chem. 67, 12366–12373. doi:10.1021/acs.jafc.9b04503

PubMed Abstract | CrossRef Full Text | Google Scholar

Craveiro, R., Aroso, I., Flammia, V., Carvalho, T., Viciosa, M. T., Dionísio, M., et al. (2016). Properties and Thermal Behavior of Natural Deep Eutectic Solvents. J. Mol. Liquids 215, 534–540. doi:10.1016/j.molliq.2016.01.038

CrossRef Full Text | Google Scholar

Dai, Y., van Spronsen, J., Witkamp, G.-J., Verpoorte, R., and Choi, Y. H. (2013). Natural Deep Eutectic Solvents as New Potential Media for Green Technology. Analytica Chim. Acta 766, 61–68. doi:10.1016/j.aca.2012.12.019

CrossRef Full Text | Google Scholar

Dai, Y., Witkamp, G.-J., Verpoorte, R., and Choi, Y. H. (2015). Tailoring Properties of Natural Deep Eutectic Solvents with Water to Facilitate Their Applications. Food Chem. 187, 14–19. doi:10.1016/j.foodchem.2015.03.123

PubMed Abstract | CrossRef Full Text | Google Scholar

Dodda, L. S., Cabeza de Vaca, I., Tirado-Rives, J., and Jorgensen, W. L. (2017a). LigParGen Web Server: An Automatic OPLS-AA Parameter Generator for Organic Ligands. Nucleic Acids Res. 45, W331–W336. doi:10.1093/nar/gkx312

PubMed Abstract | CrossRef Full Text | Google Scholar

Dodda, L. S., Vilseck, J. Z., Tirado-Rives, J., and Jorgensen, W. L. (2017b). 1.14*CM1A-LBCC: Localized Bond-Charge Corrected CM1A Charges for Condensed-phase Simulations. J. Phys. Chem. B 121, 3864–3870. doi:10.1021/acs.jpcb.7b00272

PubMed Abstract | CrossRef Full Text | Google Scholar

Duarte, A. R. C., Ferreira, A. S. D., Barreiros, S., Cabrita, E., Reis, R. L., and Paiva, A. (2017). A Comparison between Pure Active Pharmaceutical Ingredients and Therapeutic Deep Eutectic Solvents: Solubility and Permeability Studies. Eur. J. Pharmaceutics Biopharmaceutics 114, 296–304. doi:10.1016/j.ejpb.2017.02.003

CrossRef Full Text | Google Scholar

Evers, H., Von Dardel, O., Juhlin, L., Ohlsén, L., and Vinnars, E. (1985). Dermal Effects of Compositions Based on the Eutectic Mixture of Lignocaine and Prilocaine (EMLA). Br. J. Anaesth. 57, 997–1005. doi:10.1093/bja/57.10.997

PubMed Abstract | CrossRef Full Text | Google Scholar

Farazi, A., Shafaat, O., Sofian, M., and Kahbazi, M. (2015). Arginine Adjunctive Therapy Tuberculosis. Tuberc. Res. Treat, 2015, 1–5. doi:10.1155/2015/205016

CrossRef Full Text

Fetisov, E. O., Harwood, D. B., Kuo, I.-F. W., Warrag, S. E. E., Kroon, M. C., Peters, C. J., et al. (2018). First-Principles Molecular Dynamics Study of a Deep Eutectic Solvent: Choline Chloride/Urea and its Mixture with Water. J. Phys. Chem. Bacs.jpcb.7b10422 122, 1245. doi:10.1021/acs.jpcb.7b10422

PubMed Abstract | CrossRef Full Text | Google Scholar

Fitzgerald, D. W., Sterling, T. R., and Haas, D. W. (2015). Mycobacterium tuberculosis. Mycobacterium tuberculosis. Mand. Douglas, Bennett’s Princ. Pract. Infect. Dis. 2, 2512787–2512818. doi:10.1016/B978-1-4557-4801-3.00251-4

CrossRef Full Text | Google Scholar

Forbes, B. A., Hall, G. S., Miller, M. B., Novak, S. M., Rowlinson, M.-C., Salfinger, M., et al. (2018). Practical Guidance for Clinical Microbiology Laboratories: Mycobacteria. Clin. Microbiol. Rev. 31, 1–66. doi:10.1128/CMR.00038-17

CrossRef Full Text | Google Scholar

Francisco, M., van den Bruinhorst, A., Zubeir, L. F., Peters, C. J., and Kroon, M. C. (2013). A New Low Transition Temperature Mixture (LTTM) Formed by Choline Chloride+lactic Acid: Characterization as Solvent for CO2 Capture. Fluid Phase Equilibria 340, 77–84. doi:10.1016/j.fluid.2012.12.001

CrossRef Full Text | Google Scholar

Gutiérrez, A., Aparicio, S., and Atilhan, M. (2019). Design of Arginine-Based Therapeutic Deep Eutectic Solvents as Drug Solubilization Vehicles for Active Pharmaceutical Ingredients. Phys. Chem. Chem. Phys. 21, 10621–10634. doi:10.1039/c9cp01408j

PubMed Abstract | CrossRef Full Text | Google Scholar

Gutiérrez, A., Atilhan, M., and Aparicio, S. (2018). A Theoretical Study on Lidocaine Solubility in Deep Eutectic Solvents. Phys. Chem. Chem. Phys. 20, 27464–27473. doi:10.1039/c8cp05641b

PubMed Abstract | CrossRef Full Text | Google Scholar

Hameed, H. M. A., Islam, M. M., Chhotaray, C., Wang, C., Liu, Y., Tan, Y., et al. (2018). Molecular Targets Related Drug Resistance Mechanisms in MDR-, XDR-, and TDR-Mycobacterium tuberculosis Strains. Front. Cel. Infect. Microbiol. 8. doi:10.3389/fcimb.2018.00114

PubMed Abstract | CrossRef Full Text | Google Scholar

Hess, B., Bekker, H., Berendsen, H. J. C., and Fraaije, J. G. E. M. (1997). LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 18, 1463–1472. doi:10.1002/(SICI)1096-987X(199709)18:12<1463::AID-JCC4>3.0.CO;2-H

CrossRef Full Text | Google Scholar

Hess, B. (2008). P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. J. Chem. Theor. Comput. 4, 116–122. doi:10.1021/ct700200b

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoover, W. G. (1985). Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A. 31, 1695–1697. doi:10.1103/PhysRevA.31.1695

PubMed Abstract | CrossRef Full Text | Google Scholar

Humphrey, W., Dalke, A., and Schulten, K. (1996). VMD: Visual Molecular Dynamics. J. Mol. Graphics 14, 33–38. doi:10.1016/0263-7855(96)00018-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Jorgensen, W. L., and Tirado-Rives, J. (2005). Potential Energy Functions for Atomic-Level Simulations of Water and Organic and Biomolecular Systems. Proc. Natl. Acad. Sci. 102, 6665–6670. doi:10.1073/pnas.0408037102

PubMed Abstract | CrossRef Full Text | Google Scholar

Jorgensen, W. L., and Tirado-Rives, J. (1988). The OPLS [optimized Potentials for Liquid Simulations] Potential Functions for Proteins, Energy Minimizations for Crystals of Cyclic Peptides and Crambin. J. Am. Chem. Soc. 110, 1657–1666. doi:10.1021/ja00214a001

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumari, P., Shobhna, S., and Kaur, H. K. (2018). Influence of Hydration on the Structure of Reline Deep Eutectic Solvent: A Molecular Dynamics Study. ACS Omega 3, 15246–15255. doi:10.1021/acsomega.8b02447

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., and Li, C. (2010). Multiple Ligand Simultaneous Docking: Orchestrated Dancing of Ligands in Binding Sites of Protein. J. Comput. Chem. 31, 2014–2022. doi:10.1002/jcc.21486

PubMed Abstract | CrossRef Full Text | Google Scholar

Miyamoto, S., and Kollman, P. A. (1992). Settle: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 13, 952–962. doi:10.1002/jcc.540130805

CrossRef Full Text | Google Scholar

Morris, G. M., Huey, R., Lindstrom, W., Sanner, M. F., Belew, R. K., Goodsell, D. S., et al. (2009). AutoDock4 and AutoDockTools4: Automated Docking with Selective Receptor Flexibility. J. Comput. Chem. 30, 2785–2791. doi:10.1002/jcc.21256

PubMed Abstract | CrossRef Full Text | Google Scholar

Nosé, S., and Klein, M. L. (1983). Constant Pressure Molecular Dynamics for Molecular Systems. Mol. Phys. 50, 1055–1076. doi:10.1080/00268978300102851

CrossRef Full Text | Google Scholar

Nosé, S. U. I. (2002). A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 100, 191–198. doi:10.1080/00268970110089108

CrossRef Full Text | Google Scholar

Páll, S., and Hess, B. (2013). A Flexible Algorithm for Calculating Pair Interactions on SIMD Architectures. Computer Phys. Commun. 184, 2641–2650. doi:10.1016/j.cpc.2013.06.003

CrossRef Full Text | Google Scholar

Pronk, S., Páll, S., Schulz, R., Larsson, P., Bjelkmar, P., Apostolov, R., et al. (2013). GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 29, 845–854. doi:10.1093/bioinformatics/btt055

PubMed Abstract | CrossRef Full Text | Google Scholar

Rempel, S., Gati, C., Nijland, M., Thangaratnarajah, C., Karyolaimos, A., de Gier, J. W., et al. (2020). A Mycobacterial ABC Transporter Mediates the Uptake of Hydrophilic Compounds. Nature 580, 409–412. doi:10.1038/s41586-020-2072-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Rendon, A., Centis, R., D’Ambrosio, L., and Migliori, G. B. (2017). WHO Strategies for the Management of Drug-Resistant Tuberculosis. Archivos de Bronconeumología (English Edition) 53, 95–97. doi:10.1016/j.arbr.2016.07.027

CrossRef Full Text | Google Scholar

Roda, A., Santos, F., Chua, Y. Z., Kumar, A., Do, H. T., Paiva, A., et al. (2021). Unravelling the Nature of Citric Acid:l-Arginine:water Mixtures: the Bifunctional Role of Water. Phys. Chem. Chem. Phys. 23, 1706–1717. doi:10.1039/d0cp04992a

PubMed Abstract | CrossRef Full Text | Google Scholar

Santana, A. P. R., Mora-Vargas, J. A., Guimarães, T. G. S., Amaral, C. D. B., Oliveira, A., and Gonzalez, M. H. (2019). Sustainable Synthesis of Natural Deep Eutectic Solvents (NADES) by Different Methods. J. Mol. Liquids 293, 111452. doi:10.1016/j.molliq.2019.111452

CrossRef Full Text | Google Scholar

Santos, F., P.S. Leitão, M., and C. Duarte, A. (2018). Properties of Therapeutic Deep Eutectic Solvents of L-Arginine and Ethambutol for Tuberculosis Treatment. Molecules 24, 55–12. doi:10.3390/molecules24010055

PubMed Abstract | CrossRef Full Text | Google Scholar

Sapir, L., and Harries, D. (2020). Restructuring a Deep Eutectic Solvent by Water: The Nanostructure of Hydrated Choline Chloride/Urea. J. Chem. Theor. Comput. 16, 3335. doi:10.1021/acs.jctc.0c00120

PubMed Abstract | CrossRef Full Text | Google Scholar

Sartain, M. J., Dick, D. L., Rithner, C. D., Crick, D. C., and Belisle, J. T. (2011). Lipidomic Analyses of Mycobacterium tuberculosis Based on Accurate Mass Measurements and the Novel "Mtb LipidDB". J. Lipid Res. 52, 861–872. doi:10.1194/jlr.M010363

PubMed Abstract | CrossRef Full Text | Google Scholar

Schön, T., Idh, J., Westman, A., Elias, D., Abate, E., Diro, E., et al. (2011). Effects of a Food Supplement Rich in Arginine in Patients with Smear Positive Pulmonary Tuberculosis - A Randomised Trial. Tuberculosis 91, 370–377. doi:10.1016/

PubMed Abstract | CrossRef Full Text | Google Scholar

Silva, C., Perdigao, J., Alverca, E., De Matos, A. P. A., Carvalho, P. A., Portugal, I., et al. (2013). Exploring the Contribution of Mycobacteria Characteristics in Their Interaction with Human Macrophages. Microsc. Microanal. 19, 1159–1169. doi:10.1017/S1431927613001906

PubMed Abstract | CrossRef Full Text | Google Scholar

Sotgiu, G., Sulis, G., and Matteelli, A. (2017). Tuberculosis-a World Health Organization Perspective. Microbiol. Spectr. 5, 1–18. doi:10.1128/microbiolspec.TNMI7-0036-2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Stott, P., Williams, A. C., and Barry, B. W. (1998). Transdermal Delivery from Eutectic Systems: Enhanced Permeation of a Model Drug, Ibuprofen. J. Control Release 50, 297–308. doi:10.1016/S0168-3659(97)00153-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Torfs, E., Piller, T., Cos, P., and Cappoen, D. (2019). Opportunities for Overcoming mycobacterium Tuberculosis Drug Resistance: Emerging Mycobacterial Targets and Host-Directed Therapy. Ijms 20, 2868. doi:10.3390/ijms20122868

PubMed Abstract | CrossRef Full Text | Google Scholar

Viveiros, M., Leandro, C., and Amaral, L. (2003). Mycobacterial Efflux Pumps and Chemotherapeutic Implications. Int. J. Antimicrob. Agents 22, 274–278. doi:10.1016/S0924-8579(03)00208-5

PubMed Abstract | CrossRef Full Text | Google Scholar

WHO (2020). Global Tuberculosis Report 2020. Glob. Tuberc. Rep. 2020.

Google Scholar

WHO (2015). The End TB Strategy. World Heal. Organ.

Keywords: molecular dynamics, therapeutic liquid mixtures, mycobacteria, ethambuthol, tuberculosis

Citation: Monteiro H, Santos F, Paiva A, Duarte ARC and Ferreira RJ (2021) Molecular Dynamics Studies of Therapeutic Liquid Mixtures and Their Binding to Mycobacteria. Front. Pharmacol. 12:626735. doi: 10.3389/fphar.2021.626735

Received: 06 November 2020; Accepted: 06 April 2021;
Published: 20 April 2021.

Edited by:

Diana Machado, New University of Lisbon, Portugal

Reviewed by:

Yun Qian, Shanghai Jiao Tong University, China
Chiara Borsari, University of Basel, Switzerland

Copyright © 2021 Monteiro, Santos, Paiva, Duarte and Ferreira. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Ana Rita C. Duarte,; Ricardo J. Ferreira,