Dialkyl Substituent Effects on the Interaction of Imidazolium Ionic Liquids with POPC Membranes by Molecular Dynamics

1. Introduction

Ionic liquids (ILs) are widely recognized as ionic compounds that remain liquid at relatively low temperatures below 100°C [1,2,3,4]. In addition to their characteristic thermal stability and low volatility, these compounds permit to tune their physicochemical properties by structural modifications in both, cations and anions, as well as by combining diverse cations and anions [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21]. This structural adaptability has facilitated the development of ILs for more sophisticated applications [22]. Whereas early-generation ILs were primarily employed for industrial purposes, nowadays, ILs are increasingly been explored for biological applications [22]. Thefore, it becomes necessary to address significant challenges related with their potentially toxic interactions with living organisms [23,24]. The growing interest in biological applications of ILs has emerged in part as a response to environmental concerns associated with the ecotoxicity of certain IL classes [23,24]. Particular attention has been given to the cytotoxicity of imidazolium-based ILs, which stems from their strong interactions with phospholipid membranes in cellular structures [25,26,27].

The cytotoxic effects of imidazolium ILs appear directly correlated with their membrane-disrupting capabilities. Experimental studies have demonstrated that these cations intercalate into lipid bilayers, inducing perturbations that range from moderate fluidity alterations to complete membrane disintegration, depending on their specific molecular architecture [25,26,27,28,29,30]. A comprehensive understanding of these interaction mechanisms is essential for developing safe pharmaceutical and medical applications of ILs [28,29]. Cytotoxicity assessments across various biological systems, including antifungal and bactericidal tests as well as mammalian cell exposure studies, have revealed that both, cation-anion combinations and alkyl chain lengths, play crucial roles in determining toxicity profiles [1,31,32]. These observed structure-activity relationships highlight the need for more detailed mechanistic studies on IL-membrane interactions.

Recent advances in computational chemistry have provided valuable theoretical insights into the interactions between imidazolium ILs and phospholipid bilayers [33,34,35]. Current computational models primarily focus on two aspects: the structural modifications induced in the bilayer and the molecular mechanisms governing cation insertion and stabilization within membranes. The stabilization of incorporated cations appears to involve both electrostatic interactions with phosphate groups and hydrophobic effects mediated by alkyl side chains [33,34,35]. However, these studies describe in general dialkylated cations with a fixed methyl group [36]. Dialkyl-substituted variants remain relatively unexplored. The presence of two longer alkyl chains may significantly alter both, the insertion dynamics and the membrane perturbation patterns, due to additional steric constraints in the polar headgroup region.

Herein, we employ atomistic molecular dynamics simulations to systematically investigate the interaction of 1-hexadecyl-3-nalkyl-imidazolium cations (where n represents the number saturated carbons in the second alkyl chain) with 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayers, focusing particularly on the structural perturbations induced in the membrane. These simulations provide atomic-level insights into IL-membrane interactions, serving as a tool for elucidating the fundamental molecular mechanisms underlying in IL cytotoxicity.

2. Materials and Methods

We performed molecular dynamics simulations on systems containing a hydrated POPC phospholipid bilayer composed of 128 lipids (64 per leaflet) incorporating 1-hexadecyl-3-nalkyl-imidazolium cations using the GROMACS software package [37]. Each simulated system contained the solvated bilayer along with four ion pairs, with the cationic component consisting of [C₁₆MIM] (1-hexadecyl-3-methylimidazolium), [C₁₆BMIM] (1-hexadecyl-3-butylimidazolium), or [C₁₆OMIM] (1-hexadecyl-3-octylimidazolium), and the chloride (Cl⁻) anion. The molecular structures of the utilized cations are presented in Figure 1. The pre-equilibrated POPC bilayer system has been obtained from the Slipids database [38], with the ionic liquid structures and force fields generated following the AMBER protocol methodology [39,40]. This approach ensured proper system preparation while maintaining consistency with established computational chemistry standards for membrane simulations.

The simulations were conducted in the isothermal-isobaric (NPT) ensemble. Temperature was maintained at 303 K using the velocity-rescaling algorithm [41]. The pressure was regulated at 1 bar by the Berendsen barostat [42] operating in semi-isotropic mode. Three-dimensional periodic boundary conditions were implemented with a 1.4 nm cutoff radius for intermolecular interactions. Long-range electrostatic interactions were calculated using the Particle Mesh Ewald (PME) method [43] with an identical 1.4 nm cutoff distance. The LINCS algorithm [44] was employed to maintain all chemical bond constraints throughout the simulations.

The AMBER force field was selected for molecular representations in the simulations with water molecules modeled described by the TIP3P potential [45]. Ionic pairs were randomly inserted into the hydrated bilayer simulation boxes. Prior to production runs, all systems underwent potential energy minimization using the steepest descent algorithm [46]. Subsequently, simulations were carried out for 100 ns of production time with four cations and chlorides inserted into the bilayer systems. A reference system consisting solely of a hydrated phospholipid bilayer without ionic liquids was subjected to identical simulation protocols and analysis procedures for comparative purposes.

3. Results and Discussions

