Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
Soft Matter View Article Online
PAPER
Cite this: Soft Matter, 2014, 10, 8641
View Journal | View Issue
Amphiphilic interactions of ionic liquids with lipid biomembranes: a molecular simulation study† Brian Yoo,a Jindal K. Shah,ab Yingxi Zhua and Edward J. Maginn*a Current bottlenecks in the large-scale commercial use of many ionic liquids (ILs) include their high costs, low biodegradability, and often unknown toxicities. As a proactive effort to better understand the molecular mechanisms of ionic liquid toxicities, the work herein presents a comprehensive molecular simulation
study
on
the
interactions
of
1-n-alkyl-3-methylimidazolium-based
ILs
with
a
phosphatidylcholine (PC) lipid bilayer. We explore the effects of increasing alkyl chain length (n ¼ 4, 8, and 12) in the cation and anion hydrophobicity on the interactions with the lipid bilayer. Bulk atomistic molecular dynamics (MD) simulations performed at millimolar (mM) IL concentrations show spontaneous insertion of cations into the lipid bilayer regardless of the alkyl chain length and a favorable orientational preference once a cation is inserted. Cations also exhibit the ability to “flip” inside the lipid bilayer (as is common for amphiphiles) if partially inserted with an unfavorable orientation. Moreover, structural analysis of the lipid bilayer show that cationic insertion induces roughening of the bilayer surface, which may be a precursor to bilayer disruption. To overcome the limitation in the timescale of our simulations, free energies for a single IL cation and anion insertion have been determined based on potential of mean force calculations. These results show a decrease in free energy in response to both short and long alkyl Received 11th July 2014 Accepted 1st September 2014
chain IL cation insertion, and likewise for a single hydrophobic anion insertion, but an increase in free energy for the insertion of a hydrophilic chloride anion. Both bulk MD simulations and free energy
DOI: 10.1039/c4sm01528b www.rsc.org/softmatter
calculations suggest that toxicity mechanisms toward biological systems are likely caused by ILs behaving as ionic surfactants. [Yoo et al., Soft Matter, 2014].
1. Introduction Molecular simulations of synthetic molecules with model biomembranes have become a common method used to explore the physical mechanisms of molecular cytotoxicity. With recent developments in the design of new nanomaterials, there has been increasing concern for the safety of their widespread use, as it has been claimed that such synthetic materials can cause adverse effects toward biological systems.1 Many atomistic and coarse grained molecular simulation studies have provided extensive insight into the interactions that possibly govern the cytotoxic interactions of such materials. In particular, these studies have investigated the interactions between nanomaterials such as carbon nanotubes,2–4 fullerene,2,5–7 and functionalized metal-based nanoparticles8–12 with a hydrated lipid bilayer. Such studies have revealed spontaneous insertion or translocation through the lipid bilayer, provided insight into a
Department of Chemical and Biomolecular Engineering, University of Notre Dame, 182 Fitzpatrick Hall, Notre Dame, Indiana 46556, USA. E-mail: [email protected]; Fax: +1-574-631-8366; Tel: +1-574-631-5687
b
Center for Research Computing, University of Notre Dame, 111 Information Technology Center, Notre Dame, Indiana 46556, USA † Electronic supplementary information (ESI) available: Simulation input les are provided. See DOI: 10.1039/c4sm01528b
This journal is © The Royal Society of Chemistry 2014
the dominating energetic driving forces for molecular insertion or translocation while identifying the subsequent changes in the mechanical properties and structure (i.e. pore formation or rupture) of the lipid bilayer, and determined the likelihood for lipid peroxidation. Among these simulation studies, few have looked into the toxic interactions of a class of recently engineered ionic uids commonly known as ionic liquids (ILs).13–15 Within the past decade, the study of ILs has exploded as these liquids have been touted as promising designer solvents or additives for a wide variety of commercial applications.16–28 ILs are well known for their favorable characteristics: melting points lower than 100 C, good thermal stability, ionic conductivity, high chemical tunability, and negligible vapor pressure. ILs are oen called “green solvents” because their extremely low vapor pressures eliminate the possibility of their release into the atmosphere, unlike conventional volatile solvents. ILs still can nd their way into the environment, however, most likely via aqueous waste streams. Their low volatility makes them odorless, so that leaks are much harder to detect than with conventional solvents. Some ILs have been shown to be toxic toward several types of microorganisms,29,30 mammalian cell lines,31,32 and aquatic organisms33 at concentrations signicantly lower than those of conventional solvents. Thus “green solvent” is somewhat of a
Soft Matter, 2014, 10, 8641–8651 | 8641
View Article Online
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
Soft Matter
misnomer if it is applied to all ILs. At present, the toxicity mechanisms by which ILs cause such adversities toward biological systems are not well understood, in which case molecular simulations can be quite useful. Along with relieving environmental and safety concerns, better understanding of the molecular toxicity mechanisms of ionic liquids provides substantial insight in alternative energy research. Specically, in recent research related to enhanced biofuels production, ILs have shown signicant advantages in the pretreatment of lignocellulosic biomass in comparison to other conventional pretreatment solvent candidates because of their ability to overcome biomass recalcitrance at ambient conditions.17,18 Discoveries made recently by Ruegg et al.34 have shown that bacteria can be engineered to become more IL tolerant, a desirable attribute (for economic reasons) in the development of a “one-pot” process35,36 for fermentation-based biofuels production. For such a case, clear understanding of the molecular toxicity mechanisms of ILs with a cellular membrane could provide enhancements in the biological design choices for the bacterial membrane or in the molecular design of ILs. Recently, Bingham and Ballone13 used classical molecular dynamics (MD) simulations to explore the cytotoxic interactions between solutions of imidazolium-based ILs containing a 1-nbutyl-3-methylimidazolium ([C4mim]+) cation with a zwitterionic 1-palmiltoyl-2-oleoylphosphatidylcholine (POPC) lipid bilayer, a major component in the outer lipid bilayer leaet of eukaryotic cell membranes.37 Their atomistic simulations (over timescales on the order of 100 nanoseconds) have shown that the dominant interactions are likely due to the cations inserting into the lipid bilayer. Hydrophilic anions such as [Cl] and [PF6] tended to freely disperse in the aqueous environment, while in the case of ILs containing a very hydrophobic anion such as bis(triuoromethanesulfonyl)imide ([NTf2]), both the cation and the anion could insert into the lipid bilayer. Such attributes could provide an explanation for the increased toxicities observed for the hydrophobic imidazolium-based ionic liquids toward biological systems. The authors further observed that some of the commonly determined structural properties for the lipid bilayer systems such as the area and volume per lipid were not signicantly altered aer spontaneous insertion of the cation or IL pair. They also did not observe signicant morphology changes or disruption of the lipid bilayer over the length of their simulations. Another recent study by Kl¨ ahn and Zacharias15 has compared the interactions of 1-n-octyl-3-methylimidazolium chloride ([C8mim][Cl]) with model eukaryotic and cancerous cell membranes as means to understand the behavior of ionic liquids as antimicrobial agents. They determined the potential of mean force (PMF) or the free energy for insertion of a single [C8mim]+ into the biomembranes. The membrane models in their work were comprised of varying compositions of phosphatidylcholine, phosphatidylserine (negatively charged lipid), and cholesterol. Kl¨ ahn and Zacharias found that the free energy change for insertion of the [C8mim]+ cation into these membranes is negative and thus thermodynamically favorable. The presence of charged lipids and cholesterol in the bilayer could alter the free energy for the cation insertion. Their
8642 | Soft Matter, 2014, 10, 8641–8651
Paper
calculated free energies showed good agreement with those based on the octanol–water partition constant measured experimentally.38 As a comprehensive simulation study, in this work we have investigated the effect of cation alkyl chain length on the interactions of ILs of varied concentrations with a model lipid bilayer. We report the results of classical MD simulations of aqueous 1-n-alkyl-3-methylimidazolium-based ionic liquids of alkyl chain with n ¼ 4, 8 and 12 in contact with a POPC lipid bilayer. We also examine the effects associated with IL concentration and the hydrophobicity of the anion. The present work is motivated by recent experimental studies where several biological toxicity assays for imidazolium-based ILs have shown that increasing the length of the cation alkyl chain leads to signicantly higher toxicities, while the impact of varying the anion is more modest.30,31,39–41 In an independent study, we have conducted uorescence microscopy experiments on synthetic aPC bilayers in IL solutions and have observed that morphological changes in the bilayer depend on the concentration of the IL, the length of the alkyl chain and the type of anion paired with the cation. These results will be presented in a forthcoming paper. Our simulations are carried out at the nominal concentrations of the microscopy studies. We have also used umbrella sampling42 and the weighted histogram analysis method (WHAM)43 method to compute the PMF44 or free energy change (used interchangeably hereaer) for the insertion of the [C4mim]+ and 1-n-dodecyl-3-methylimidazolium ([C12mim]+) cations as well as the [Cl] and [NTf2] anions into the bilayer at low IL concentration.
2. 2.1
Simulation methods Classical MD simulation
Atomistic molecular dynamics simulations were performed for aqueous solutions of 1-n-alkyl-3-methylimidazolium ionic liquids with a POPC planar lipid bilayer patch using Gromacs 4.5.5.45,46 Cation alkyl chains of length n ¼ 4, 8 and 12 were simulated. Counter anions [Cl], [BF4], [PF6], and [NTf2] were selected based on their varying degree of hydrophobicity. Force eld parameters were taken from Liu et al.47 and Zhong et al.,48 and consist of united atom representations for each methyl group along the alkyl chains and imidazolium ring. The POPC force eld is based on the united atom Berger49 lipid model, which utilizes a combination of GROMOS87 and OPLS parameters with partial charges from ab initio quantum mechanical calculations.50 The model also utilizes a Ryckaert– Bellemans dihedral potential for the torsion angles along the sn1 and sn-2 lipid tails. The POPC lipid bilayer consisted of 128 lipids, with 64 lipids on the upper and lower leaet. The starting conguration for a hydrated lipid bilayer was provided elsewhere.51 Fig. 1(a) shows the different molecules and ions simulated, while Fig. 1(b) shows a snapshot of a representative IL/POPC/water system. For nonbonded interactions between the IL, water, and lipids, we used the Lorentz–Berthelot mixing rule. To account for the differences in 1–4 and charge interaction scaling between the ionic liquid (scaling factor ¼ 0.5) and the lipid
This journal is © The Royal Society of Chemistry 2014
View Article Online
Paper
Soft Matter
LX LY NPOPC =2
(1)
Vsys Nw Vw NPOPC
(2)
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
APOPC ¼ VPOPC ¼
Fig. 1 (a) Representation of the molecules and ions used in this study. Cations are in order of increasing alkyl side chain length: [C4mim]+, [C8mim]+, [C12mim]+ (top to bottom); anions are in order of increasing hydrophobicity: [Cl], [BF4], [PF6], and [NTf2] (left to right). The POPC lipid head (orange) and lipid tail (gray) are also shown. (b) The IL/ POPC/water bilayer system (water omitted for clarity). Box dimensions are roughly 6.3 6.3 17.5 nm in the xyz directions, respectively.
(scaling factor ¼ 1.0) models, we used the “half energy double pair-list” method for Gromacs, in which case the ionic liquid and lipid 1–4 nonbonded interactions were scaled by 0.5, while a pair list index for the lipid charge interactions were doubled. Periodic boundary conditions were used and electrostatic interactions were calculated using the particle mesh Ewald (PME) algorithm.52 The starting congurations of the IL/POPC/water systems were created by rst equilibrating the fully hydrated POPC bilayer conguration (128 POPC lipids solvated by 17 000 SPC water molecules) for 100 ns in the semi-isotropic NPT ensemble corresponding to near zero surface tension. A Nos´ e–Hoover53,54 thermostat was used to maintain the temperature at 298 K, which is above the gelation temperature (Tm ¼ 271 K) of POPC. A Parrinello–Rahman55 barostat was used to maintain the pressure at 1 bar. Ionic liquid ions were then placed roughly 2 nm above and below each leaet of the equilibrated POPC– water conguration aer removal of a few water molecules.56 Each system was subsequently run for another 100 ns. To compare the effect of the increasing alkyl side chain length, we ran simulations with [C4mim][Cl], [C8mim][Cl], and [C12mim][Cl]. To examine the effect of anion hydrophobicity, we ran simulations with [Cl], [BF4], [PF6], and [NTf2] paired with the [C4mim]+ cation. Finally, different amounts of cations and anions were used to study the effect of IL concentration. In all cases, concentration were below the experimental solubility limit of [C4mim][NTf2], which among the studied ILs is the least soluble in water.57
2.2
Lipid bilayer analysis
Equilibrium properties for the lipid bilayer were determined based on block averaging the last 20 ns of each independently run MD simulation.58 The area per lipid (APOPC) and volume per lipid (VPOPC) for the hydrated lipid bilayer were calculated via the following expressions
This journal is © The Royal Society of Chemistry 2014
where LX and LY and the box lengths in the X and Y dimensions respectively; NPOPC and Nw are the number of lipids and waters within the system, respectively, and Vw and Vsys are the volume occupied by a single water molecule and the total volume of the system. An independent simulation of pure water was performed to determine Vw using the same simulation conditions as the IL/POPC/water systems. Additionally, the bilayer thickness (h) was determined for each IL/POPC/water system using the GridMAT-MD analysis tool.59 The phosphorous atom in the POPC head group was selected as a reference atom, upon which a grid (150 149 points) was used to determine their locations and to calculate the local thicknesses. Total bilayer thicknesses were taken by averaging data from the last 20 ns every 200 ps (100 frames in total). Subsequently, they were used to set the integration bounds in the density proles when calculating the number of inserted ILs (see Table 1). Lateral self diffusion coefficients (DL) were determined by the Einstein relation: E 1 D 2 2 DL ¼ lim (3) jxðt0 Þ xðt0 þ tÞj þ jyðt0 Þ yðt0 þ tÞj t/N 2dt t0 where d is the dimensionality (d ¼ 2). Artifacts from the lateral center of mass motion of the bilayer leaets were subtracted from the lipid displacements. To determine DL, each of the 100 ns simulation trajectories was split into ve separate 20 ns trajectories. Slopes of the mean square displacement were taken from 2 to 8 ns and averaged from each of the 20 ns trajectories. This method has been shown to provide reasonable estimates of diffusion coefficients for a bilayer when compared to experimental uorescence correlation spectrometry measurements.60–62 Deuterium order parameters, used to determine the liquid ordered Lo or liquid disordered Ld phase of the lipid bilayer, were determined for individual carbon atoms along the sn-1 and sn-2 tails using the second order Legendre polynomial SCD ¼
1 2 3 cos qCD 1 2
(4)
where qCD is the angle between the carbon–deuterium CD bond (hydrogen for the simulations) and the membrane normal. A value of SCD < 0.35 typically indicates the bilayer is in the Ld uid phase.63 To determine the positions of the hydrogens along the saturated united atom carbons of the lipid tails, a tetrahedral geometry between the central and adjacent carbon atoms was assumed. For unsaturated carbons along the sn-2 tail of the POPC lipid, the vector normal to the plane of the unsaturated carbon and adjacent carbons was used to determine the hydrogen position. It should be noted here that the utility in Gromacs 4.5.5 used to determine the order parameter in this work had been patched to x a bug which assumed an incorrect
Soft Matter, 2014, 10, 8641–8651 | 8643
View Article Online
Soft Matter
Paper
Structural and dynamic properties determined for each IL/POPC/water system over the course of the last 20 ns including: the number of IL ion pairs (nIL), number of water molecules nH2O, concentration of IL (cIL), area (APOPC) per lipid, bilayer thickness (h), lipid lateral diffusion coefficient (DL), average order parameters (hSCDi) for the sn-1 and sn-2 tails, and the number of cations inserted (nIL+). (A nIL+ as a non-integer indicates that the IL is partially inserted.)
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
Table 1
System
nIL
nH2O
cIL/mM
APOPC/nm2
h/nm
POPC–water [C4mim][Cl]
— 5 15 30 15 2 4 5 5 5
17 152 16 877 16 877 16 811 16 877 17 014 16 611 16 877 16 877 16 877
— 16 50 100 50 7 13 16 16 16
0.615 0.612 0.609 0.614 0.604 0.604 0.616 0.586 0.598 0.612
3.885 3.907 3.975 4.005 3.988 3.944 3.957 4.076 3.967 3.939
[C8mim][Cl] [C12mim][Cl] [C4mim][BF4] [C4mim][PF6] [C4mim][NTf2]
0.008 0.006 0.004 0.003 0.002 0.002 0.007 0.001 0.006 0.006
0.172 0.198 0.289 0.277 0.195 0.262 0.265 0.196 0.184 0.262
geometry when determining hydrogen positions for the unsaturated carbons. Electrostatic potential proles were based on integrations of the Poisson equation d2 jðzÞ rc ðzÞ ¼ dz2 30
(5)
where j(z) is the electrostatic potential along the z-direction of the lipid bilayer, rc(z) is the charge density and 30 is the vacuum permittivity. In order to remove artifacts arising from the bilayer center of mass motion, a position correctional term was subtracted from the Z distances.64
2.3 Potential of mean force calculation using umbrella sampling To determine the PMFs for inserting a single [C4mim]+, [C12mim]+, [Cl] and [NTf2] into the POPC bilayer, we used the umbrella sampling and WHAM analysis method implemented in Gromacs.65 For each cation or anion of interest, we allowed for its counterion to freely disperse in the aqueous environment, while implementing a restrictive harmonic potential (in the Z direction of the bilayer plane) on one of the terminal atoms in the cation or anion. To be consistent with our bulk simulations, we paired a [Cl] anion with the cations [C4mim]+ or [C12mim]+ to determine the free energy change for cationic insertion. Likewise, we paired a [C4mim]+ with the [Cl] and [NTf2] anions when determining the anion PMFs. Simulation parameters and the number of water molecules in the systems were the same as those of our bulk simulations. A simulation time of 80 ns for each of the 40–50 frames, with 0.1 nm spacing, was used to obtain better convergence and overlap of energies and to allow for substantial structural relaxation time for the lipid bilayer. More details on the simulation parameters for both bulk MD simulations and umbrella sampling-based PMF calculations can be found in the Supplemental Information. All necessary input les to reproduce our simulations are also provided. Visualizations for our systems have been created with VMD.66
8644 | Soft Matter, 2014, 10, 8641–8651
3. 3.1
DL/108 cm2 s1
hSCDisn-1
hSCDisn-2
nIL+ inserted
1.20 0.41 3.62 0.47 3.88 0.71 3.70 0.45 4.63 1.32 1.17 0.19 1.17 0.17 3.90 0.79 4.72 0.66 3.89 0.35
0.184 0.034 0.185 0.032 0.194 0.039 0.197 0.038 0.201 0.042 0.193 0.034 0.186 0.032 0.211 0.038 0.197 0.030 0.196 0.033
0.146 0.049 0.148 0.052 0.153 0.049 0.156 0.051 0.155 0.054 0.154 0.048 0.149 0.050 0.169 0.053 0.151 0.045 0.153 0.050
— 2.0 6.6 9.5 5.0 2.0 2.7 2.0 1.0 4.0
Results and discussion Bulk MD simulation
3.1.1 Fully hydrated POPC bilayer validation. To draw comparison between the lipid bilayer structural properties for the IL/POPC/water MD simulations, a fully hydrated POPC lipid bilayer system was simulated and equilibrium properties were compared with those from the literature.13,60,67–69 The average area per lipid, bilayer thickness, lateral diffusion coefficient, and average order parameters are listed in Table 1. We obtain an average area per lipid of 0.615 nm2 with bilayer thicknesses of 3.885 nm, which is in close agreement with the simulation results of Mukhopadhyay et al.67 The average area per lipid is slightly lower than the experimental value of 0.64 0.01 nm2 at 298 K.69 The volume per lipid and volume per water are 1.205 nm3 and 0.0306 nm3, respectively, are also in close agreement with previous simulation results and are in reasonable agreement with experimental studies (1.256 nm3 at 303 K and 0.030 nm3 at 298 K, respectively).70 Lateral diffusion coefficients of the lipids obtained in this work are roughly a factor of 3 lower than those calculated by B¨ ockmann et al.,60 who obtained a lateral diffusion coefficient of 3.9 108 cm2 s1 at 300 K using the same lipid and water force elds. This disparity is likely due to the differences in our simulation conditions. The van der Waals and electrostatic cutoffs in our study (rcut ¼ 1.2 nm) are slightly longer than those used by B¨ ockmann et al.60 (rcut ¼ 1.0 nm). A longer cutoff was used in the present work in order to be consistent with the cutoff used in parameterizing the IL force eld. Another difference between the present simulations and those of B¨ ockmann et al.60 is that we used a Nos´ e–Hoover thermostat, while B¨ ockmann et al. used a Berendsen thermostat. It is known that a Berendsen thermostat incorrectly samples canonical ensemble uctuations, which in turn can lead to faster dynamics.71,72 3.1.2 IL–lipid bilayer interaction and the effect of IL cation alkyl chain length. Independent simulations at concentrations of roughly 50 mM [C4mim][Cl], 50 mM [C8mim][Cl], and two systems at 7 and 13 mM [C12mim][Cl], were performed to observe the effect of the increased cation alkyl side chain length. For the [C12mim][Cl] systems, we did not run a full 100
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
Paper
ns simulation at a 50 mM concentration (corresponding to 15 IL pairs) because the cations tended to aggregate into micelles. For this system our simulation timescales are too short (even aer 60 ns) to observe these aggregates insert into the lipid bilayer. For the other systems, we observed that a fraction of the cations spontaneously inserted into the lipid bilayer, while all of the [Cl] anions remain freely dispersed in water. Cations that insert into the bilayer during the simulation tend to favor orientations in which the alkyl side chain (“IL tail”) inserts deeply into the hydrophobic region of the lipid bilayer, while the imidazolium ring (“IL head”) interacts with the lipid head groups. (The terms “IL tail” and “IL head” are used to draw the similarities of these ILs with surfactants.) For some of the IL cations, we observe partial insertion of the IL head, and eventual reorientation or “ipping”, such that the IL tail is inserted into the hydrophobic regime of the bilayer, as shown in Fig. 2. Also indicated in the gure are the angles qNC between the membrane leaet normal and the vector between a nitrogen from the imidazolium ring and the terminal methyl group of the alkyl side chain. When qNC is approximately greater than 120 , the tail can be assumed to be completely inserted into the bilayer. For the 50 mM [C4mim][Cl], 50 mM [C8mim][Cl] and 13 mM [C12mim][Cl] systems, we observe at least one partially inserted IL ip, such that the IL tail is inserted deeper
Soft Matter
into the lipid bilayer, while the IL head remains near the lipid and water interface. However, the time it takes for a partially inserted IL cation to ip is signicantly longer than the time it takes for a cation to insert directly in a “tail rst” orientation. Normalized density proles of the orientation-favored insertion of the cations for the 50 mM [C4mim][Cl], 50 mM [C8mim][Cl], and 7 mM [C12mim][Cl] systems are plotted in Fig. 3(a)–(c). The density proles have been averaged based on the last 20 ns of the simulation run. For the longer alkyl chain cations, the density proles show deeper penetration of the alkyl tails into the lipid bilayer with density maxima closer to the lipid bilayer center, while the imidazolium rings stay at the water–lipid interface. The [Cl] anion density proles show no indication of penetration into the bilayer or strong reassociation with the inserted cations. The results of the SCD order parameter calculations are shown in Fig. 4. Cation insertion tends to increase the local ordering of the sn-1 and sn-2 lipid tails. The carbons in the lipid tails near the lipid head groups (cn with n ¼ 2 to 5) exhibit slightly increased ordering compared with those at the end of the tails, as indicated by the increase in SCD for the [C4mim][Cl] system. For the 50 mM [C8mim][Cl] system, insertion of the cations causes a greater increase in SCD than that of [C4mim][Cl], particularly for the sn-1 carbons cn with n ranging from 2 to 12. For the 7 and 13
Fig. 2 Snapshots of partially inserted IL cation “flipping” into the bilayer during the course of the simulation run for systems: (a) 50 mM [C4mim][Cl], (b) 50 mM [C8mim][Cl], and (c) 13 mM [C12mim][Cl]. The angles qNC between the membrane leaflet normal and vector between a nitrogen from the imidazolium ring and the terminal methyl group of the alkyl side chain (NC) are indicated. Arrows pointing in the direction of the NC-vector are shown for guidance. The rest of the anions, cations, and water molecules in each system are omitted for clarity.
This journal is © The Royal Society of Chemistry 2014
Soft Matter, 2014, 10, 8641–8651 | 8645
View Article Online
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
Soft Matter
Paper
Fig. 3 Normalized density r* profiles along the Z axis for lipid bilayer systems in aqueous solutions of (a) 50 mM [C4mim][Cl], (b) 50 mM [C8mim][Cl], (c) 7 mM [C12mim][Cl], (d) 16 mM [C4mim][BF4], (e) 16 mM [C4mim][PF6], and (f) 16 mM [C4mim][NTf2]. Density profiles for the POPC bilayer (black) has been multiplied by 2 for clarity. For almost all cases, cations and anions are freely disassociated. Density profiles for the IL cation head (blue) and tail (cyan) are plotted separately. Those for anions are also plotted (green).
Fig. 4 (a) Deuterium order parameter for the sn-1 (solid) and sn-2 (dashed) tails for the POPC–water system (black), the 50 mM [C4mim] [Cl] (blue) and the 50 mM [C8mim][Cl] (red) system.
mM [C12mim][Cl] systems, we observe a similar trend in the sn-1 and sn-2 order parameters, although a shi in the order parameters is not as noticeable for these systems, due to the lower number of inserted cations. For the [C4mim][Cl] and [C8mim][Cl] systems, lateral diffusion coefficients of the lipids tend to increase by roughly a factor of 2–4 as compared to the lipid without ILs. An increase in the lateral diffusion coefficient and increase in the sn-1 and sn-2 order parameters may suggest that insertion of the IL cations induces a tendency for greater collective lateral motion, in which case the lipid head groups in these systems are possibly conned locally due to the presence of a nearby inserted cation and forced to displace toward a nearby lipid.
8646 | Soft Matter, 2014, 10, 8641–8651
For the [C12mim][Cl] systems, we do not observe a similar increase in the lipid lateral diffusion coefficient. This can be attributed to the deeper insertion of the longer alkyl chain into the bilayer, which may interact much like a lipid tail group. Recently, it has been shown that a large electrostatic potential gradient induced by the presence of ions and cationic nanoparticles can drive pore formation within a double bilayer.11 It should be noted that we observe slight changes to the electrostatic potential aer cation insertion for these systems, although it is difficult to determine if such deviations were directly caused by cation insertion. Previous studies have shown that slight deviations in the electrostatic potential can be caused by uctuations and averaging over an inadequate simulation time (4.2 nm) occur more frequently with an increased number of inserted cations. The positions of the inserted cations (represented by the magenta dots in Fig. 6) do not appear to strongly correlate with the locations of the regions of high or low bilayer thicknesses. Hill et al.75 have shown that bilayer thickness deformation, particularly when the thicknesses decreases signicantly (< 0.8 nm) in comparison to an equilibrated hydrated lipid bilayer, tends to eventually lead to formation of a porous channel in the low thickness regions, allowing for the transport of water or ions across the bilayer. However, the reason for the pore formation in their simulations was due to strong charge interactions between aggregates of the inserted cytotoxic molecules with the charged head groups of the nearby lipids in the bilayer hydrophobic region. We do not observe any pore formation or diffusion of inserted cations into large aggregates for the IL/ water/POPC systems within the simulated timescale, though the locations of the inserted cations in the 100 mM system are positioned relatively close to one another. It is thus suggested that the ILs denitely cause increased surface roughness on 100 ns timescales, and this could be a precursor to bilayer disruption. As was the case with the increasing chain length, the diffusion coefficients are roughly 3 times larger for all the systems when compared to the hydrated lipid bilayer. Interestingly, although the lateral diffusion coefficient is higher for the cases when ILs insert into the bilayer, it does not appear to change as the number of inserted ions increases. Given that a small number of systems were studied, it is not possible to know if this is a general phenomena. Likewise, order parameters increase very slightly with increasing IL concentration, while the average overall bilayer thickness, based on the locations of the phosphorus atoms in the lipid head group, is insensitive to IL concentration.
3.2
Free energy change upon cation or anion insertion
The Gibbs free energy change DG for a single cation or anion insertion from bulk aqueous solution was determined from potential of mean force calculations in the semi-isotropic NPT ensemble. From our observation of the orientational preference of inserted IL cations in our bulk simulations, we decided to select as the reaction coordinate for [C12mim]+ and [C4mim]+
Soft Matter, 2014, 10, 8641–8651 | 8647
View Article Online
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
Soft Matter
Paper
Fig. 6 Snapshots of configurations (above) and bilayer thickness contour plots (h) (below). The contour plots were averaged over the last 20 ns of the simulation. (a) The hydrated lipid bilayer, (b) 16 mM [C4mim][Cl], (c) 50 mM [C4mim][Cl], and (d) 100 mM [C4mim][Cl]. Locations of inserted cation atoms from each frame are indicated by the dots in magenta.
the distance along the Z direction between the center of mass of the lipid bilayer and the terminal methyl group on the alkyl chain, as shown in Fig. 7(a). We found that the selection of this reaction coordinate allows for quicker convergence in the free energy prole due to the relatively long ipping time of the aliphatic [C12mim]+. For [NTf2], one of the united atom CF3 groups was chosen as the reference. The center of mass was chosen as the reference for the other anions. Fig. 7(b) and (c) plots the free energy proles for insertion of a single [C4mim]+, [C12mim]+, [Cl] and [NTf2] from bulk solution. The free energies of both [C4mim]+ and [C12mim]+ (paired with a [Cl]) monotonically decrease from bulk solution as the cations approach the lipid bilayer. The insertion process is essentially barrierless, indicating that the insertion process is thermodynamically spontaneous. The minimum in the free energy for the [C4mim]+ insertion is roughly 27 kJ mol1 at z z 1.2 nm, while for [C12mim]+ the minimum is 37 kJ mol1 at z z 0.07 nm. The location of each free energy minimum along the reaction coordinate corresponds
to the cation head interacting strongly with the lipid head groups and the depth of the cation tail insertion. The sharp increase in the free energy prole starting at about z < 0.5 nm for [C4mim]+ and the slight increase at z < 0.07 nm for [C12mim]+ indicate that the cation head groups prefer to reside at the bilayer interface with the tail group extended into the hydrophobic regime of the lipid bilayer. Thus the cations do not cross the lipid bilayer, but rather remain inserted in the bilayer. Kl¨ ahn and Zacharias have reported a minimum free energy for [C8mim]+ insertion into a POPC bilayer of 15 kJ mol1 at about z z 0.6 nm using a similar reaction coordinate. Their value is roughly a factor of 2 smaller in magnitude than those of our [C4mim]+ and [C12mim]+ estimates. Based on the qualitative trends of our free energy minimum, we would expect our values to be relatively close to those of Kl¨ ahn and Zacharias. The minimum free energy for [C12mim]+ should be less, while that for the [C4mim]+ should be greater. The most likely explanation for this discrepancy is in the difference in sampling time. While Kl¨ ahn and Zacharias obtained a minimum free energy close in
Fig. 7 Free energy profiles DG along the reaction coordinate z. (a) The reaction coordinates z for cation or anion are defined based on the Z distance between the reference atom and the center of mass of the lipid bilayer (z ¼ 0). Free energy profiles for the insertion of a single (b) [C4mim]+ (black) and [C12mim]+ (blue), and (c) [Cl] (green) and [NTf2] (red) anion into a POPC lipid bilayer. Uncertainties are indicated by the shaded regions.
8648 | Soft Matter, 2014, 10, 8641–8651
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
Paper
value to experimental studies, they used a sampling time of 10 ns per conguration. We nd that although we obtain a free energy minimum close in value with theirs within 10–20 ns, the free energies do not show good convergence as they continue to decrease until about 50–80 ns. A longer sampling time per conguration allows for longer structural relaxation or equilibration of the bilayer in response to the inserted IL. Another reason for the discrepancy can arise from the differences in the lipid and IL force elds, as their IL force eld partial charges were reparameterized to account for the solvation effects of water. Nevertheless, the shape in their free energy curve and the location of the minimum free energy obtained are consistent and as expected with those of our work. For the most hydrophobic [NTf2] anion the free energy minimum is approximately 10 kJ mol1 at z z 0.6 nm, while for the hydrophilic [Cl] anion there is only a shallow minimum free energy at the bilayer surface. As the [Cl] anion moves into the bilayer, the free energy is strongly repulsive. The PMF results conrm the results from the bulk MD simulations. There is a larger thermodynamic driving force for the longer alkyl chain cations to insert into the bilayer, though even the short chain cations have negative free energies for insertion. The hydrophobic [NTf2] anion has a moderate tendency to insert into the bilayer, while the hydrophilic [Cl] shows no tendency to independently insert into the bilayer. In addition, we nd that the free energy change in the insertion of [NTf2] is independent of the [C4mim]+ cation insertion. Note that the PMF calculations were all performed on a single ion. The bulk MD simulations suggest that the presence of an inserted cation could lower the free energy barrier for subsequent anion insertion, but this was not examined in the present work.
4. Conclusions Molecular dynamics simulations and potential of mean force calculations were used to study the interactions of a range of different imidazolium-based ionic liquids with a POPC lipid bilayer. The simulations show that the ionic liquid cations tend to spontaneously insert into the lipid bilayer, while exhibiting strong amphiphilic interactions. Specically, we nd that as the cations insert, regardless of the alkyl chain length, cations favor orientation such that the imidazolium ring and the alkyl chain are strongly associated with the lipid head and tail groups, respectively. IL cations that have partially inserted initially with the imidazolium ring in the bilayer demonstrate the ability to ip such that the alkyl chain is completely inserted into the hydrophobic region of the bilayer. This ipping time tends to be much longer for the more aliphatic imidazolium cations of increasing alkyl chain length ([C8mim][Cl] and [C12mim][Cl]). Insertion of the cations tends to change the structural and dynamic properties of the bilayer. In particular, the ionic liquids tend to roughen the surface of the bilayer, with higher concentrations of ionic liquids leading to greater surface roughness. On the timescales of the simulation (100 ns) we do not observe direct disruption or pore formation in the bilayer, but this roughening could be a precursor for bilayer disruption, a potential cause of ionic liquid toxicity.
This journal is © The Royal Society of Chemistry 2014
Soft Matter
The effect of the anion hydrophobicity of the ionic liquid was studied by pairing the same cation ([C4mim]+) with [Cl], [BF4], [PF6], and [NTf2]. The latter three anions are more hydrophobic than [Cl] and tend to have a slightly greater effect on the structural properties of the bilayer than that of [Cl]. Within the timescales of the simulations, the [NTf2] anion is the only hydrophobic anion that inserts into the bilayer and reassociates with an inserted cation. Unlike the cations, the degree of hydrophobicity in the anion does not correlate well with the changes in the structural and dynamic properties of the bilayer. Umbrella sampling potential of mean force calculations were performed to compute the free energies for inserting various ions into the bilayer. The results conrm the ndings of the bulk MD simulations. Insertion of the cations is thermodynamically favorable, where the free energy change for insertion of [C4mim]+ is 27 kJ mol1 and is 37 kJ mol1 for the insertion of [C12mim]+. The [C12mim]+ cation also inserts more deeply into the bilayer than does [C4mim]+. The free energy change for insertion of the hydrophobic [NTf2] anion is also thermodynamically favorable at 10 kJ mol1. The [Cl] anion exhibits a positive free energy as it approaches the lipid bilayer, suggesting that it is more thermodynamically favorable for it to remain hydrated in the aqueous phase or to reside near the charged head groups of the bilayer. This study demonstrates that the deep insertion of imidazolium cations with long alkyl chains into the lipid bilayer is a likely mechanism for greater cytotoxicity toward biological systems. We also conrm insertion of hydrophobic anions with other computational work.13 Although we do not observe direct pore or channel formation, previous simulation studies have shown that the presence of charged molecule aggregates (particularly those that are cationic) or deeply inserted nonionic surfactants into a lipid bilayer can cause such structural damage.11,75,76 The results of this study suggest that the underlying toxicity mechanisms of the ionic liquids are governed by the ionic liquids essentially behaving as ionic surfactants. This classication of ionic liquids as ionic surfactants is by no means a novel discovery, as many studies have shown that ILs with longer alkyl chains (typically for [Cnmim]+ with n > 7) exhibit critical micelle concentrations,77–84 a common characteristic of surfactants.85 Within the context of understanding IL toxicity mechanisms toward biological systems, however, this study may suggest that assessments on the environmental safety concerns for imidazolium-based ionic liquids (including those with shorter alkyl chains such as [C4mim]+) can be made similarly to that of more commonly used ionic surfactants in industry, although the biodegradability of these ILs should still be considered.86 Further studies such as those modeling larger systems or longer simulation times, with either atomistic simulations or coarse grained simulations, need to be performed to validate such a mechanistic basis.
Acknowledgements Funding for this project was provided by the National Science Foundation (CBET-1134238). Computational resources were
Soft Matter, 2014, 10, 8641–8651 | 8649
View Article Online
Soft Matter
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
provided by Notre Dame's Center for Research Computing. We are also grateful for Prof. Peter Tieleman for providing the lipid force eld and bilayer initial conguration les.
References 1 A. Nel, T. Xia, L. M¨ adler and N. Li, Science, 2006, 311, 622– 627. 2 L. Monticelli, E. Salonen, P. C. Ke and I. Vattulainen, So Matter, 2009, 5, 4433–4445. 3 J. Liu and A. J. Hopnger, Chem. Res. Toxicol., 2008, 21, 459– 466. 4 S. Pogodin and V. A. Baulin, ACS Nano, 2010, 4, 5293–5300. 5 R. Qiao, A. P. Roberts, A. S. Mount, S. J. Klaine and P. C. Ke, Nano Lett., 2007, 7, 614–619. 6 J. Wong-Ekkabut, S. Baoukina, W. Triampo, I.-M. Tang, D. P. Tieleman and L. Monticelli, Nat. Nanotechnol., 2008, 3, 363–368. 7 R. DeVane, A. Jusu, W. Shinoda, C.-c. Chiu, S. O. Nielsen, P. B. Moore and M. L. Klein, J. Phys. Chem. B, 2010, 114, 16364–16372. 8 K. Yang and Y.-Q. Ma, Nat. Nanotechnol., 2010, 5, 579–583. 9 V. V. Ginzburg and S. Balijepalli, Nano Lett., 2007, 7, 3716– 3722. 10 R. C. Van Lehn and A. Alexander-Katz, So Matter, 2011, 7, 11392–11404. 11 J. Lin and A. Alexander-Katz, ACS Nano, 2013, 7, 10799– 10808. 12 R. C. Van Lehn and A. Alexander-Katz, So Matter, 2014, 10, 648–658. 13 R. J. Bingham and P. Ballone, J. Phys. Chem. B, 2012, 116, 11205–11216. 14 S. Cromie, M. Del P´ opolo and P. Ballone, J. Phys. Chem. B, 2009, 113, 11642–11648. 15 M. Kl¨ ahn and M. Zacharias, Phys. Chem. Chem. Phys., 2013, 15, 14427–14441. 16 H. Wu and E. J. Maginn, Fluid Phase Equilib., 2014, 368, 72– 79. 17 S. Singh, B. A. Simmons and K. P. Vogel, Biotechnol. Bioeng., 2009, 104, 68–75. 18 M. Francisco, A. N. Mlinar, B. Yoo, A. T. Bell and J. M. Prausnitz, Chem. Eng. J., 2011, 172, 184–190. 19 H. Rodr´ıguez, S. Padmanabhan, G. Poon and J. M. Prausnitz, Bioresour. Technol., 2011, 102, 7946–7952. 20 L. D. McIntosh, T. Kubo and T. P. Lodge, Macromolecules, 2014, 47, 1090–1098. 21 A. S. Pensado, F. Malberg, M. F. C. Gomes, A. A. H. Padua, J. Fernandez and B. Kirchner, RSC Adv., 2014, 4, 18017– 18024. 22 W. Afzal, B. Yoo and J. M. Prausnitz, Ind. Eng. Chem. Res., 2012, 51, 4433–4439. 23 K. A. Maerzke, G. S. Goff, W. H. Runde, W. F. Schneider and E. J. Maginn, J. Phys. Chem. B, 2013, 117, 10852–10868. 24 X. Kang, J. Zhang, W. Shang, T. Wu, P. Zhang, B. Han, Z. Wu, G. Mo and X. Xing, J. Am. Chem. Soc., 2014, 136, 3768–3771. 25 J. Stoimenovski, D. R. MacFarlane, K. Bica and R. D. Rogers, Pharm. Res., 2010, 27, 521–526.
8650 | Soft Matter, 2014, 10, 8641–8651
Paper
26 R. Ferraz, L. C. Branco, C. Prudencio, J. P. Noronha and ˇ. Petrovski, ChemMedChem, 2011, 6, 975–985. Z 27 S. Li, H. Gao and J. M. Shreeve, Angew. Chem., Int. Ed., 2014, 53, 2969–2972. 28 B. Yoo, W. Afzal and J. M. Prausnitz, Ind. Eng. Chem. Res., 2012, 51, 9913–9917. 29 R. P. Swatloski, J. D. Holbrey, S. B. Memon, G. A. Caldwell, K. A. Caldwell and R. D. Rogers, Chem. Commun., 2004, 668–669. 30 D. J. Couling, R. J. Bernot, K. M. Docherty, J. K. Dixon and E. J. Maginn, Green Chem., 2006, 8, 82–90. 31 R. A. Kumar, N. Papa¨ıconomou, J.-M. Lee, J. Salminen, D. S. Clark and J. M. Prausnitz, Environ. Toxicol., 2009, 24, 388–395. 32 J. Ranke, A. M¨ uller, U. Bottin-Weber, F. Stock, S. Stolte, J. Arning, R. St¨ ormann and B. Jastorff, Ecotoxicol. Environ. Saf., 2007, 67, 430–438. 33 A. Latała, P. Stepnowski, M. Nedzi and W. Mrozik, Aquat. Toxicol., 2005, 73, 91–98. 34 T. L. Ruegg, E.-M. Kim, B. A. Simmons, J. D. Keasling, S. W. Singer, T. S. Lee and M. P. Thelen, Nat. Commun., 2014, 5, DOI: 10.1038/ncomms4490. 35 H. Zhao, S. Xia and P. Ma, J. Chem. Technol. Biotechnol., 2005, 80, 1089–1096. 36 P. W. Wolski, D. S. Clark and H. W. Blanch, Green Chem., 2011, 13, 3107–3110. 37 G. Van Meer, D. R. Voelker and G. W. Feigenson, Nat. Rev. Mol. Cell Biol., 2008, 9, 112–124. 38 M. Matzke, S. Stolte, K. Thiele, T. Juffernholz, J. Arning, J. Ranke, U. Welz-Biermann and B. Jastorff, Green Chem., 2007, 9, 1198–1207. 39 P. Stepnowski, A. Skladanowski, A. Ludwiczak and E. Laczy´ nska, Hum. Exp. Toxicol., 2004, 23, 513–517. 40 B. Jastorff, R. St¨ ormann, J. Ranke, K. M¨ olter, F. Stock, B. Oberheitmann, W. Hoffmann, J. Hoffmann, M. N¨ uchter, B. Ondruschka, et al., Green Chem., 2003, 5, 136–142. 41 F. Stock, J. Hoffmann, J. Ranke, R. St¨ ormann, B. Ondruschka and B. Jastorff, Green Chem., 2004, 6, 286–290. 42 B. Roux, Comput. Phys. Commun., 1995, 91, 275–282. 43 S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen and P. A. Kollman, J. Comput. Chem., 1992, 13, 1011–1021. 44 D. Chandler, Introduction to Modern Statistical Mechanics, by David Chandler, pp. 288. Foreword by David Chandler, Sep 1987, ISBN-10: 0195042778. ISBN-13: 9780195042771, Oxford University Press, 1987, vol. 1, pp. 201–202. 45 E. Lindahl, B. Hess and D. van der Spoel, Molecular Modeling Annual, 2001, 7, 306–317. 46 B. Hess, C. Kutzner, D. Van Der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447. 47 Z. Liu, X. Wu and W. Wang, Phys. Chem. Chem. Phys., 2006, 8, 1096–1104. 48 X. Zhong, Z. Liu and D. Cao, J. Phys. Chem. B, 2011, 115, 10027–10040. 49 O. Berger, O. Edholm and F. J¨ ahnig, Biophys. J., 1997, 72, 2002–2013. 50 S.-W. Chiu, M. Clark, V. Balaji, S. Subramaniam, H. L. Scott and E. Jakobsson, Biophys. J., 1995, 69, 1230–1245.
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
Paper
51 http://wcm.ucalgary.ca/tieleman/downloads, accessed 07-012013. 52 T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092. 53 S. Nos´ e, J. Chem. Phys., 1984, 81, 511–519. 54 W. G. Hoover, Phys. Rev. A, 1985, 31, 1695–1697. 55 M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182– 7190. 56 L. Mart´ınez, R. Andrade, E. G. Birgin and J. M. Mart´ınez, J. Comput. Chem., 2009, 30, 2157–2164. 57 M. G. Freire, P. J. Carvalho, R. L. Gardas, I. M. Marrucho, L. M. Santos and J. A. Coutinho, J. Phys. Chem. B, 2008, 112, 1604–1610. 58 M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford university press, 1989. 59 W. J. Allen, J. A. Lemkul and D. R. Bevan, J. Comput. Chem., 2009, 30, 1952–1958. 60 R. A. B¨ ockmann, A. Hac, T. Heimburg and H. Grubm¨ uller, Biophys. J., 2003, 85, 1647–1655. 61 S. W. Siu, R. V´ acha, P. Jungwirth and R. A. B¨ ockmann, J. Chem. Phys., 2008, 128, 125103. 62 A. A. Gurtovenko and I. Vattulainen, J. Phys. Chem. B, 2008, 112, 1953–1962. 63 O. G. Mouritsen and K. Jrgensen, Chem. Phys. Lipids, 1994, 73, 3–25. 64 A. A. Gurtovenko and I. Vattulainen, J. Chem. Phys., 2009, 130, 215107. 65 J. S. Hub, B. L. de Groot and D. van der Spoel, J. Chem. Theory Comput., 2010, 6, 3713–3720. 66 W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38. 67 P. Mukhopadhyay, H. J. Vogel and D. P. Tieleman, Biophys. J., 2004, 86, 337–345. 68 J. Seelig and N. Waespe-Sarcevic, Biochemistry, 1978, 17, 3310–3315.
This journal is © The Royal Society of Chemistry 2014
Soft Matter
69 B. K¨ onig, U. Dietrich and G. Klose, Langmuir, 1997, 13, 525– 532. 70 A. I. Greenwood, S. Tristram-Nagle and J. F. Nagle, Chem. Phys. Lipids, 2006, 143, 1–10. 71 M. D'Alessandro, A. Tenenbaum and A. Amadei, J. Phys. Chem. B, 2002, 106, 5050–5057. 72 D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, Academic press, 2001, vol. 1. 73 C. Hofs¨ ab, E. Lindahl and O. Edholm, Biophys. J., 2003, 84, 2192–2206. 74 P. Yee, J. K. Shah and E. J. Maginn, J. Phys. Chem. B, 2013, 117, 12556–12566. 75 E. H. Hill, K. Stratton, D. G. Whitten and D. G. Evans, Langmuir, 2012, 28, 14849–14854. 76 R. Groot and K. Rabone, Biophys. J., 2001, 81, 725–736. 77 M. Blesic, M. H. Marques, N. V. Plechkova, K. R. Seddon, L. P. N. Rebelo and A. Lopes, Green Chem., 2007, 9, 481–490. 78 J. Łuczak, C. Jungnickel, M. Markiewicz and J. Hupka, J. Phys. Chem. B, 2013, 117, 5653–5658. 79 C. Jungnickel, J. Łuczak, J. Ranke, J. F. Fern´ andez, A. M¨ uller and J. Th¨ oming, Colloids Surf., A, 2008, 316, 278–284. 80 T. Inoue, H. Ebina, B. Dong and L. Zheng, J. Colloid Interface Sci., 2007, 314, 236–241. 81 J. Łuczak, J. Hupka, J. Th¨ oming and C. Jungnickel, Colloids Surf., A, 2008, 329, 125–133. 82 C. Patrascu, F. Gauffre, F. Nallet, R. Bordes, J. Oberdisse, N. de Lauth-Viguerie and C. Mingotaud, ChemPhysChem, 2006, 7, 99–101. 83 B. Bhargava and M. L. Klein, Mol. Phys., 2009, 107, 393–401. 84 B. Bhargava and M. L. Klein, J. Phys. Chem. A, 2008, 113, 1898–1904. 85 J. N. Israelachvili, Intermolecular and Surface Forces, Academic Press, Revised Third Edition, 2011. 86 M. J. Rosen and J. T. Kunjappu, Surfactants and Interfacial Phenomena, John Wiley & Sons, 2012.
Soft Matter, 2014, 10, 8641–8651 | 8651
Soft Matter View Article Online
PAPER
Cite this: Soft Matter, 2014, 10, 8641
View Journal | View Issue
Amphiphilic interactions of ionic liquids with lipid biomembranes: a molecular simulation study† Brian Yoo,a Jindal K. Shah,ab Yingxi Zhua and Edward J. Maginn*a Current bottlenecks in the large-scale commercial use of many ionic liquids (ILs) include their high costs, low biodegradability, and often unknown toxicities. As a proactive effort to better understand the molecular mechanisms of ionic liquid toxicities, the work herein presents a comprehensive molecular simulation
study
on
the
interactions
of
1-n-alkyl-3-methylimidazolium-based
ILs
with
a
phosphatidylcholine (PC) lipid bilayer. We explore the effects of increasing alkyl chain length (n ¼ 4, 8, and 12) in the cation and anion hydrophobicity on the interactions with the lipid bilayer. Bulk atomistic molecular dynamics (MD) simulations performed at millimolar (mM) IL concentrations show spontaneous insertion of cations into the lipid bilayer regardless of the alkyl chain length and a favorable orientational preference once a cation is inserted. Cations also exhibit the ability to “flip” inside the lipid bilayer (as is common for amphiphiles) if partially inserted with an unfavorable orientation. Moreover, structural analysis of the lipid bilayer show that cationic insertion induces roughening of the bilayer surface, which may be a precursor to bilayer disruption. To overcome the limitation in the timescale of our simulations, free energies for a single IL cation and anion insertion have been determined based on potential of mean force calculations. These results show a decrease in free energy in response to both short and long alkyl Received 11th July 2014 Accepted 1st September 2014
chain IL cation insertion, and likewise for a single hydrophobic anion insertion, but an increase in free energy for the insertion of a hydrophilic chloride anion. Both bulk MD simulations and free energy
DOI: 10.1039/c4sm01528b www.rsc.org/softmatter
calculations suggest that toxicity mechanisms toward biological systems are likely caused by ILs behaving as ionic surfactants. [Yoo et al., Soft Matter, 2014].
1. Introduction Molecular simulations of synthetic molecules with model biomembranes have become a common method used to explore the physical mechanisms of molecular cytotoxicity. With recent developments in the design of new nanomaterials, there has been increasing concern for the safety of their widespread use, as it has been claimed that such synthetic materials can cause adverse effects toward biological systems.1 Many atomistic and coarse grained molecular simulation studies have provided extensive insight into the interactions that possibly govern the cytotoxic interactions of such materials. In particular, these studies have investigated the interactions between nanomaterials such as carbon nanotubes,2–4 fullerene,2,5–7 and functionalized metal-based nanoparticles8–12 with a hydrated lipid bilayer. Such studies have revealed spontaneous insertion or translocation through the lipid bilayer, provided insight into a
Department of Chemical and Biomolecular Engineering, University of Notre Dame, 182 Fitzpatrick Hall, Notre Dame, Indiana 46556, USA. E-mail: [email protected]; Fax: +1-574-631-8366; Tel: +1-574-631-5687
b
Center for Research Computing, University of Notre Dame, 111 Information Technology Center, Notre Dame, Indiana 46556, USA † Electronic supplementary information (ESI) available: Simulation input les are provided. See DOI: 10.1039/c4sm01528b
This journal is © The Royal Society of Chemistry 2014
the dominating energetic driving forces for molecular insertion or translocation while identifying the subsequent changes in the mechanical properties and structure (i.e. pore formation or rupture) of the lipid bilayer, and determined the likelihood for lipid peroxidation. Among these simulation studies, few have looked into the toxic interactions of a class of recently engineered ionic uids commonly known as ionic liquids (ILs).13–15 Within the past decade, the study of ILs has exploded as these liquids have been touted as promising designer solvents or additives for a wide variety of commercial applications.16–28 ILs are well known for their favorable characteristics: melting points lower than 100 C, good thermal stability, ionic conductivity, high chemical tunability, and negligible vapor pressure. ILs are oen called “green solvents” because their extremely low vapor pressures eliminate the possibility of their release into the atmosphere, unlike conventional volatile solvents. ILs still can nd their way into the environment, however, most likely via aqueous waste streams. Their low volatility makes them odorless, so that leaks are much harder to detect than with conventional solvents. Some ILs have been shown to be toxic toward several types of microorganisms,29,30 mammalian cell lines,31,32 and aquatic organisms33 at concentrations signicantly lower than those of conventional solvents. Thus “green solvent” is somewhat of a
Soft Matter, 2014, 10, 8641–8651 | 8641
View Article Online
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
Soft Matter
misnomer if it is applied to all ILs. At present, the toxicity mechanisms by which ILs cause such adversities toward biological systems are not well understood, in which case molecular simulations can be quite useful. Along with relieving environmental and safety concerns, better understanding of the molecular toxicity mechanisms of ionic liquids provides substantial insight in alternative energy research. Specically, in recent research related to enhanced biofuels production, ILs have shown signicant advantages in the pretreatment of lignocellulosic biomass in comparison to other conventional pretreatment solvent candidates because of their ability to overcome biomass recalcitrance at ambient conditions.17,18 Discoveries made recently by Ruegg et al.34 have shown that bacteria can be engineered to become more IL tolerant, a desirable attribute (for economic reasons) in the development of a “one-pot” process35,36 for fermentation-based biofuels production. For such a case, clear understanding of the molecular toxicity mechanisms of ILs with a cellular membrane could provide enhancements in the biological design choices for the bacterial membrane or in the molecular design of ILs. Recently, Bingham and Ballone13 used classical molecular dynamics (MD) simulations to explore the cytotoxic interactions between solutions of imidazolium-based ILs containing a 1-nbutyl-3-methylimidazolium ([C4mim]+) cation with a zwitterionic 1-palmiltoyl-2-oleoylphosphatidylcholine (POPC) lipid bilayer, a major component in the outer lipid bilayer leaet of eukaryotic cell membranes.37 Their atomistic simulations (over timescales on the order of 100 nanoseconds) have shown that the dominant interactions are likely due to the cations inserting into the lipid bilayer. Hydrophilic anions such as [Cl] and [PF6] tended to freely disperse in the aqueous environment, while in the case of ILs containing a very hydrophobic anion such as bis(triuoromethanesulfonyl)imide ([NTf2]), both the cation and the anion could insert into the lipid bilayer. Such attributes could provide an explanation for the increased toxicities observed for the hydrophobic imidazolium-based ionic liquids toward biological systems. The authors further observed that some of the commonly determined structural properties for the lipid bilayer systems such as the area and volume per lipid were not signicantly altered aer spontaneous insertion of the cation or IL pair. They also did not observe signicant morphology changes or disruption of the lipid bilayer over the length of their simulations. Another recent study by Kl¨ ahn and Zacharias15 has compared the interactions of 1-n-octyl-3-methylimidazolium chloride ([C8mim][Cl]) with model eukaryotic and cancerous cell membranes as means to understand the behavior of ionic liquids as antimicrobial agents. They determined the potential of mean force (PMF) or the free energy for insertion of a single [C8mim]+ into the biomembranes. The membrane models in their work were comprised of varying compositions of phosphatidylcholine, phosphatidylserine (negatively charged lipid), and cholesterol. Kl¨ ahn and Zacharias found that the free energy change for insertion of the [C8mim]+ cation into these membranes is negative and thus thermodynamically favorable. The presence of charged lipids and cholesterol in the bilayer could alter the free energy for the cation insertion. Their
8642 | Soft Matter, 2014, 10, 8641–8651
Paper
calculated free energies showed good agreement with those based on the octanol–water partition constant measured experimentally.38 As a comprehensive simulation study, in this work we have investigated the effect of cation alkyl chain length on the interactions of ILs of varied concentrations with a model lipid bilayer. We report the results of classical MD simulations of aqueous 1-n-alkyl-3-methylimidazolium-based ionic liquids of alkyl chain with n ¼ 4, 8 and 12 in contact with a POPC lipid bilayer. We also examine the effects associated with IL concentration and the hydrophobicity of the anion. The present work is motivated by recent experimental studies where several biological toxicity assays for imidazolium-based ILs have shown that increasing the length of the cation alkyl chain leads to signicantly higher toxicities, while the impact of varying the anion is more modest.30,31,39–41 In an independent study, we have conducted uorescence microscopy experiments on synthetic aPC bilayers in IL solutions and have observed that morphological changes in the bilayer depend on the concentration of the IL, the length of the alkyl chain and the type of anion paired with the cation. These results will be presented in a forthcoming paper. Our simulations are carried out at the nominal concentrations of the microscopy studies. We have also used umbrella sampling42 and the weighted histogram analysis method (WHAM)43 method to compute the PMF44 or free energy change (used interchangeably hereaer) for the insertion of the [C4mim]+ and 1-n-dodecyl-3-methylimidazolium ([C12mim]+) cations as well as the [Cl] and [NTf2] anions into the bilayer at low IL concentration.
2. 2.1
Simulation methods Classical MD simulation
Atomistic molecular dynamics simulations were performed for aqueous solutions of 1-n-alkyl-3-methylimidazolium ionic liquids with a POPC planar lipid bilayer patch using Gromacs 4.5.5.45,46 Cation alkyl chains of length n ¼ 4, 8 and 12 were simulated. Counter anions [Cl], [BF4], [PF6], and [NTf2] were selected based on their varying degree of hydrophobicity. Force eld parameters were taken from Liu et al.47 and Zhong et al.,48 and consist of united atom representations for each methyl group along the alkyl chains and imidazolium ring. The POPC force eld is based on the united atom Berger49 lipid model, which utilizes a combination of GROMOS87 and OPLS parameters with partial charges from ab initio quantum mechanical calculations.50 The model also utilizes a Ryckaert– Bellemans dihedral potential for the torsion angles along the sn1 and sn-2 lipid tails. The POPC lipid bilayer consisted of 128 lipids, with 64 lipids on the upper and lower leaet. The starting conguration for a hydrated lipid bilayer was provided elsewhere.51 Fig. 1(a) shows the different molecules and ions simulated, while Fig. 1(b) shows a snapshot of a representative IL/POPC/water system. For nonbonded interactions between the IL, water, and lipids, we used the Lorentz–Berthelot mixing rule. To account for the differences in 1–4 and charge interaction scaling between the ionic liquid (scaling factor ¼ 0.5) and the lipid
This journal is © The Royal Society of Chemistry 2014
View Article Online
Paper
Soft Matter
LX LY NPOPC =2
(1)
Vsys Nw Vw NPOPC
(2)
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
APOPC ¼ VPOPC ¼
Fig. 1 (a) Representation of the molecules and ions used in this study. Cations are in order of increasing alkyl side chain length: [C4mim]+, [C8mim]+, [C12mim]+ (top to bottom); anions are in order of increasing hydrophobicity: [Cl], [BF4], [PF6], and [NTf2] (left to right). The POPC lipid head (orange) and lipid tail (gray) are also shown. (b) The IL/ POPC/water bilayer system (water omitted for clarity). Box dimensions are roughly 6.3 6.3 17.5 nm in the xyz directions, respectively.
(scaling factor ¼ 1.0) models, we used the “half energy double pair-list” method for Gromacs, in which case the ionic liquid and lipid 1–4 nonbonded interactions were scaled by 0.5, while a pair list index for the lipid charge interactions were doubled. Periodic boundary conditions were used and electrostatic interactions were calculated using the particle mesh Ewald (PME) algorithm.52 The starting congurations of the IL/POPC/water systems were created by rst equilibrating the fully hydrated POPC bilayer conguration (128 POPC lipids solvated by 17 000 SPC water molecules) for 100 ns in the semi-isotropic NPT ensemble corresponding to near zero surface tension. A Nos´ e–Hoover53,54 thermostat was used to maintain the temperature at 298 K, which is above the gelation temperature (Tm ¼ 271 K) of POPC. A Parrinello–Rahman55 barostat was used to maintain the pressure at 1 bar. Ionic liquid ions were then placed roughly 2 nm above and below each leaet of the equilibrated POPC– water conguration aer removal of a few water molecules.56 Each system was subsequently run for another 100 ns. To compare the effect of the increasing alkyl side chain length, we ran simulations with [C4mim][Cl], [C8mim][Cl], and [C12mim][Cl]. To examine the effect of anion hydrophobicity, we ran simulations with [Cl], [BF4], [PF6], and [NTf2] paired with the [C4mim]+ cation. Finally, different amounts of cations and anions were used to study the effect of IL concentration. In all cases, concentration were below the experimental solubility limit of [C4mim][NTf2], which among the studied ILs is the least soluble in water.57
2.2
Lipid bilayer analysis
Equilibrium properties for the lipid bilayer were determined based on block averaging the last 20 ns of each independently run MD simulation.58 The area per lipid (APOPC) and volume per lipid (VPOPC) for the hydrated lipid bilayer were calculated via the following expressions
This journal is © The Royal Society of Chemistry 2014
where LX and LY and the box lengths in the X and Y dimensions respectively; NPOPC and Nw are the number of lipids and waters within the system, respectively, and Vw and Vsys are the volume occupied by a single water molecule and the total volume of the system. An independent simulation of pure water was performed to determine Vw using the same simulation conditions as the IL/POPC/water systems. Additionally, the bilayer thickness (h) was determined for each IL/POPC/water system using the GridMAT-MD analysis tool.59 The phosphorous atom in the POPC head group was selected as a reference atom, upon which a grid (150 149 points) was used to determine their locations and to calculate the local thicknesses. Total bilayer thicknesses were taken by averaging data from the last 20 ns every 200 ps (100 frames in total). Subsequently, they were used to set the integration bounds in the density proles when calculating the number of inserted ILs (see Table 1). Lateral self diffusion coefficients (DL) were determined by the Einstein relation: E 1 D 2 2 DL ¼ lim (3) jxðt0 Þ xðt0 þ tÞj þ jyðt0 Þ yðt0 þ tÞj t/N 2dt t0 where d is the dimensionality (d ¼ 2). Artifacts from the lateral center of mass motion of the bilayer leaets were subtracted from the lipid displacements. To determine DL, each of the 100 ns simulation trajectories was split into ve separate 20 ns trajectories. Slopes of the mean square displacement were taken from 2 to 8 ns and averaged from each of the 20 ns trajectories. This method has been shown to provide reasonable estimates of diffusion coefficients for a bilayer when compared to experimental uorescence correlation spectrometry measurements.60–62 Deuterium order parameters, used to determine the liquid ordered Lo or liquid disordered Ld phase of the lipid bilayer, were determined for individual carbon atoms along the sn-1 and sn-2 tails using the second order Legendre polynomial SCD ¼
1 2 3 cos qCD 1 2
(4)
where qCD is the angle between the carbon–deuterium CD bond (hydrogen for the simulations) and the membrane normal. A value of SCD < 0.35 typically indicates the bilayer is in the Ld uid phase.63 To determine the positions of the hydrogens along the saturated united atom carbons of the lipid tails, a tetrahedral geometry between the central and adjacent carbon atoms was assumed. For unsaturated carbons along the sn-2 tail of the POPC lipid, the vector normal to the plane of the unsaturated carbon and adjacent carbons was used to determine the hydrogen position. It should be noted here that the utility in Gromacs 4.5.5 used to determine the order parameter in this work had been patched to x a bug which assumed an incorrect
Soft Matter, 2014, 10, 8641–8651 | 8643
View Article Online
Soft Matter
Paper
Structural and dynamic properties determined for each IL/POPC/water system over the course of the last 20 ns including: the number of IL ion pairs (nIL), number of water molecules nH2O, concentration of IL (cIL), area (APOPC) per lipid, bilayer thickness (h), lipid lateral diffusion coefficient (DL), average order parameters (hSCDi) for the sn-1 and sn-2 tails, and the number of cations inserted (nIL+). (A nIL+ as a non-integer indicates that the IL is partially inserted.)
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
Table 1
System
nIL
nH2O
cIL/mM
APOPC/nm2
h/nm
POPC–water [C4mim][Cl]
— 5 15 30 15 2 4 5 5 5
17 152 16 877 16 877 16 811 16 877 17 014 16 611 16 877 16 877 16 877
— 16 50 100 50 7 13 16 16 16
0.615 0.612 0.609 0.614 0.604 0.604 0.616 0.586 0.598 0.612
3.885 3.907 3.975 4.005 3.988 3.944 3.957 4.076 3.967 3.939
[C8mim][Cl] [C12mim][Cl] [C4mim][BF4] [C4mim][PF6] [C4mim][NTf2]
0.008 0.006 0.004 0.003 0.002 0.002 0.007 0.001 0.006 0.006
0.172 0.198 0.289 0.277 0.195 0.262 0.265 0.196 0.184 0.262
geometry when determining hydrogen positions for the unsaturated carbons. Electrostatic potential proles were based on integrations of the Poisson equation d2 jðzÞ rc ðzÞ ¼ dz2 30
(5)
where j(z) is the electrostatic potential along the z-direction of the lipid bilayer, rc(z) is the charge density and 30 is the vacuum permittivity. In order to remove artifacts arising from the bilayer center of mass motion, a position correctional term was subtracted from the Z distances.64
2.3 Potential of mean force calculation using umbrella sampling To determine the PMFs for inserting a single [C4mim]+, [C12mim]+, [Cl] and [NTf2] into the POPC bilayer, we used the umbrella sampling and WHAM analysis method implemented in Gromacs.65 For each cation or anion of interest, we allowed for its counterion to freely disperse in the aqueous environment, while implementing a restrictive harmonic potential (in the Z direction of the bilayer plane) on one of the terminal atoms in the cation or anion. To be consistent with our bulk simulations, we paired a [Cl] anion with the cations [C4mim]+ or [C12mim]+ to determine the free energy change for cationic insertion. Likewise, we paired a [C4mim]+ with the [Cl] and [NTf2] anions when determining the anion PMFs. Simulation parameters and the number of water molecules in the systems were the same as those of our bulk simulations. A simulation time of 80 ns for each of the 40–50 frames, with 0.1 nm spacing, was used to obtain better convergence and overlap of energies and to allow for substantial structural relaxation time for the lipid bilayer. More details on the simulation parameters for both bulk MD simulations and umbrella sampling-based PMF calculations can be found in the Supplemental Information. All necessary input les to reproduce our simulations are also provided. Visualizations for our systems have been created with VMD.66
8644 | Soft Matter, 2014, 10, 8641–8651
3. 3.1
DL/108 cm2 s1
hSCDisn-1
hSCDisn-2
nIL+ inserted
1.20 0.41 3.62 0.47 3.88 0.71 3.70 0.45 4.63 1.32 1.17 0.19 1.17 0.17 3.90 0.79 4.72 0.66 3.89 0.35
0.184 0.034 0.185 0.032 0.194 0.039 0.197 0.038 0.201 0.042 0.193 0.034 0.186 0.032 0.211 0.038 0.197 0.030 0.196 0.033
0.146 0.049 0.148 0.052 0.153 0.049 0.156 0.051 0.155 0.054 0.154 0.048 0.149 0.050 0.169 0.053 0.151 0.045 0.153 0.050
— 2.0 6.6 9.5 5.0 2.0 2.7 2.0 1.0 4.0
Results and discussion Bulk MD simulation
3.1.1 Fully hydrated POPC bilayer validation. To draw comparison between the lipid bilayer structural properties for the IL/POPC/water MD simulations, a fully hydrated POPC lipid bilayer system was simulated and equilibrium properties were compared with those from the literature.13,60,67–69 The average area per lipid, bilayer thickness, lateral diffusion coefficient, and average order parameters are listed in Table 1. We obtain an average area per lipid of 0.615 nm2 with bilayer thicknesses of 3.885 nm, which is in close agreement with the simulation results of Mukhopadhyay et al.67 The average area per lipid is slightly lower than the experimental value of 0.64 0.01 nm2 at 298 K.69 The volume per lipid and volume per water are 1.205 nm3 and 0.0306 nm3, respectively, are also in close agreement with previous simulation results and are in reasonable agreement with experimental studies (1.256 nm3 at 303 K and 0.030 nm3 at 298 K, respectively).70 Lateral diffusion coefficients of the lipids obtained in this work are roughly a factor of 3 lower than those calculated by B¨ ockmann et al.,60 who obtained a lateral diffusion coefficient of 3.9 108 cm2 s1 at 300 K using the same lipid and water force elds. This disparity is likely due to the differences in our simulation conditions. The van der Waals and electrostatic cutoffs in our study (rcut ¼ 1.2 nm) are slightly longer than those used by B¨ ockmann et al.60 (rcut ¼ 1.0 nm). A longer cutoff was used in the present work in order to be consistent with the cutoff used in parameterizing the IL force eld. Another difference between the present simulations and those of B¨ ockmann et al.60 is that we used a Nos´ e–Hoover thermostat, while B¨ ockmann et al. used a Berendsen thermostat. It is known that a Berendsen thermostat incorrectly samples canonical ensemble uctuations, which in turn can lead to faster dynamics.71,72 3.1.2 IL–lipid bilayer interaction and the effect of IL cation alkyl chain length. Independent simulations at concentrations of roughly 50 mM [C4mim][Cl], 50 mM [C8mim][Cl], and two systems at 7 and 13 mM [C12mim][Cl], were performed to observe the effect of the increased cation alkyl side chain length. For the [C12mim][Cl] systems, we did not run a full 100
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
Paper
ns simulation at a 50 mM concentration (corresponding to 15 IL pairs) because the cations tended to aggregate into micelles. For this system our simulation timescales are too short (even aer 60 ns) to observe these aggregates insert into the lipid bilayer. For the other systems, we observed that a fraction of the cations spontaneously inserted into the lipid bilayer, while all of the [Cl] anions remain freely dispersed in water. Cations that insert into the bilayer during the simulation tend to favor orientations in which the alkyl side chain (“IL tail”) inserts deeply into the hydrophobic region of the lipid bilayer, while the imidazolium ring (“IL head”) interacts with the lipid head groups. (The terms “IL tail” and “IL head” are used to draw the similarities of these ILs with surfactants.) For some of the IL cations, we observe partial insertion of the IL head, and eventual reorientation or “ipping”, such that the IL tail is inserted into the hydrophobic regime of the bilayer, as shown in Fig. 2. Also indicated in the gure are the angles qNC between the membrane leaet normal and the vector between a nitrogen from the imidazolium ring and the terminal methyl group of the alkyl side chain. When qNC is approximately greater than 120 , the tail can be assumed to be completely inserted into the bilayer. For the 50 mM [C4mim][Cl], 50 mM [C8mim][Cl] and 13 mM [C12mim][Cl] systems, we observe at least one partially inserted IL ip, such that the IL tail is inserted deeper
Soft Matter
into the lipid bilayer, while the IL head remains near the lipid and water interface. However, the time it takes for a partially inserted IL cation to ip is signicantly longer than the time it takes for a cation to insert directly in a “tail rst” orientation. Normalized density proles of the orientation-favored insertion of the cations for the 50 mM [C4mim][Cl], 50 mM [C8mim][Cl], and 7 mM [C12mim][Cl] systems are plotted in Fig. 3(a)–(c). The density proles have been averaged based on the last 20 ns of the simulation run. For the longer alkyl chain cations, the density proles show deeper penetration of the alkyl tails into the lipid bilayer with density maxima closer to the lipid bilayer center, while the imidazolium rings stay at the water–lipid interface. The [Cl] anion density proles show no indication of penetration into the bilayer or strong reassociation with the inserted cations. The results of the SCD order parameter calculations are shown in Fig. 4. Cation insertion tends to increase the local ordering of the sn-1 and sn-2 lipid tails. The carbons in the lipid tails near the lipid head groups (cn with n ¼ 2 to 5) exhibit slightly increased ordering compared with those at the end of the tails, as indicated by the increase in SCD for the [C4mim][Cl] system. For the 50 mM [C8mim][Cl] system, insertion of the cations causes a greater increase in SCD than that of [C4mim][Cl], particularly for the sn-1 carbons cn with n ranging from 2 to 12. For the 7 and 13
Fig. 2 Snapshots of partially inserted IL cation “flipping” into the bilayer during the course of the simulation run for systems: (a) 50 mM [C4mim][Cl], (b) 50 mM [C8mim][Cl], and (c) 13 mM [C12mim][Cl]. The angles qNC between the membrane leaflet normal and vector between a nitrogen from the imidazolium ring and the terminal methyl group of the alkyl side chain (NC) are indicated. Arrows pointing in the direction of the NC-vector are shown for guidance. The rest of the anions, cations, and water molecules in each system are omitted for clarity.
This journal is © The Royal Society of Chemistry 2014
Soft Matter, 2014, 10, 8641–8651 | 8645
View Article Online
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
Soft Matter
Paper
Fig. 3 Normalized density r* profiles along the Z axis for lipid bilayer systems in aqueous solutions of (a) 50 mM [C4mim][Cl], (b) 50 mM [C8mim][Cl], (c) 7 mM [C12mim][Cl], (d) 16 mM [C4mim][BF4], (e) 16 mM [C4mim][PF6], and (f) 16 mM [C4mim][NTf2]. Density profiles for the POPC bilayer (black) has been multiplied by 2 for clarity. For almost all cases, cations and anions are freely disassociated. Density profiles for the IL cation head (blue) and tail (cyan) are plotted separately. Those for anions are also plotted (green).
Fig. 4 (a) Deuterium order parameter for the sn-1 (solid) and sn-2 (dashed) tails for the POPC–water system (black), the 50 mM [C4mim] [Cl] (blue) and the 50 mM [C8mim][Cl] (red) system.
mM [C12mim][Cl] systems, we observe a similar trend in the sn-1 and sn-2 order parameters, although a shi in the order parameters is not as noticeable for these systems, due to the lower number of inserted cations. For the [C4mim][Cl] and [C8mim][Cl] systems, lateral diffusion coefficients of the lipids tend to increase by roughly a factor of 2–4 as compared to the lipid without ILs. An increase in the lateral diffusion coefficient and increase in the sn-1 and sn-2 order parameters may suggest that insertion of the IL cations induces a tendency for greater collective lateral motion, in which case the lipid head groups in these systems are possibly conned locally due to the presence of a nearby inserted cation and forced to displace toward a nearby lipid.
8646 | Soft Matter, 2014, 10, 8641–8651
For the [C12mim][Cl] systems, we do not observe a similar increase in the lipid lateral diffusion coefficient. This can be attributed to the deeper insertion of the longer alkyl chain into the bilayer, which may interact much like a lipid tail group. Recently, it has been shown that a large electrostatic potential gradient induced by the presence of ions and cationic nanoparticles can drive pore formation within a double bilayer.11 It should be noted that we observe slight changes to the electrostatic potential aer cation insertion for these systems, although it is difficult to determine if such deviations were directly caused by cation insertion. Previous studies have shown that slight deviations in the electrostatic potential can be caused by uctuations and averaging over an inadequate simulation time (4.2 nm) occur more frequently with an increased number of inserted cations. The positions of the inserted cations (represented by the magenta dots in Fig. 6) do not appear to strongly correlate with the locations of the regions of high or low bilayer thicknesses. Hill et al.75 have shown that bilayer thickness deformation, particularly when the thicknesses decreases signicantly (< 0.8 nm) in comparison to an equilibrated hydrated lipid bilayer, tends to eventually lead to formation of a porous channel in the low thickness regions, allowing for the transport of water or ions across the bilayer. However, the reason for the pore formation in their simulations was due to strong charge interactions between aggregates of the inserted cytotoxic molecules with the charged head groups of the nearby lipids in the bilayer hydrophobic region. We do not observe any pore formation or diffusion of inserted cations into large aggregates for the IL/ water/POPC systems within the simulated timescale, though the locations of the inserted cations in the 100 mM system are positioned relatively close to one another. It is thus suggested that the ILs denitely cause increased surface roughness on 100 ns timescales, and this could be a precursor to bilayer disruption. As was the case with the increasing chain length, the diffusion coefficients are roughly 3 times larger for all the systems when compared to the hydrated lipid bilayer. Interestingly, although the lateral diffusion coefficient is higher for the cases when ILs insert into the bilayer, it does not appear to change as the number of inserted ions increases. Given that a small number of systems were studied, it is not possible to know if this is a general phenomena. Likewise, order parameters increase very slightly with increasing IL concentration, while the average overall bilayer thickness, based on the locations of the phosphorus atoms in the lipid head group, is insensitive to IL concentration.
3.2
Free energy change upon cation or anion insertion
The Gibbs free energy change DG for a single cation or anion insertion from bulk aqueous solution was determined from potential of mean force calculations in the semi-isotropic NPT ensemble. From our observation of the orientational preference of inserted IL cations in our bulk simulations, we decided to select as the reaction coordinate for [C12mim]+ and [C4mim]+
Soft Matter, 2014, 10, 8641–8651 | 8647
View Article Online
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
Soft Matter
Paper
Fig. 6 Snapshots of configurations (above) and bilayer thickness contour plots (h) (below). The contour plots were averaged over the last 20 ns of the simulation. (a) The hydrated lipid bilayer, (b) 16 mM [C4mim][Cl], (c) 50 mM [C4mim][Cl], and (d) 100 mM [C4mim][Cl]. Locations of inserted cation atoms from each frame are indicated by the dots in magenta.
the distance along the Z direction between the center of mass of the lipid bilayer and the terminal methyl group on the alkyl chain, as shown in Fig. 7(a). We found that the selection of this reaction coordinate allows for quicker convergence in the free energy prole due to the relatively long ipping time of the aliphatic [C12mim]+. For [NTf2], one of the united atom CF3 groups was chosen as the reference. The center of mass was chosen as the reference for the other anions. Fig. 7(b) and (c) plots the free energy proles for insertion of a single [C4mim]+, [C12mim]+, [Cl] and [NTf2] from bulk solution. The free energies of both [C4mim]+ and [C12mim]+ (paired with a [Cl]) monotonically decrease from bulk solution as the cations approach the lipid bilayer. The insertion process is essentially barrierless, indicating that the insertion process is thermodynamically spontaneous. The minimum in the free energy for the [C4mim]+ insertion is roughly 27 kJ mol1 at z z 1.2 nm, while for [C12mim]+ the minimum is 37 kJ mol1 at z z 0.07 nm. The location of each free energy minimum along the reaction coordinate corresponds
to the cation head interacting strongly with the lipid head groups and the depth of the cation tail insertion. The sharp increase in the free energy prole starting at about z < 0.5 nm for [C4mim]+ and the slight increase at z < 0.07 nm for [C12mim]+ indicate that the cation head groups prefer to reside at the bilayer interface with the tail group extended into the hydrophobic regime of the lipid bilayer. Thus the cations do not cross the lipid bilayer, but rather remain inserted in the bilayer. Kl¨ ahn and Zacharias have reported a minimum free energy for [C8mim]+ insertion into a POPC bilayer of 15 kJ mol1 at about z z 0.6 nm using a similar reaction coordinate. Their value is roughly a factor of 2 smaller in magnitude than those of our [C4mim]+ and [C12mim]+ estimates. Based on the qualitative trends of our free energy minimum, we would expect our values to be relatively close to those of Kl¨ ahn and Zacharias. The minimum free energy for [C12mim]+ should be less, while that for the [C4mim]+ should be greater. The most likely explanation for this discrepancy is in the difference in sampling time. While Kl¨ ahn and Zacharias obtained a minimum free energy close in
Fig. 7 Free energy profiles DG along the reaction coordinate z. (a) The reaction coordinates z for cation or anion are defined based on the Z distance between the reference atom and the center of mass of the lipid bilayer (z ¼ 0). Free energy profiles for the insertion of a single (b) [C4mim]+ (black) and [C12mim]+ (blue), and (c) [Cl] (green) and [NTf2] (red) anion into a POPC lipid bilayer. Uncertainties are indicated by the shaded regions.
8648 | Soft Matter, 2014, 10, 8641–8651
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
Paper
value to experimental studies, they used a sampling time of 10 ns per conguration. We nd that although we obtain a free energy minimum close in value with theirs within 10–20 ns, the free energies do not show good convergence as they continue to decrease until about 50–80 ns. A longer sampling time per conguration allows for longer structural relaxation or equilibration of the bilayer in response to the inserted IL. Another reason for the discrepancy can arise from the differences in the lipid and IL force elds, as their IL force eld partial charges were reparameterized to account for the solvation effects of water. Nevertheless, the shape in their free energy curve and the location of the minimum free energy obtained are consistent and as expected with those of our work. For the most hydrophobic [NTf2] anion the free energy minimum is approximately 10 kJ mol1 at z z 0.6 nm, while for the hydrophilic [Cl] anion there is only a shallow minimum free energy at the bilayer surface. As the [Cl] anion moves into the bilayer, the free energy is strongly repulsive. The PMF results conrm the results from the bulk MD simulations. There is a larger thermodynamic driving force for the longer alkyl chain cations to insert into the bilayer, though even the short chain cations have negative free energies for insertion. The hydrophobic [NTf2] anion has a moderate tendency to insert into the bilayer, while the hydrophilic [Cl] shows no tendency to independently insert into the bilayer. In addition, we nd that the free energy change in the insertion of [NTf2] is independent of the [C4mim]+ cation insertion. Note that the PMF calculations were all performed on a single ion. The bulk MD simulations suggest that the presence of an inserted cation could lower the free energy barrier for subsequent anion insertion, but this was not examined in the present work.
4. Conclusions Molecular dynamics simulations and potential of mean force calculations were used to study the interactions of a range of different imidazolium-based ionic liquids with a POPC lipid bilayer. The simulations show that the ionic liquid cations tend to spontaneously insert into the lipid bilayer, while exhibiting strong amphiphilic interactions. Specically, we nd that as the cations insert, regardless of the alkyl chain length, cations favor orientation such that the imidazolium ring and the alkyl chain are strongly associated with the lipid head and tail groups, respectively. IL cations that have partially inserted initially with the imidazolium ring in the bilayer demonstrate the ability to ip such that the alkyl chain is completely inserted into the hydrophobic region of the bilayer. This ipping time tends to be much longer for the more aliphatic imidazolium cations of increasing alkyl chain length ([C8mim][Cl] and [C12mim][Cl]). Insertion of the cations tends to change the structural and dynamic properties of the bilayer. In particular, the ionic liquids tend to roughen the surface of the bilayer, with higher concentrations of ionic liquids leading to greater surface roughness. On the timescales of the simulation (100 ns) we do not observe direct disruption or pore formation in the bilayer, but this roughening could be a precursor for bilayer disruption, a potential cause of ionic liquid toxicity.
This journal is © The Royal Society of Chemistry 2014
Soft Matter
The effect of the anion hydrophobicity of the ionic liquid was studied by pairing the same cation ([C4mim]+) with [Cl], [BF4], [PF6], and [NTf2]. The latter three anions are more hydrophobic than [Cl] and tend to have a slightly greater effect on the structural properties of the bilayer than that of [Cl]. Within the timescales of the simulations, the [NTf2] anion is the only hydrophobic anion that inserts into the bilayer and reassociates with an inserted cation. Unlike the cations, the degree of hydrophobicity in the anion does not correlate well with the changes in the structural and dynamic properties of the bilayer. Umbrella sampling potential of mean force calculations were performed to compute the free energies for inserting various ions into the bilayer. The results conrm the ndings of the bulk MD simulations. Insertion of the cations is thermodynamically favorable, where the free energy change for insertion of [C4mim]+ is 27 kJ mol1 and is 37 kJ mol1 for the insertion of [C12mim]+. The [C12mim]+ cation also inserts more deeply into the bilayer than does [C4mim]+. The free energy change for insertion of the hydrophobic [NTf2] anion is also thermodynamically favorable at 10 kJ mol1. The [Cl] anion exhibits a positive free energy as it approaches the lipid bilayer, suggesting that it is more thermodynamically favorable for it to remain hydrated in the aqueous phase or to reside near the charged head groups of the bilayer. This study demonstrates that the deep insertion of imidazolium cations with long alkyl chains into the lipid bilayer is a likely mechanism for greater cytotoxicity toward biological systems. We also conrm insertion of hydrophobic anions with other computational work.13 Although we do not observe direct pore or channel formation, previous simulation studies have shown that the presence of charged molecule aggregates (particularly those that are cationic) or deeply inserted nonionic surfactants into a lipid bilayer can cause such structural damage.11,75,76 The results of this study suggest that the underlying toxicity mechanisms of the ionic liquids are governed by the ionic liquids essentially behaving as ionic surfactants. This classication of ionic liquids as ionic surfactants is by no means a novel discovery, as many studies have shown that ILs with longer alkyl chains (typically for [Cnmim]+ with n > 7) exhibit critical micelle concentrations,77–84 a common characteristic of surfactants.85 Within the context of understanding IL toxicity mechanisms toward biological systems, however, this study may suggest that assessments on the environmental safety concerns for imidazolium-based ionic liquids (including those with shorter alkyl chains such as [C4mim]+) can be made similarly to that of more commonly used ionic surfactants in industry, although the biodegradability of these ILs should still be considered.86 Further studies such as those modeling larger systems or longer simulation times, with either atomistic simulations or coarse grained simulations, need to be performed to validate such a mechanistic basis.
Acknowledgements Funding for this project was provided by the National Science Foundation (CBET-1134238). Computational resources were
Soft Matter, 2014, 10, 8641–8651 | 8649
View Article Online
Soft Matter
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
provided by Notre Dame's Center for Research Computing. We are also grateful for Prof. Peter Tieleman for providing the lipid force eld and bilayer initial conguration les.
References 1 A. Nel, T. Xia, L. M¨ adler and N. Li, Science, 2006, 311, 622– 627. 2 L. Monticelli, E. Salonen, P. C. Ke and I. Vattulainen, So Matter, 2009, 5, 4433–4445. 3 J. Liu and A. J. Hopnger, Chem. Res. Toxicol., 2008, 21, 459– 466. 4 S. Pogodin and V. A. Baulin, ACS Nano, 2010, 4, 5293–5300. 5 R. Qiao, A. P. Roberts, A. S. Mount, S. J. Klaine and P. C. Ke, Nano Lett., 2007, 7, 614–619. 6 J. Wong-Ekkabut, S. Baoukina, W. Triampo, I.-M. Tang, D. P. Tieleman and L. Monticelli, Nat. Nanotechnol., 2008, 3, 363–368. 7 R. DeVane, A. Jusu, W. Shinoda, C.-c. Chiu, S. O. Nielsen, P. B. Moore and M. L. Klein, J. Phys. Chem. B, 2010, 114, 16364–16372. 8 K. Yang and Y.-Q. Ma, Nat. Nanotechnol., 2010, 5, 579–583. 9 V. V. Ginzburg and S. Balijepalli, Nano Lett., 2007, 7, 3716– 3722. 10 R. C. Van Lehn and A. Alexander-Katz, So Matter, 2011, 7, 11392–11404. 11 J. Lin and A. Alexander-Katz, ACS Nano, 2013, 7, 10799– 10808. 12 R. C. Van Lehn and A. Alexander-Katz, So Matter, 2014, 10, 648–658. 13 R. J. Bingham and P. Ballone, J. Phys. Chem. B, 2012, 116, 11205–11216. 14 S. Cromie, M. Del P´ opolo and P. Ballone, J. Phys. Chem. B, 2009, 113, 11642–11648. 15 M. Kl¨ ahn and M. Zacharias, Phys. Chem. Chem. Phys., 2013, 15, 14427–14441. 16 H. Wu and E. J. Maginn, Fluid Phase Equilib., 2014, 368, 72– 79. 17 S. Singh, B. A. Simmons and K. P. Vogel, Biotechnol. Bioeng., 2009, 104, 68–75. 18 M. Francisco, A. N. Mlinar, B. Yoo, A. T. Bell and J. M. Prausnitz, Chem. Eng. J., 2011, 172, 184–190. 19 H. Rodr´ıguez, S. Padmanabhan, G. Poon and J. M. Prausnitz, Bioresour. Technol., 2011, 102, 7946–7952. 20 L. D. McIntosh, T. Kubo and T. P. Lodge, Macromolecules, 2014, 47, 1090–1098. 21 A. S. Pensado, F. Malberg, M. F. C. Gomes, A. A. H. Padua, J. Fernandez and B. Kirchner, RSC Adv., 2014, 4, 18017– 18024. 22 W. Afzal, B. Yoo and J. M. Prausnitz, Ind. Eng. Chem. Res., 2012, 51, 4433–4439. 23 K. A. Maerzke, G. S. Goff, W. H. Runde, W. F. Schneider and E. J. Maginn, J. Phys. Chem. B, 2013, 117, 10852–10868. 24 X. Kang, J. Zhang, W. Shang, T. Wu, P. Zhang, B. Han, Z. Wu, G. Mo and X. Xing, J. Am. Chem. Soc., 2014, 136, 3768–3771. 25 J. Stoimenovski, D. R. MacFarlane, K. Bica and R. D. Rogers, Pharm. Res., 2010, 27, 521–526.
8650 | Soft Matter, 2014, 10, 8641–8651
Paper
26 R. Ferraz, L. C. Branco, C. Prudencio, J. P. Noronha and ˇ. Petrovski, ChemMedChem, 2011, 6, 975–985. Z 27 S. Li, H. Gao and J. M. Shreeve, Angew. Chem., Int. Ed., 2014, 53, 2969–2972. 28 B. Yoo, W. Afzal and J. M. Prausnitz, Ind. Eng. Chem. Res., 2012, 51, 9913–9917. 29 R. P. Swatloski, J. D. Holbrey, S. B. Memon, G. A. Caldwell, K. A. Caldwell and R. D. Rogers, Chem. Commun., 2004, 668–669. 30 D. J. Couling, R. J. Bernot, K. M. Docherty, J. K. Dixon and E. J. Maginn, Green Chem., 2006, 8, 82–90. 31 R. A. Kumar, N. Papa¨ıconomou, J.-M. Lee, J. Salminen, D. S. Clark and J. M. Prausnitz, Environ. Toxicol., 2009, 24, 388–395. 32 J. Ranke, A. M¨ uller, U. Bottin-Weber, F. Stock, S. Stolte, J. Arning, R. St¨ ormann and B. Jastorff, Ecotoxicol. Environ. Saf., 2007, 67, 430–438. 33 A. Latała, P. Stepnowski, M. Nedzi and W. Mrozik, Aquat. Toxicol., 2005, 73, 91–98. 34 T. L. Ruegg, E.-M. Kim, B. A. Simmons, J. D. Keasling, S. W. Singer, T. S. Lee and M. P. Thelen, Nat. Commun., 2014, 5, DOI: 10.1038/ncomms4490. 35 H. Zhao, S. Xia and P. Ma, J. Chem. Technol. Biotechnol., 2005, 80, 1089–1096. 36 P. W. Wolski, D. S. Clark and H. W. Blanch, Green Chem., 2011, 13, 3107–3110. 37 G. Van Meer, D. R. Voelker and G. W. Feigenson, Nat. Rev. Mol. Cell Biol., 2008, 9, 112–124. 38 M. Matzke, S. Stolte, K. Thiele, T. Juffernholz, J. Arning, J. Ranke, U. Welz-Biermann and B. Jastorff, Green Chem., 2007, 9, 1198–1207. 39 P. Stepnowski, A. Skladanowski, A. Ludwiczak and E. Laczy´ nska, Hum. Exp. Toxicol., 2004, 23, 513–517. 40 B. Jastorff, R. St¨ ormann, J. Ranke, K. M¨ olter, F. Stock, B. Oberheitmann, W. Hoffmann, J. Hoffmann, M. N¨ uchter, B. Ondruschka, et al., Green Chem., 2003, 5, 136–142. 41 F. Stock, J. Hoffmann, J. Ranke, R. St¨ ormann, B. Ondruschka and B. Jastorff, Green Chem., 2004, 6, 286–290. 42 B. Roux, Comput. Phys. Commun., 1995, 91, 275–282. 43 S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen and P. A. Kollman, J. Comput. Chem., 1992, 13, 1011–1021. 44 D. Chandler, Introduction to Modern Statistical Mechanics, by David Chandler, pp. 288. Foreword by David Chandler, Sep 1987, ISBN-10: 0195042778. ISBN-13: 9780195042771, Oxford University Press, 1987, vol. 1, pp. 201–202. 45 E. Lindahl, B. Hess and D. van der Spoel, Molecular Modeling Annual, 2001, 7, 306–317. 46 B. Hess, C. Kutzner, D. Van Der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447. 47 Z. Liu, X. Wu and W. Wang, Phys. Chem. Chem. Phys., 2006, 8, 1096–1104. 48 X. Zhong, Z. Liu and D. Cao, J. Phys. Chem. B, 2011, 115, 10027–10040. 49 O. Berger, O. Edholm and F. J¨ ahnig, Biophys. J., 1997, 72, 2002–2013. 50 S.-W. Chiu, M. Clark, V. Balaji, S. Subramaniam, H. L. Scott and E. Jakobsson, Biophys. J., 1995, 69, 1230–1245.
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 02 September 2014. Downloaded by State University of New York at Stony Brook on 24/10/2014 23:28:38.
Paper
51 http://wcm.ucalgary.ca/tieleman/downloads, accessed 07-012013. 52 T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092. 53 S. Nos´ e, J. Chem. Phys., 1984, 81, 511–519. 54 W. G. Hoover, Phys. Rev. A, 1985, 31, 1695–1697. 55 M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182– 7190. 56 L. Mart´ınez, R. Andrade, E. G. Birgin and J. M. Mart´ınez, J. Comput. Chem., 2009, 30, 2157–2164. 57 M. G. Freire, P. J. Carvalho, R. L. Gardas, I. M. Marrucho, L. M. Santos and J. A. Coutinho, J. Phys. Chem. B, 2008, 112, 1604–1610. 58 M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford university press, 1989. 59 W. J. Allen, J. A. Lemkul and D. R. Bevan, J. Comput. Chem., 2009, 30, 1952–1958. 60 R. A. B¨ ockmann, A. Hac, T. Heimburg and H. Grubm¨ uller, Biophys. J., 2003, 85, 1647–1655. 61 S. W. Siu, R. V´ acha, P. Jungwirth and R. A. B¨ ockmann, J. Chem. Phys., 2008, 128, 125103. 62 A. A. Gurtovenko and I. Vattulainen, J. Phys. Chem. B, 2008, 112, 1953–1962. 63 O. G. Mouritsen and K. Jrgensen, Chem. Phys. Lipids, 1994, 73, 3–25. 64 A. A. Gurtovenko and I. Vattulainen, J. Chem. Phys., 2009, 130, 215107. 65 J. S. Hub, B. L. de Groot and D. van der Spoel, J. Chem. Theory Comput., 2010, 6, 3713–3720. 66 W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38. 67 P. Mukhopadhyay, H. J. Vogel and D. P. Tieleman, Biophys. J., 2004, 86, 337–345. 68 J. Seelig and N. Waespe-Sarcevic, Biochemistry, 1978, 17, 3310–3315.
This journal is © The Royal Society of Chemistry 2014
Soft Matter
69 B. K¨ onig, U. Dietrich and G. Klose, Langmuir, 1997, 13, 525– 532. 70 A. I. Greenwood, S. Tristram-Nagle and J. F. Nagle, Chem. Phys. Lipids, 2006, 143, 1–10. 71 M. D'Alessandro, A. Tenenbaum and A. Amadei, J. Phys. Chem. B, 2002, 106, 5050–5057. 72 D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, Academic press, 2001, vol. 1. 73 C. Hofs¨ ab, E. Lindahl and O. Edholm, Biophys. J., 2003, 84, 2192–2206. 74 P. Yee, J. K. Shah and E. J. Maginn, J. Phys. Chem. B, 2013, 117, 12556–12566. 75 E. H. Hill, K. Stratton, D. G. Whitten and D. G. Evans, Langmuir, 2012, 28, 14849–14854. 76 R. Groot and K. Rabone, Biophys. J., 2001, 81, 725–736. 77 M. Blesic, M. H. Marques, N. V. Plechkova, K. R. Seddon, L. P. N. Rebelo and A. Lopes, Green Chem., 2007, 9, 481–490. 78 J. Łuczak, C. Jungnickel, M. Markiewicz and J. Hupka, J. Phys. Chem. B, 2013, 117, 5653–5658. 79 C. Jungnickel, J. Łuczak, J. Ranke, J. F. Fern´ andez, A. M¨ uller and J. Th¨ oming, Colloids Surf., A, 2008, 316, 278–284. 80 T. Inoue, H. Ebina, B. Dong and L. Zheng, J. Colloid Interface Sci., 2007, 314, 236–241. 81 J. Łuczak, J. Hupka, J. Th¨ oming and C. Jungnickel, Colloids Surf., A, 2008, 329, 125–133. 82 C. Patrascu, F. Gauffre, F. Nallet, R. Bordes, J. Oberdisse, N. de Lauth-Viguerie and C. Mingotaud, ChemPhysChem, 2006, 7, 99–101. 83 B. Bhargava and M. L. Klein, Mol. Phys., 2009, 107, 393–401. 84 B. Bhargava and M. L. Klein, J. Phys. Chem. A, 2008, 113, 1898–1904. 85 J. N. Israelachvili, Intermolecular and Surface Forces, Academic Press, Revised Third Edition, 2011. 86 M. J. Rosen and J. T. Kunjappu, Surfactants and Interfacial Phenomena, John Wiley & Sons, 2012.
Soft Matter, 2014, 10, 8641–8651 | 8651