Along the simulations of the bilayer systems containing Ils, we observed in any case the insertion of the cations from the water phase into the phospholipid bilayer. Afterwards, we extended the simulations up to 100 ns monitoring the convergence of the area per phospholipid (APL) and interactions energies between cations and bilayer molecules. The final 10 ns of each trajectory have been utilized to compute the results presented in the following. A snapshot of the final configuration for each system is presented in Figure 2.

At a glance, we observe the configurations of cations inserted into the bilayer. The [C₁₆OMIM] cation inserts both alkyl chains into the hydrophobic environment of the bilayer. The [C₁₆BMIM] cation, containing the shorter butyl chain (4 carbons), can not fully insert the second substituent into the hydrophobic domain, leaving it closer to the polar headgroups at the membrane surface. As exptected, [C₁₆MIM] directs only the long alkyl chain towards the membrane’s center.

Radial pair distribution functions (RDFs) for distances between the center of the imidazolium ring and three distinct atoms of POPC’s polar head group have been computed. We have chosen the nitrogen atom of the choline group (N), the phosphorous (P), and the innermost oxygen (O) of the glycerol backbone to monitor changes in the coordination of the imidazolium ring within the membrane’s polar region. The RDFs are illustrated in Figure 3.

In general, peak positions of these RDFs are maintained. Both, the imidazolium ring and choline group contain positive charge distributions. Thus, not surprisingly, the RDF with choline’s nitrogen presents maximum amplitudes at larger distances than the other RDFs and appear only slightly affected by the second alkyl substitution in the cations. The [C₁₆MIM] cation presents the most intense peak with the negatively charged phosphate group. Increasing the second alkyl chain, the amplitude of this peak is decreased accompanied by larger amplitudes in the RDF with the oxygen of POPC’s glycerol unit. This observation indicates that longer second alkyl chains favor a deaper insertion of the cation into the membrane.

Figure 4 presents RDFs between water molecules and the aforementioned headgroup atoms of the bilayer. Comparing these functions demonstrates no significant changes in the hydration patterns of the membrane surface. These results indicate that the cation insertion does not promote water intrusion into the bilayer’s hydrophobic domain, preserving the membrane’s fundamental barrier properties despite the structural perturbations caused by the ionic liquids.

The interaction energy between the cations and the phospholipid molecules has been separated into contributions stemming from Lennard-Jones and Coulomb interactions as depicted in Figure 5. More negative values in this analysis indicate stronger contributions to the overall interaction energy within the system. The total interaction energy is increased by adding more atoms to the cations as one might expect. The electrostatic contributions are mostly due to interactions of the cation’s ring and the head group atoms of the membrane and are almost the same for the three cations. On the other side, increasing the second alkyl chain in the cations comes along with enhanced Lennard-Jones interactions. A large fraction of these van der Waals type interactions reflects the hydrophobic interactions between the alkyl chains of the cations and the lipid tails and, therefore, grows strongly with larger alkyl groups in the cations.

In Table 1, we have summarized several structural data for the bilayer in comparison with the unperturbed POPC bilayer. We focused on the area per phospholipid (APL) and the membrane’s thickness computed with the APL@Voro post-processing software [48].

The APL and thickness of the unperturbed bilayer are in agreement with literature data [49]. Our POPC membrane contains 128 phospholipid molecules and four cations which represents a molar fraction of approximately 0.03. We observed slightly increased APLs due to the presence of the cations without a clear tendency. The computed differences in the APLs are smaller than the statistical uncertainty of the numerical values. This holds also for the membrane’s thickness exhibiting numerical values within the error range for the unperturbed bilayer. As outlined by Figure 2, in the case of the [C₁₆BMIM] cation, all the four cations inserted into the same leaflet of the membrane, whereas in the other systems, each leaflet contains two cations causing the apparentely different trends in the APLs and thicknesses. The membrane structural results obtained for the C16MIM cation are in good agreement with values reported in the literature [36] for studies involving smaller monoalkylated cations such as C12MIM, showing that the new data generated for dialkylated cations have good reliability.

In Figure 6, we illustrate the deuterium order parameter for the unsaturated (left panel) and saturated (right panel) alkyl chains of the POPC molecules calculated with the gmx order module of the GROMACS software. The presence of the cations quite generally reduces slightly the order parameter affecting more the carbons close to the membrane’s head group than those close to the membrane’s center

4. Conclusions

Ionic liquids based on the imidazolium cation containing a fixed hexadecyl substituent and a substituent ranging from methyl, butyl and octyl with chloride as the counterion were inserted into hydrated POPC-type phospholipid bilayers. MD simulations of 100 nanoseconds were performed to evaluate the structural effects on the bilayer promoted by the presence of the cations. The stabilization of the cations in the membrane occurs through nonpolar interactions within the bilayer between carbon chains of the cations and lipids and through polar interactions on the membrane surface between the head groups of the phospholipids and the rings present in the cations. Increasing the second alkyl chain in the cations favor configurations with deeper inserted cations. Four inserted cations do not promote significant changes in the overall structure of the membrane’s (area per lipid and thickness), but decrease the deuterium order parameter for POPC chain carbons close to the head groups. In any case, the simulated systems describe stable bilayers without changes in the hydration patterns.