Structure and dynamics of POPC bilayers in water solutions of room temperature ionic liquids.

Structure and dynamics of POPC bilayers in water solutions of room temperature ionic liquids Antonio Benedetto, Richard J. Bingham, and Pietro Ballone Citation: The Journal of Chemical Physics 142, 124706 (2015); doi: 10.1063/1.4915918 View online: http://dx.doi.org/10.1063/1.4915918 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/142/12?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Nanomechanical properties of lipid bilayer: Asymmetric modulation of lateral pressure and surface tension due to protein insertion in one leaflet of a bilayer J. Chem. Phys. 138, 065101 (2013); 10.1063/1.4776764 Diffusion of water and selected atoms in DMPC lipid bilayer membranes J. Chem. Phys. 137, 204910 (2012); 10.1063/1.4767568 Shape transformation of lipid vesicles induced by diffusing macromolecules J. Chem. Phys. 134, 024110 (2011); 10.1063/1.3530069 Simulations of the dynamics of thermal undulations in lipid bilayers in the tensionless state and under stress J. Chem. Phys. 125, 234905 (2006); 10.1063/1.2402919 Structure and dynamics of water at the interface with phospholipid bilayers J. Chem. Phys. 123, 224702 (2005); 10.1063/1.2132277

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.179.16.201 On: Thu, 04 Jun 2015 13:24:40

THE JOURNAL OF CHEMICAL PHYSICS 142, 124706 (2015)

Structure and dynamics of POPC bilayers in water solutions of room temperature ionic liquids Antonio Benedetto,1,2,a) Richard J. Bingham,3 and Pietro Ballone4,5 1

School of Physics, University College Dublin, Dublin 4, Ireland Laboratory for Neutron Scattering and Imaging, Paul Scherrer Institut, 5232 Villigen, Switzerland 3 York Centre for Complex Systems Analysis, University of York, York YO10 5GE, United Kingdom 4 Center for Life Nano Science @Sapienza, Istituto Italiano di Tecnologia (IIT), 00185 Roma, Italy 5 Department of Physics, Università di Roma “La Sapienza,” 00185 Roma, Italy 2

(Received 19 January 2015; accepted 27 February 2015; published online 27 March 2015) Molecular dynamics simulations in the NPT ensemble have been carried out to investigate the effect of two room temperature ionic liquids (RTILs), on stacks of phospholipid bilayers in water. We consider RTIL compounds consisting of chloride ([bmim][Cl]) and hexafluorophosphate ([bmim][PF6]) salts of the 1-buthyl-3-methylimidazolium ([bmim]+) cation, while the phospholipid bilayer is made of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC). Our investigations focus on structural and dynamical properties of phospholipid and water molecules that could be probed by inelastic and quasi-elastic neutron scattering measurements. The results confirm the fast incorporation of [bmim]+ into the lipid phase already observed in previous simulations, driven by the Coulomb attraction of the cation for the most electronegative oxygens in the POPC head group and by sizeable dispersion forces binding the neutral hydrocarbon tails of [bmim]+ and of POPC. The [bmim]+ absorption into the bilayer favours the penetration of water into POPC, causes a slight but systematic thinning of the bilayer, and further stabilises hydrogen bonds at the lipid/water interface that already in pure samples (no RTIL) display a lifetime much longer than in bulk water. On the other hand, the effect of RTILs on the diffusion constant of POPC (DPOPC) does not reveal a clearly identifiable trend, since DPOPC increases upon addition of [bmim][Cl] and decreases in the [bmim][PF6] case. Moreover, because of screening, the electrostatic signature of each bilayer is only moderately affected by the addition of RTIL ions in solution. The analysis of long wavelength fluctuations of the bilayers shows that RTIL sorption causes a general decrease of the lipid/water interfacial tension and bending rigidity, pointing to the destabilizing effect of RTILs on lipid bilayers. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4915918] I. INTRODUCTION

Water molecules at the interface between biomolecules and their physiological environment play a role far exceeding the volume and mass they encompass.1 Protein folding, for instance, often relies on the interplay of hydrophobic and hydrophilic domains.2 Water in close proximity of nucleic acids screens their large electrostatic charge, affecting the structure and properties of DNA and RNA.3 Last but not least, bio-membranes that are made primarily of phospholipids are stabilised and made functional by their interaction with the water environment in which they reside.2 Our study is concerned with this last class of systems, i.e., hydrated phospholipid bilayers. We explore, in particular, the possibility of affecting the properties of phospholipid bilayers in water by the addition of organic electrolyte species belonging to the so-called room temperature ionic liquids (RTILs) family.4 In our study based on molecular dynamics (MD) simulations, we consider bilayers made of 1-palmitoyl-2-oleoylsn-glycero-3-phosphocholine (POPC), whose structural fora)Author to whom correspondence should be addressed. Electronic mail:

[email protected] 0021-9606/2015/142(12)/124706/21/$30.00

mula is shown in Fig. 1. POPC is a neutral (zwitterionic) molecule that represents the most abundant phospholipid in the membrane of animal cells.5 Extensive investigations by experiments and simulations6–8 have made POPC a benchmark and a basic model system in phospholipid and bio-membrane studies. Simple POPC bi-layers and mono-layers deposited on solid conducting surfaces might also find applications in electrochemistry and in nanotechnology.9 Phospholipid bilayers in living cells are solvated by water solutions of electrolyte species,10 whose permeability across the bilayer has major physiological relevance. Crossing the phospholipid barrier typically occurs through transmembrane proteins, whose concentration, distribution, and structure are therefore crucial. Equally important, however, are the properties of the local environment. These, in particular, affect the rate of trans-membrane processes through the formation of the electrostatic double layer and through changes in the water structure and dynamics at the interface. Thus, the control of the quality and concentration of ions in solution provides a suitable way to tune a variety of bilayer or membrane properties through direct lipid-ion coupling, as well as indirect ion-water interactions.10,11 Experiments and simulations have shown that, not surprisingly, the strongest lipid-ion interactions concern ionic lipids,

142, 124706-1

© 2015 AIP Publishing LLC


124706-2

Benedetto, Bingham, and Ballone

FIG. 1. Structural formula of (a) the POPC phospholipid molecule, RTIL ions [bmim]+ (b), and [PF6]− (c). Numbers in part (a) identify the carbonyl oxygens (1 and 2) and the non-ester phosphate oxygens (3 and 4).

while neutral lipids are far less affected.12 However, di-cation species such as Mg++ and Ca++ are able to change the properties of even zwitterionic species such as POPC, mainly because of their small size and high charge.13 In recent years, the class of ionic compounds of interest for the soft-matter research community has swelled to enormous proportions with the inclusion of the so-called RTILs.4 By definition, RTILs are ionic compounds that in their homogeneous phase are liquid below the conventional limit of 100 ◦C. We consider, in particular, two relatively simple systems, based on the prototypical 1-butyl-3-methyl imidazolium cation ([bmim]+) neutralised either by [Cl]− or [PF6]−. Our interest in RTILs parallels that of a large portion of the chemical physics community, fuelled, in turn, by the number and relevance of the applications proposed for these systems.14 RTIL ions, typically bulkier than their inorganic counterparts, interact with bio-molecules through their charge but also through neutral dispersion forces. Moreover, their hydrophilic/hydrophobic character can be controlled by selecting the radius of their ionic moieties, as well as the size of their neutral parts. As a result, their ability to enter and modify the water/bio-molecule interface is enhanced with respect to simpler electrolytes and can be tuned by a suitable choice of cation and anion, opening new vistas for biological manipulations, including a wide range of pharmaceutical applications. On the other hand, the affinity of RTILs for bio-molecules could also cause unexpected effects, with potential harmful consequences on health and on the environment.15 These considerations have motivated an increasing number of investigations, concerning biological effects of RTIL compounds, analysing, in turn, their interactions with a wide variety of bio-molecules, ranging from proteins16 to nucleic acids17 such as DNA and RNA. The broadest impact, however, could arguably be expected from the experimental investigation of lipid-RTIL interactions, since the affinity and permeability of membranes with respect to RTIL ions are the first pre-requisite for an important RTIL effect on cells.18 Ingenious luminescence measurements have revealed the tendency of imidazolium cations to be incorporated in bilayers made of phospholipids such as 1,2-Dioleoyl-sn-glycero-3phosphocholine (DOPC) and 1,2-di[trans-9-octadecenoyl]-snglycero-3-phosphocholine (DEPC),19 which are found in both

J. Chem. Phys. 142, 124706 (2015)

natural and model membranes. The picture is confirmed and completed by quartz microbalance measurements20 and AFM images,21 these last showing nanometer-scale holes in DEPC bilayers resulting from the sorption of 1-octyl-3-methyl imidazolium bis(trifluoromethylsulfonyl)imide ([omim][Tf 2N]). All these measurements point to an increasing effect of RTIL with increased length of the saturated and neutral alkane tail in the 1position of imidazolium, emphasizing the role that dispersion forces play in the RTIL/phospholipid interaction in a water environment. Recent neutron reflectometry measurements22 on POPC and 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) bilayers deposited on a silica substrate and hydrated by either pure water or RTIL ([bmim][Cl] and choline chloride, [Chol][Cl]) water solutions at 0.5M concentration have investigated the ab- and ad-sorption of cations in/on the lipid bilayer. Furthermore, the insertion of RTIL cations into phospholipid monolayers adsorbed on the Hg electrode has been investigated by electrochemical methods,23 finding results compatible with those obtained in other experimental studies on bilayers. These experimental findings have already motivated a few computational studies24–28 that, however, still represent only the first preliminary explorations of a vast new subject. The results of these studies, all concerned with imidazolium-based RTILs, agree on the general picture of easy solvation of the RTIL cation by lipids, driven by a combination of Coulomb and dispersion forces. In agreement with experimental findings, the solubility of [Cn mim]+ into the lipid phase increases with increasing n.27 More specific aspects, concerning the effect of inserting RTIL cations into the plasma membrane of cancerous cells, have been investigated by molecular dynamics in Ref. 29. The present study aims at improving and extending the results of Refs. 24 and 25, using larger systems, and, more importantly, carrying out a more extensive analysis of trajectories, focusing on dynamical properties and hydrogen bond (HB) dynamics before and after the addition of RTILs in solution. We focus, in particular, on aspects that could be probed by neutron reflectometry, elastic, and inelastic neutron scattering experiments that are already programmed. Since these measurements require intact bilayers, homogeneous along their reference plane, we considered short tail imidazolium derivatives that are adsorbed without compromising the basic integrity of the lipid bilayer.19 This choice has the additional advantage of fairly high mobility and relatively short relaxation times for the added species, compatible with the simulation times scales. The fairly large size of our samples allows us to consider multiple bilayer in each simulation cell. First of all, this reflects the setup used in many experimental studies, in which measurements are carried out on stacks of up to 1000 bilayers, separated by water films. The thickness of these lasts, and thus, the periodicity of the stack is determined by equilibration with respect to a water vapour phase of known thermodynamic conditions.30,31 More importantly, this multi-layer configuration allows us to investigate the subtle correlations in the motion of independent lipid bilayers separated by nanometric water films, both in pure systems and in the presence of RTIL molecules.


124706-3


J. Chem. Phys. 142, 124706 (2015)

II. MODEL AND METHOD

The model and method applied in the present study closely follow those used in our previous investigation of similar systems,25 apart from some specific choices suggested by the aims of the present investigation. The force field underlying our simulations has been described in detail in Ref. 25. For this reason, we summarise here only the information required to reproduce our computations, without discussing the choice or the details of the model. The potential energy model that we use consists of a force field taken from the literature for the lipids (POPC) and for the imidazolium-based ionic liquids, together with the well known Extended-Simple Point Charge (SPC/E) model of water.32 In our computations, POPC is modelled at the united atom level, replacing CH2 and CH3 aliphatic groups by a single particle, resulting in the representation of POPC (134 atoms) by only 52 particles. Parameters for bonded and nonbonded interactions have been taken from the Gromos53a6 force field,33 incorporating Berger’s model.34 Atomic charges for particles representing POPC have been taken from Ref. 35. RTIL ions are represented by the all-atom potential of Ref. 36, conforming to the AMBER/OPLS (Assisted Model Building with Energy Refinement/Optimized Potentials for Liquid Simulations) model. To improve the description of the ion dynamics in pure RTIL systems, it has been suggested37 to rescale the atomic charges provided in Ref. 36, reducing the electrostatic charge of each ion to ±0.8e. We retain the original charges, summing up to ±1e per ion, because while the rescaling might be appropriate for bulk RTILs, it is less justified and certainly less tested for ions in an inhomogeneous environment, in which the ion-ion contact is replaced by ionwater and ion-lipid interaction. The most delicate feature in the model concerns the cross interaction between molecules whose force field has been developed independently. In all cases, we adopt the geometric-average Berthelot rule, ϵij = σi j =

√ √

ϵ ii ϵ j j ,

(1)

σii σ j j .

(2)

A few further technical remarks on the matching of the AMBER/OPLS and Gromos force fields used for RTILs and lipids, respectively, are given in Ref. 25. The development of lipid and RTIL force fields is a very active research field, still producing a stream of new and better models. Potentials introduced in the last few years (Ref. 38 for lipids and Ref. 39 for RTILs) might indeed provide a better overall description of the systems that we studied. We retain the same basic model of Ref. 25 because the results published until now show that properties computed with the old and new force fields are still very close, and the difference does not compensate for the loss of an immediate and quantitative comparison with our previous computations. Molecular dynamics simulations in the NPT ensemble have been carried out using the GROMACS package,40 version 4.5.5, running in parallel on the Stokes and Fionn computers of the Irish Centre for High-End Computing (ICHEC). The simulated sample is enclosed into an orthorhombic box

FIG. 2. Lipid and water domains identified on a single frame from our simulations of POPC in pure water at T = 300 K. Points of the regular mesh (see text) are represented by full dots. Red: lipid bilayers; blue: water inter-layers. The lipid and water domains have been shifted by L x /4 (along x) with respect to each other to show the rugged surface in between lipids and water.

periodically repeated in space. Long range (Coulomb) forces are computed through particle mesh Ewald. Short range potentials, including the real-space part of the Ewald sum, are cut at 1.2 nm. The LINCS algorithm41 has been used to keep all bond lengths constant. The integration of the equations of motion has been carried out by the Verlet algorithm, with a time step of 1 fs. An integral part of our approach is represented by the method we used to define the instantaneous interface between lipids and water and to compute the intrinsic profile of all chemical species in the system. A variety of methods have been proposed, see, for example, Refs. 42 and 43. In the present study, we superimpose an orthorhombic 3D grid to our systems, consisting of n x × n y × nz points uniformly spaced in each direction and nearly equally spaced among different directions. The number of non-hydrogen atoms within a distance Rcut from each point is determined, and the mesh point is attributed to the lipid or water phase according to a majority rule applied to the number of neighbours. Ions are treated as guest species and are not considered at this stage. The choice Rcut = 6 Å used throughout our analysis already limits the role of short wavelength fluctuations and ensures that the surfaces dividing lipids from water are regularly spaced. To prevent the formation of isolated pockets of one phase within the other, we apply a second round of smoothing, resorting again to a sort of majority rule, this time applied among mesh points: any point surrounded by less than two neighbours of the same lipid or water type is re-attributed to the other type. We obtain in these way the subdivision of the whole system into two (or four) POPC bilayers (see below) and two (or four) distinct water inter-layers. Each of these domains is simply connected, delimited by single-valued surfaces approximately oriented along the (x, y) plane. A representative example is shown in Fig. 2. The unambiguous identification of these surfaces provides (i) a first basis for attributing a volume to lipid and water phases, (ii) a geometric representation for the instantaneous interface between lipids and water, and (iii) a reference frame to compute intrinsic profiles for all species, and, in particular, for the RTIL cations and anions. For our systems made of ∼150 000 atoms, the analysis of one frame requires some 15 s on a laptop when considering a mesh made of 128 × 64 × 64 or, equivalently, 64 × 64 × 128 points for the two sample geometries detailed below.


124706-4


J. Chem. Phys. 142, 124706 (2015)

FIG. 3. Schematic view of the two sample configurations used in our simulations, obtained by the l × m × n replication of a square bilayer patch consisting of 340 POPC molecules. (a) 2 × 1 × 2 replicas; (b) 1 × 1 × 4 replicas. The gray areas represent POPC domains, red areas represent the inter-layer water films. The orientation of the axes is shown on the side of each sample. In both cases, the z axis is normal to the nominal bilayer planes.

III. RESULTS

We aim at modelling oriented stacks of lipid bilayers that in experiments are usually prepared on a solid planar substrate44 and routinely probed by X-ray diffraction and neutron scattering.7,31 We consider two different sample configurations, the first made of two bilayers and the second one by four bilayers (see Fig. 3). Each sample has been prepared starting from the same 340-lipid POPC bilayer equilibrated in water45 already used for the simulations of Ref. 25. The bilayer, approximatively parallel to the (x, y) plane and covered by a thin film of closely bound water molecules, has been replicated 2 × 1 × 2 and 1 × 1 × 4 times in the x, y, and z directions, respectively, giving in both cases samples of 1360 POPC molecules. These orthorhombic samples, in turn, are periodically replicated in space to mimic macroscopically extended systems. The spacing of the bilayers in the z direction has been set to approach the basic periodicity D measured

in X-ray scattering experiments6 on POPC samples at full hydration,30 i.e., in equilibrium with saturated water vapour, at (P, T) conditions close to the P = 1 atm, T = 300 K of our simulations. The D values reported by different papers, however, are somewhat scattered, mainly because of the steep dependence of D on hydration close to the 100% limit, and for this reason, our choice has a slight margin of uncertainty. We took as a reference the D = 6.4 nm (at P = 1 atm, T = 303 K) given in Ref. 6, measured on multi-lamellar vesicles (MLV) in excess water, and we decreased it to D = 6.0 nm, partly to account for slight differences between MLV and oriented lipid stack6 but mainly to enhance the bilayer-bilayer interaction that is one of the targets of our investigation. This choice of D leaves a gap slightly wider than 2 nm between consecutive bilayers that has been filled with water molecules using the genbox utility that is part of Gromacs. We remark that the 2 nm width of the water interlayer corresponds to the range of hydration forces in phospholipid multilayers, which has been estimated by experiments and simulations.46,47 For the lipids in pure water, this procedure gave an hydration level of n w = 20 water molecules per POPC molecule to be compared with the n w′ = 9.4 water molecules directly hydrating the lipid headgroup6,48 and with the n w = 30 estimated in experiments.6 We remark that the n w difference between simulation and experiments is equally due to the reduction of the periodicity D from 6.4 nm to 6.0 nm and to the discrepancy in the surface area per lipid that is discussed below. The lipid fraction thus accounts for two thirds of the system mass, and after equilibration, it occupies a similar proportion of the system volume. The geometrical parameters of the two pure-bilayer samples (#1 and #2) are listed in Table I. The 2 × 1 × 2 sample (#1), in particular, allows us to probe fairly long wavelength fluctuations in the bilayer plane, enhancing the accuracy of surface tension and bending rigidity estimated from fluctuations, as well as providing a clear picture of the bilayer roughness. The 1 × 1 × 4 sample (#2), on the other hand, offers a less constrained view of correlations in the motion of different bilayers separated by a relatively thin hydration interlayer. Simulations in the NPT ensemble have been carried out at P = 1 atm and T = 300 K, with a few test simulations performed at T = 250 K and T = 225 K. At T = 300 K, samples #1 and #2 have been simulated for slightly more than 100 ns. The first 30 ns have been considered equilibration, and statistics has been collected over the remaining 70 ns.

TABLE I. Samples simulated in the NPT ensemble at T = 300 K, P = 1 atm. N w is the number of water molecules; N+ and N− are the number of cations and anions, respectively. When present, cations are [bmim]+. In samples #3 and #4, anions are [Cl]−, while in samples #5 and #6, they are [PF6]−. Each sample contains Nlip = 1360 POPC molecules, arranged on Nb bilayers. Lengths are in nm and volumes in nm3. Error bars are equal to the standard deviation of values measured every 1 ns. #

Nb

Nw

N+ = N −

⟨L x ⟩

⟨L y ⟩

⟨L z ⟩

⟨V ⟩

1 2 3 4 5 6

2 4 2 4 2 4

26 816 26 816 23 862 24 040 23 263 23 421

0 0 224 224 224 224

20.19 ± 0.05 10.09 ± 0.02 20.23 ± 0.06 10.10 ± 0.02 20.26 ± 0.03 10.12 ± 0.02

10.01 ± 0.02 10.00 ± 0.02 10.02 ± 0.03 10.01 ± 0.02 10.04 ± 0.02 10.03 ± 0.02

12.07 ± 0.06 24.17 ± 0.09 11.93 ± 0.06 23.97 ± 0.09 11.89 ± 0.04 23.86 ± 0.10

2440 ± 1.3 2440 ± 1.6 2419 ± 1.5 2424 ± 1.6 2418 ± 1.5 2423 ± 1.7


124706-5


POPC is a singly unsaturated lipid, whose experimental49 main transition temperature Tm is 270.5 K and our simulations at T = 300 K reproduce the disordered liquid phase (L d ) of the bilayer (see the discussion of diffusion and membrane dynamics reported below). Our choice of [bmim][Cl] and [bmim][PF6] for the RTIL component of our samples has been determined by the fact that these are natural benchmark systems, extensively investigated by experiments and simulations. Moreover, [bmim][Cl] is fully miscible with water;50 [bmim][PF6], instead, is beyond its miscibility limit at the 0.5M concentration of our simulations,51 but the sizeable sorption of [bmim]+ and [PF6]− ions into/onto the POPC bilayers prevents the precipitation of RTIL droplets in the water portion of our samples. In each of the two pure lipid/water samples (#1 and #2) already equilibrated during 30 ns, [bmim][Cl] and [bmim][PF6] RTIL molecules have been inserted into the water layers, giving origin to samples #3–#6. More precisely, we generated random positions within the simulation cell, retaining only those whose minimum distance from lipid atoms or other RTIL ions was more than 6 Å. Individual ions have been inserted at these positions, removing the closest water molecules. To approach the 0.5M concentration that we already used in Ref. 25, we inserted 224 ion pairs both in the case of [bmim][Cl] and [bmim][PF6]. To maintain a nearly constant volume, we decreased somewhat the number of water molecules and, therefore, also the hydration level that however remains in the range 17 ≤ n w ≤ 18 waters per lipid molecule. The size of all the simulated systems is listed in Table I. Samples #3, . . . , #6 prepared in this way have been reequilibrated for 30 ns, and also in this case, statistics has been collected over an additional 70 ns. Visual analysis of trajectories shows that over this time span and for all compositions, the bilayers retain their overall integrity. The first close contact of RTIL ions, or, more precisely, cations, with POPC takes place within 1 ns, i.e., in the early stages of equilibration, and after 2 ns, we identify [bmim]+ ions absorbed into the bilayers. No obvious change in the lipid aggregation state L d has been observed upon RTIL addition, and after equilibration, no apparent drift is detected on the main system properties such as volume, area, and RTIL absorption into POPC. Needless to say, we cannot exclude that the system might still evolve even drastically on significantly longer time scales, from micro- to milli-seconds. The in-plane (x, y) distribution of the RTIL ions absorbed into POPC, solvated in water, or sitting at the lipid/water interface appears to be homogeneous. On average, therefore, the structural and electrostatic properties discussed in Subsection III B are functions of the z coordinate only. A. Area and volume analysis

The average volume and aspect ratios of the six samples considered in our simulation are reported in Table I. All systems have nearly the same volume, as intended. The average cross sectional area per lipid A1 is listed in Table II for all the samples considered in our simulations. The results agree with those of Ref. 25 and are close to those of

J. Chem. Phys. 142, 124706 (2015) TABLE II. Cross sectional area per lipid A1, surface compressibility modulus K A, and isothermal (volume) compressibility χ T of the samples listed in Table I. All quantities estimated from the fluctuations of the simulation box. NPT simulation at T = 300 K and P = 1 atm. We report one digit more than consistent with the statistical error. #

A 1 (nm2)

K A (mN/m)

1010 χ T (Pa−1)

1 2 3 4 5 6

0.594 ± 0.01 0.594 ± 0.01 0.597 ± 0.01 0.598 ± 0.01 0.598 ± 0.01 0.598 ± 0.01

550 ± 3 581 ± 2 619 ± 2 646 ± 3 899 ± 4 942 ± 2

3.03 ± 0.01 3.06 ± 0.02 3.16 ± 0.015 3.19 ± 0.02 3.23 ± 0.02 3.28 ± 0.02

Ref. 27 for the same system, computed with a similar but distinct force field model. We first discuss the results for the pure samples that provide a comparison with experiments and previous simulations. Experimental values for A1 range from 66 Å2 at T = 310 K in Ref. 52 to 62 Å2 at T = 323 K in Ref. 53. One measurement6 reports 68 Å2 at T = 303 K. Despite the spread of data, and as already observed in Ref. 25 discussing simulation results obtained using the same potential, the computational value that we obtain seems to underestimate somewhat the experimental value. Nevertheless, the POPC bilayers in our simulations at T = 300 K are in the liquid disordered phase, as already mentioned and as discussed again below when presenting our results for the diffusion coefficient and the deuterium order parameter SC D of POPC. We note, in passing, that the experimental determination of the area per lipid A1 in multilayers is based on the accurate determination (by diffraction) of the primary repetition length D and of the number n w of hydration water molecules. Then, A1 is computed from6,54 A1 D = 2(Vlip + n w vw ),

(3)

where Vlip and vw are the lipid and water molecule volumes, respectively, determined by independent measurement. Our results, however, show that the water and lipid volumes are not exactly additive (see below), and the deviation from non-additivity might depend on the hydration level. This observation adds a further (although relatively minor) uncertainty to the comparison of experimental and simulation results. Addition of RTIL at 0.5M concentration causes a slight expansion of the cross sectional area. Assuming, to first approximation, that the POPC volume remains constant, the change in cross sectional area should be accompanied by a correspondingly slight (δd ∼ −0.2 Å) thinning of the bilayers, taking place through the mechanism suggested in Ref. 55. This variation agrees qualitatively but only semi-quantitatively with the results of recent neutron reflectometry measurements on the same (POPC + RTIL) and closely related (DMPC + RTIL) systems,22 showing a more substantial thinning of up to 1 Å. The area compressibility can be obtained from a fluctuation expression, derived from the fluctuation-dissipation formulation of linear response theory.56 Following the standard practice, we report its inverse, representing the surface


124706-6


J. Chem. Phys. 142, 124706 (2015)

compressibility modulus, whose definition in terms of thermodynamic quantities is ( ) ∂γ K A = ⟨A⟩ , (4) ∂⟨A⟩ T while in our computations is estimated from the fluctuation relation KA = 

⟨A⟩k BT ⟨A2⟩ − ⟨A⟩2

.

(5)

In the two equations above, A is the cross sectional area covered by the lipid bilayer and γ is the interfacial tension, with averages and partial derivatives computed at constant temperature. The results from our simulations are listed in Table II. The experimental surface compressibility modulus of 12 diacyl-PC lipid bilayers of different chain un-saturation has been determined by micropipette pressurisation of giant micelles.57 The results are centred on an average of 243 mN/m that is somewhat less than 50% of the simulation result. This tells us that the computational model underestimates the fluctuations in the bilayer cross sectional area, and this might be consistent with the underestimation by our model of the area per lipid. The difference, however, could also be due to other reasons, such as the different hydration level, or the fact that the experiment of Ref. 57 probes single lipid bilayers in bulk water, curved into closed, nearly spherical micelles. Our simulations, instead, concern stacks of planar bilayers interacting with each other at close range. Moreover, we did not apply the correction for undulation discussed in Ref. 58, because we think that the overall (model + simulation) uncertainty on this quantity is still larger than the undulation correction. The effect of this correction would be to decrease somewhat the computational estimate of K A reported in Table II. Turning now to volume properties of the simulated samples, we first discuss here bulk properties of the composite water/lipid samples, even though these quantities are seldom reported or discussed in experimental papers. The isothermal compressibility of each sample is proportional to its volume fluctuations,56,59 according to χT =

⟨V 2⟩ − ⟨V ⟩2 , k BT⟨V ⟩

(6)

where V is the fluctuating system volume and k B is the Boltzmann constant. The results for χT are reported in Table II. We observe that the systems with only two bilayers are slightly but systematically less compressible than those with four bilayers. In other terms, and perhaps not surprisingly, the volume of the sample with four bilayers fluctuates slightly more than that of the sample of equal composition and only two independent bilayers. This slight difference between systems that are nominally equivalent apparently is the effect of correlations in the volume fluctuations of neighbouring lipid bilayers. Compared to the compressibility of other molecular liquids, the isothermal compressibility of bulk water is relatively low.60 At T = 300 K and ambient pressure, the experimental value for water is χT = 4.53 × 10−10 Pa−1, reproduced almost

exactly ( χT = 4.63 × 10−10 Pa−1 at T = 298 K, p = 1 atm) by the SPC/E model.61 Incidentally, and apparently fortuitously, the experimental compressibility of water is comparable to that of room temperature ionic liquids62 such as [Cn mim][PF6] (n = 4 and 8) or [Cn mim][BF4] (n = 8), whose cohesive energy per molecule, however, is an order of magnitude higher than that of water. Remarkably, the compressibility of the layered lipid/water system, reported in Table II, is significantly lower than that of water and of homogeneous RTILs. The reasons for this low compressibility of layered lipid/water systems are discussed below. The compressibility of each sample increases moderately but systematically upon addition of the RTILs. This is consistent with the idea that the inclusion of RTILs somewhat decreases the overall stability of the ordered bilayer stacks. The relation between the volume isothermal compressibility χT of the lipid/water composite sample and the surface compressibility modulus K A of the lipid bilayers is not so obvious because of correlations in the fluctuations along the directions parallel and perpendicular to the bilayer. As a result, we observe a slight increase of K A upon adding the RTIL, implying lower surface fluctuations, while the volume compressibility also increases slightly, pointing however to a softer system and larger volume fluctuations. To characterise the anisotropy in the system response to external perturbations, we consider the fluctuations of the length L α of the sides α = x, y, and z. Assuming that the sample behaves like an elastic spring in each direction, we have 1 Eα ⟨Aα ⟩⟨(L α − ⟨L α ⟩)2⟩ 1 k BT = , 2 2 ⟨L α ⟩

(7)

where ⟨Aα ⟩ is the average area of the simulation box perpendicular to the α direction, and Eα is Young’s modulus in the α direction. The results are summarised in Table 1 of the supplementary material.94 We find that for α = x or y, Eα is of the order of 0.4 GPa, i.e., comparable to that of low-density polyethylene. In all samples, Ez is ∼1/3 of either E x or E y . In other words, again as expected, the system is softer and fluctuates more along the direction perpendicular to the bilayer than in the bilayer plane. We observe that volume fluctuations (Eq. (6)) are much less than predicted on the basis of the fluctuations of the sides of the simulation box (Eq. (7)), the difference being accounted for by correlations in the fluctuations of different sides. In other terms, a fluctuation increasing L α is likely to be accompanied by a corresponding decrease of L β × Lγ ( β, γ , α), reflecting the fact that it is energetically less costly to change the aspect ratio of the simulation box than to change its volume. This observation is not surprising, and it is generally true for most systems, but we point it out to emphasise the role of correlations among fluctuations as an important factor in determining the response to external perturbations of a composite system such as the water-lipid samples that we simulate. Addition of [bmim][Cl] or [bmim][PF6] increases slightly but systematically all the Eα ’s. Since these quantities are a measure of the system rigidity, their increase may seem at odds with the parallel increase of the compressibility. A


124706-7


J. Chem. Phys. 142, 124706 (2015)

detailed analysis shows that RTILs modify the correlation among fluctuations of the simulation box sides, changing compressibility and Young’s moduli in opposite directions. The analysis of equilibrium volume and volume fluctuations can be refined to consider the contribution of independent components, such as lipid, water, and RTIL molecules, provided we specify a rule to assign a volume to every species. In experiments, the determination of the volume of lipid molecules in lipid/water systems is based on accurate measurements of density and composition. The lipid and water contributions, then, are determined assuming that water molecules retain their bulk volume.54 As a reference, the volume vw of the water molecule at the same thermodynamic conditions has been determined using the same SPC/E model in an independent computation, giving vw = (2.997 ± 0.002) × 10−2 nm3 (already quoted in Ref. 25). More precisely, this value is obtained from the NPT simulation of a bulk water sample, and vw is given by the average volume divided the number of molecules (Nw = 8000) in the system. Assigning a reference volume to POPC is more uncertain. Pure POPC samples, in fact, might not adopt a layered structure at equilibrium, and, at ambient conditions, POPC samples are unlikely to be in the same fluid state of hydrated bilayers. In experiments, even nominally anhydrous lipid samples contain a few percent water in weight. To approach the experimental reality as far as possible, we progressively removed water from samples #1 and #2 until reaching in both cases a composition of 1360 POPC and 3650 water molecules. The samples retained their layered structure during simulations covering 10 ns, but appear to be glassy rather than fluid, possibly reflecting a hindered transition to an ordered (gel) phase. The contribution of the residual water to the average volume has been subtracted out, giving an estimate of Vlip = 1.222 ± 0.002 nm3 per POPC molecule. These reference molecular volumes allow us to verify the volume additivity assumptions used in experiments and in previous simulations to disentangle the lipid and water volumes in ordered stacks of hydrated bilayers. To this aim, we compute the volume of each chemical species in the simulation samples by resorting to the statistical Voronoi construction briefly discussed in the Appendix. The results for the molecular volume of water, POPC, and RTIL ion pairs are listed in Table III. The data for the pure water/POPC samples show that the volume per water molecules is (2.951 ± 0.02) × 10−2 nm3, i.e., 1.5% less than the reference volume of a bulk water

molecule. At the same time, the molecular volume of POPC in the hydrated bilayers turns out to be vl = 1.212 ± 0.003 nm3, i.e., nearly 1% less than the reference anhydrous value. These differences in molecular volume can be interpreted in terms of a negative non-additivity of the water and POPC volumes upon mixing. Negative non-additivity, in turn, reflects the high affinity of water and lipids that stabilises the bilayer stacks. Admittedly, the relative variation of molecular volumes is small, but in our pure samples (#1 and #2), they correspond to a negative mixing volume of ∼3 nm3, equal to the volume of 100 water molecules. Notice that although expressed in volume units, the non-additivity is in fact an interface effect, that, therefore, might not be completely negligible for stacks of closely packed lipid bilayers separated by nanometric water films, and perhaps, it could be taken into account in the interpretation of experimental data. The results for #3, . . . , #6 reported in Table III show that the inclusion of RTIL reduces somewhat (by ∼30%) the size of the lipid-water non-additivity, pointing to a lesser cohesion of these two basic components of our samples. Having identified the volume of lipids and water, we can now discuss the compressibility of each component. Applying the same fluctuation formula separately to water and lipid volumes, we obtain the isothermal compressibility values listed in Table III. The compressibility of the lipid fraction is close to that of similar phospholipids63 and of saturated or nearly saturated hydrocarbon systems.64 Somewhat surprisingly, the compressibility of the water fraction in sample #1 and sample #2 is between 1/3 and 1/2 of the compressibility of bulk SPC/E water, closer to the isothermal compressibility of ice at T ≤ 273 K (equal to 1.3 × 10−10 Pa−1)65 than to the value of the liquid phase at ambient conditions. Such a low compressibility could be due to the slow diffusive dynamics of water molecules and hydrogen bonds at the water/lipid interface that will be discussed in Subsection III D. A comparable decrease of water compressibility in proximity of proteins has been found and discussed in Ref. 66. To first approximation, the isothermal compressibility of the composite system (see Table II) is indeed the average of the lipid and water compressibilities, weighted by their relative volume fraction. In other terms, correlations between volume fluctuations of POPC and water, in this case, are not quantitatively relevant. Adding RTILs to the solution increases the compressibility of both the lipid and water fractions in all samples. The relative variation is somewhat larger for water than for POPC and is almost the same for samples added with [bmim][Cl] and

TABLE III. Volume per molecule (nm3) and isothermal compressibility (Pa−1) of water (v w , χ w ), POPC (Vlip, χ lip), and RTIL ion pairs (VRTIL, χ RTIL), respectively. All quantities estimated from the results of the statistical Voronoi construction described in the Appendix. The error bar on compressibilities is of the order of 5%. #

v w (nm3)

1010 × χ w

Vlip (nm)3

1010 × χ lip

VRTIL

1010 × χ RTIL

1 2 3 4 5 6

0.029 51 ± 0.000 02 0.029 50 ± 0.000 02 0.029 78 ± 0.000 02 0.029 77 ± 0.000 02 0.029 68 ± 0.000 02 0.029 67 ± 0.000 02

1.78 1.81 2.22 2.29 2.25 2.29

1.212 ± 0.002 1.212 ± 0.002 1.215 ± 0.002 1.215 ± 0.002 1.216 ± 0.002 1.216 ± 0.002

3.41 3.44 3.95 3.99 4.08 4.15

0.2469 ± 0.001 0.2461 ± 0.001 0.3409 ± 0.001 0.3410 ± 0.001

1.30 1.33 1.52 1.58


124706-8


J. Chem. Phys. 142, 124706 (2015)

[bmim][PF6]. The effect of the RTIL addition is stronger on the individual compressibilities of water and POPC than on the volume compressibility of the whole samples, the difference being due, apparently, to correlations in the fluctuations of POPC and water molecular volumes. Compressibility of the RTIL ions themselves can be obtained by the same fluctuation approach used for water and lipids. The resulting χT ∼ 1.5 × 10−10 Pa−1 is roughly one half the value measured in their bulk liquid phase (reduced to account for the 2D diffusion in our systems), apparently because of the tight coordination by their first hydration shell. Compressibility of water and of lipids in the samples with four layers is slightly larger than that of the samples with two layers, pointing again to the effect of weak but systematic correlations among the fluctuations of different water and lipid layers. B. Atomic and charge density profile

We now turn to the analysis of the equilibrium structure of the lipid bilayers, water films, and ionic species, when present. The first apparent result is the fact that our bilayers, whose maximum cross sectional area reaches 20 × 10 nm2 in the 2 × 1 × 2 samples, are stable during the entire simulation covering about 100 ns, including equilibration. To be precise, the density profile of non-hydrogen atoms67 for pure samples, shown in Fig. 4 for sample #1, is fairly broad. A widely accepted parameter used in the statistical mechanics to quantify the width of monotonic density profiles at surfaces and interfaces is the so called 90%-10% width, measuring the separation between the two planes at which the density profile of interest reaches 90% and 10% of the bulk value.68 In the case of pure sample #1, we estimate a fairly high 90%-10% width of 8.0 Å for lipids and 11.2 Å for water. Moreover, while the density distributions of different water inter-layers are always well separated, the lipid profiles of adjacent bilayers are slightly superimposed, almost suggesting that their identity is compromised. As already anticipated, the analysis of configurations shows that this is not the case: although the

FIG. 4. Density profile of non-hydrogen atoms along the direction z perpendicular to the POPC reference plane. Neat sample #1. The distribution of P atoms is a subset of the lipid distribution. It is shown here to emphasise the compact localisation of the polar POPC heads along z.

FIG. 5. Comparison of cation and anion density profiles in sample #3 (blue curves) containing [bmim][Cl] and in sample #5 (red curves) containing [bmim][PF6]. The [bmim]+ curves report the density of non-hydrogen atoms in the cation. The [PF6]− curve refers to the distribution of its central P atom only. For the sake of clarity, both the [Cl]− and [PF6]− curves have been multiplied by 7, and only the two central interfaces have been included in the plot. The light shaded area corresponds to the density profile of lipid atoms.

lipid and water density profiles are, as expected, blurred by large-amplitude long-wavelength oscillations, correlations in the motion of neighbouring surfaces keep the bilayers well separated and always hydrated by a water inter-lipid layer of nearly constant width. The stability of the bilayers and their subdivision into well defined leaflets are confirmed by the plot of the density profile for the P-atom in the phosphatidylcholine head of POPC, showing well localised peaks, whose full width at half maximum is only 6 Å. The density profile of individual bilayers obtained for the 2 × 1 × 2 (#1) and 1 × 1 × 4 (#2) samples can be superimposed almost exactly. The addition of RTIL ions changes slightly the lipid and water profiles and, of course, introduces new species into the picture. The detailed analysis of the relative position of RTIL ions with respect to lipids and water requires the computation of intrinsic density (or concentration) profiles that are defined and discussed below. Absolute density profiles, however, confirm what had already been seen in previous computations: cations enter the lipid bilayer, giving origin to a sizeable accumulation just within the lipid side of the interface (see the results for [bmim][Cl] in Fig. 5). The [Cl]− anion does not follow [bmim]+ into the POPC bilayer, while the [PF6]− profile is broader and extends several Å inside the lipid layer (see Fig. 5), in proximity of the POPC head. The enhanced affinity of [PF6]− for POPC is likely to be due to the combination of dispersion interactions, which increase with increasing size of the anions, and of hydration free energies, which decrease with increasing ionic radius, making it easier for [PF6]− than for [Cl]− to enter the lipid range. Setting a cutoff distance of 3.5 Å between POPC oxygens and the centre of mass and charge of the imidazolium ring, we find a sizeable association of imidazolium with the most ionic oxygens of POPC, consisting of the carbonyl oxygens 2 and especially 1 (see Fig. 1), and of the non-ester oxygens (3 and


124706-9


FIG. 6. Representative configuration of POPC and [bmim]+ from the [bmim][Cl] simulation. The green segment joins the O2 atom on POPC to the centre of the imidazolium ring of the [bmim]+ ion. The distance from the centre of mass of the imidazolium ring and one of the carbonyl oxygens of POPC is 3.2 Å.

4 in Fig. 1). On average, 30% of the [bmim]+ ions show such a short range association with POPC oxygens in the [bmim][Cl] case. The preferred docking is at the carbonyl oxygens, but also the non-ester oxygens of phosphonium show a non-negligible coordination by imidazolium. No [Cl]− is found within such a short (i.e., 3.5 Å) distance from POPC oxygens because of the mutual repulsion between negatively charged atoms. The proportion of [bmim]+ closely associated to POPC oxygens increases to 35% in the [bmim][PF6] solution, and in this case, a few [PF6]− ions fall, perhaps accidentally, within the 3.5 Å cutoff distance from POPC oxygens, confirming the somewhat higher affinity of this larger anion with POPC. The tight association of cations with oxygen atoms in the POPC polar head is reminiscent of a similar result for traditional salts like NaCl,69 but the preference for carbonyl oxygens ([bmim]+) or phosphate oxygens ([Na]+) is opposite in the two cases. The relative position and orientation of POPC and [bmim]+ arising from the interplay of dispersion and Coulomb forces has been described in Refs. 25–27 and is confirmed by the present simulations (see Fig. 6). Further analysis of snapshots shows that [bmim]+ ions absorbed just below the POPC/water interface retain at least some of the water molecules hydrating their ionic five-fold ring upon entering the lipid bilayer, thus avoiding the large free energy barrier of fully shedding the solvation shell of their most ionic moiety. Snapshots of hydration shells around POPC and [bmim]+ molecules are given as supplementary material.94 Understanding the equilibrium and dynamical properties of dipole polarisation at the lipid/water interface is important because it provides insight into the origin of the hydration force between bilayers that has been discussed in a number of paper (see Ref. 70 for a short review). Because of the intrinsic anisotropy of the interface, the orientation of the water dipole along the direction perpendicular to the bilayer is expected and observed in our samples (see Fig. [S1] in the supplementary material94 for sample #1). Our results on the water orientation agree well with those of previous simulations of the same and similar systems.71 At each interface, in particular, the net water dipole points away from the bilayer, and the peak of the dipole density profile along z is equivalent to the full dipole orientation of 8% of the water molecules at the interface. Besides water, POPC also has a non-negligible dipole moment that shows a net polarisation along z. In each leaflet,

J. Chem. Phys. 142, 124706 (2015)

the POPC dipole moment points towards the inside of the bilayer. Thus, POPC and water dipoles orient themselves in opposite directions and to a large extent compensate each other at each interface, thus avoiding the large self-energy of a planar distribution of dipole polarisation. Notice that at short range, Coulomb forces favour the opposite arrangement of dipoles, corresponding to a +-/+- layering of charges at each lipid/water interface, giving origin to a large net surface dipole that, at extended planar interfaces, is opposed by the long range part of the Coulomb energy. It would be interesting to verify whether water polarisation around lipid aggregates of nanometric diameter such as micelles is opposite from what found at extended planar bilayers. Short range dispersion forces and entropy effects could also contribute to the dipole polarisation of water and lipids at their interface, but we do not have a quantitative estimate of their effect. Analysis of simulation trajectories confirms the strong correlation in the dipole moment of all system components, apparently imposed by the domineering role of Coulomb forces in extended systems: (i) the dipole of the two POPC leaflets in each bilayer points in opposite directions and nearly cancels each other; (ii) as already mentioned, the same compensation of dipole moments occurs at the lipid-water interface; (iii) the simulations on samples with four independent bilayers show that even in this case all dipole moments are closely coupled, and correlations decay only moderately over the longest inter-bilayer distance (L z /2 ∼ 12 nm) covered by our simulations. Because of the [bmim]+ penetration into POPC, the addition of RTIL ions gives origin to a new double layer at each lipid/water interface, whose dipole is oriented from water to the inside of the lipid bilayer. The size of this dipole is not easy to assess independently because of overlapping charge distributions and large fluctuations, but its presence is made apparent by a 30% increase of the average size of the POPC dipole at each leaflet. At the same time, the peak dipole polarisation of water decreases by ∼10%, both in the [bmim][Cl] (dashed curve in Fig. [S1] of the supplementary material,94 referring to the [bmim][Cl] sample #3) and in the [bmim][PF6] cases. Since the water density profile is nearly the same in all samples, irrespective of the solvation of RTILs, the decrease of polarisation is only due to a weaker orientation of water molecules, pointing again to a (slightly) weaker coupling of lipids and water. To quantify the strength of Coulomb interactions across our samples, we re-visit the same electrostatic aspects analysed in the previous paragraphs, now referring to the charge and electrostatic potential distributions, instead of focusing on the double layers and dipoles at the interfaces. Even before the addition of dissociated species such as [bmim][Cl] and [bmim][PF6], the system displays a significant (bound) charge distribution ρQ(z), arising from the z dependence of the dipole polarisation density, according to ρQ(z) = −∇P(z).

(8)

As pointed out in Ref. 25, the water and lipid contributions to the charge density nearly but not exactly cancel each other (see Fig. [S2] of the supplementary material94 for the change density profile of the pure sample #1). Integration of the


124706-10


J. Chem. Phys. 142, 124706 (2015)

Poisson equation along z, 4π d 2ϕ = − ρQ(z), 2 ϵ0 dz

(9)

gives the electrostatic potential profile across the bilayers (shown in panel (b) of Fig. [S2] in the supplementary material94). Inside the POPC bilayer, the electrostatic potential is 0.25 V higher than in the water interlayer, and a positive unit charge (e) crossing the water to lipid interface has to overcome an electrostatic potential energy barrier of 0.46 eV. The results of the present simulations are compatible but quantitatively different from those of Ref. 25 and also from previous computations.12 In particular, both the potential energy difference and the electrostatic barrier are ∼30% less than in previous computations. The major reason for the discrepancy is certainly the lower hydration level, causing the interpenetration and partial neutralisation of electrostatic double layers. A second reason is the large cross section of our samples that blurs sharp features of the charge density profile because of long wavelength capillary oscillations. The analysis of instantaneous values and of fluctuations in the strength of electrostatic properties, not carried out in the present study, could give a far more vivid picture of the role of Coulomb interactions in these systems and could be experimentally probed by optical spectroscopy on suitable chromophore species. The addition of compounds such as [bmim][Cl] and [bmim][PF6] that dissociate in water gives origin to free charges that partially separate at the lipid/water interfaces.

FIG. 7. (a) Density profile ρ Q (z) of bound (blue, dashed line) and free charge (red, dashed-dotted line) in sample #3. The total charge density is represented by the full, black line. (b) Comparison of the electrostatic potential profile in sample #1 (no RTIL, dashed line) and in sample #3, containing [bmim][Cl]. The zero of the potential is conventional. The lipid density profile is reported in arbitrary units to locate the bilayer and water regions.

FIG. 8. Main plot: intrinsic density profiles for sample #3. Green line: POPC; blue: water; red: [bmim]+; black: [Cl]−. Only non-hydrogen atoms are counted, and for the sake of clarity, the [Cl]− profile has been multiplied by 10. Inset: intrinsic density profile for water in sample #1 (black curve) and in sample #3 (blue line). The difference, highlighted in red, points to a water excess in the sample #3 case, located at the absorption site of the ionic ring of [bmim]+.

This introduces a further and sizeable electrostatic double layer at each interface, whose amplitude is, in fact, larger than that of the bound charge in pure samples (see panel (a) in Fig. 7). The energy density of Coulomb interactions, however, is such that any perturbation of the charge density is effectively opposed by screening, and the net result is that the charge and electrostatic potential profiles along z are not drastically changed by the RTIL addition (panel (b) in Fig. 7), since the individual water and lipid contributions change significantly to compensate the effect of the free ions. In the case of the [bmim][Cl] solution (sample #3, see panel (b) in Fig. 7), for instance, the electrostatic potential inside the bilayer is 0.41 V higher than in the electrolyte solution, and the corresponding electrostatic potential energy barrier is 0.56 eV. The trend is similar to the one observed in the simulations of Ref. 12, following the addition of [Na][Cl] to the hydrating water of POPC bilayers. To compensate the effect of long-wavelength fluctuations of the POPC bilayers that blur the density profiles of lipids, water, and ions, we resort to the definition of intrinsic density profiles that are proportional to the probability that an atom of coordinates (X,Y, Z) is at a distance z i from the geometric ˆ ˆ surface Z(X,Y ) dividing lipids from water. Two such Z(X,Y ) surfaces are defined for each bilayer. The construction of ˆ Z(X,Y ) for every atomic configuration is described in Sec. II, and in our definition, the distance z i is measured along the z direction, i.e., it is not the standard distance of a point from a surface that usually is measured along the normal to the surface. Since bilayers remain nearly planar


124706-11


throughout our simulations, the distinction is not quantitatively important. Intrinsic profiles for the different chemical species are shown in Fig. 8 for sample #3, while for sample #5 is given in the supplementary material.94 For the sake of our exposition, the profiles computed with respect to the z = 0 plane of the simulation space are referred to as Cartesian profiles. Admittedly, the smoothing procedure discussed in Sec. II, together with the choice of the smoothing radius Rc = 6 Å, affects the result at short distance from the Zˆ surface, giving again fairly broad profiles. However, reducing Rc significantly, or removing the smoothing altogether, produces spikes at the interface that are likely to be spurious. Since we do not want to pick the Rc value that gives profiles that we like, we keep the safe Rc = 6 Å that later is used to analyse the membranelike dynamics of the interface. With this choice, the water profile is still monotonic, with a 90%-10% width of 8 Å. The lipid profile instead shows a characteristic structure, with a shoulder at ∼2.5 Å below the reference surface (i.e., on the lipid side), corresponding to the location of the phosphate group, and a sharper peak at ∼6 Å below the surface, around the average position of the POPC carbonyl groups. No statistically significant difference is detectable between the 2 (samples #1, #3, and #5) and the 4 bilayer (samples #2, #4, and #6) cases. Small but statistically significant differences are instead observed in the lipid and especially water profiles between pure [bmim][Cl] and [bmim][PF6] samples, as illustrated in the inset of Fig. 8 for sample #3 and sample #1. The difference points to an increase of water density a few Å below the lipid/water interface. The location of this excess water suggests that it might represent the water molecules hydrating the imidazolium moiety of [bmim]+, driven into the lipid layer by the penetration of [bmim]+. The intrinsic profiles for cations and anions shown in Fig. 9 are qualitatively similar to the Cartesian ones (see Fig. 5), but the double hump of the [bmim]+ profile is more apparent in the intrinsic than in the Cartesian profile. Integration of the profile from the surface to the inside of

J. Chem. Phys. 142, 124706 (2015)

FIG. 10. Distance dependence of the height-height correlation function along the POPC/water interfaces in sample #1 (black, full circle), sample #3 (red, empty circle), and sample #5 (blue, filled squares). The radial variable R is reported on a logarithmic scale.

the POPC bilayers confirms that roughly 1/3 of the cations are absorbed within the bilayer, both in the [bmim][Cl] and [bmim][PF6] samples, with a slight excess in the latter case. Acyl chain order parameters (SC D ),72 related to the degree of ordering of the POPC tails, and relevant for the interpretation of NMR data, show the same features seen in previous computations (see Refs. 25 and 27). For this reason, these parameters are not shown here, but given in the supplementary material.94 The first observation concerning these parameters, however, is that their values around −SC D ∼ 0.2 are fully consistent with the disordered liquid phase (L d ) for the bilayer. On each acyl chain, addition of either [bmim][Cl] or [bmim][PF6] has a relatively minor effect on SC D , the changes being primarily localised on the first few carbon atoms closest to the two carbonate groups. The absorption of [bmim]+ into POPC, in particular, has the effect of increasing the chain ordering at these first carbon atoms, as measured by a slight but nevertheless statistically significant increase of their −SC D parameter. This trend upon RTIL addition agrees very well with the results of Ref. 27. The height-height correlation function along the POPC/ water interface is defined as h(R) = ⟨[z(R + R0) − z(R0)]2⟩R0,

FIG. 9. Intrinsic density profiles of cations and anions in the 2 × 1 × 2 samples #3 and #5. Only non-hydrogen atoms are counted, and for the sake of clarity, the [Cl]− profile has been multiplied by 10. The vertical dashed line separates water (z > zsurf ) from the lipid range (z < zsurf ).

(10)

and for each interface, it has been computed by averaging over R0 on a grid of nx × n y points73 and over 70 independent configurations. The distance dependence of h(R) has been used to characterise the smooth/rough state of the interface. The results for the pure sample #1 are shown in Fig. 10. According to the accepted statistical mechanics definition,74 the POPC/water interfaces appear to be rough, since on a semilogarithmic plot h(R) shows a non-vanishing growth ∝ log(R) at the longest distances accessible to our simulations. This result is consistent with the fluid-like state of the bilayer, implying, in turn, capillary wave excitations at long wavelength.


124706-12


Adding RTIL ions increases the amplitude of the logarithmic divergence, in agreement with the general result of decreasing cohesion and stability of the lipid bilayer. The effect is stronger for [bmim][PF6] than for [bmim][Cl], consistent with the higher accumulation of ions at the interface in the [bmim][PF6] case. C. Hydrogen bonding at the phospholipid/water interface

HBs have been identified in terms of purely geometric criteria.75 In our definition, triplets of O–H–O atoms form a hydrogen bond whenever the distance between the two oxygens is less than 3.5 Å, and the O–H–O bending angle deviates less than 30◦ from linearity. Both in pure and in RTIL-added samples, water is the only proton-donating species. POPC can accept hydrogen bonds at its two carbonyl oxygen atoms (oxygens 1 and 2 in Fig. 1) and at the two non-ester phosphate oxygens (3 and 4). The other four POPC oxygen atoms (ester oxygens) are not suitable for forming HB. Let us first summarise the results for pure bilayers. Each of the non-ester phosphate oxygen atoms (3 and 4 in Fig. 1) in POPC accepts, on average, 1.38 ± 0.01 hydrogen bonds from water, providing an important contribution to the bilayer stability. The carbonyl oxygen on the oleoyl tail of POPC (oxygen 1 in Fig. 1) accepts virtually the same number (1.39 ± 0.01) of protons from water, while the carbonyl oxygen 2 on the palmitoyl tail is far less well coordinated, accepting, on average, only 0.39 ± 0.02 water protons. The reduced hydrogen bonding of this carbonyl atom is at least partly due to its location on the POPC molecule, somewhat buried below the head groups, with limited contact with water. It is tempting to speculate that the limited water solvation of this highly polar oxygen atom provides an important driving force for the absorption of the imidazolium ionic head of [bmim]+ into POPC. Water molecules share 1.272 ± 0.002 HB among themselves. In other terms, each water molecule on average donates and accepts 1.272 ± 0.002 protons to/from neighbouring water molecules. In this way, at any given time each water shares, i.e., accepts or donates, 2.543 ± 0.004 HB. This number is lower but still comparable to the instantaneous number (∼2.8) of HB shared by each water molecule in bulk water, estimated by ab-initio computations.76 The difference is partly due to the water confinement in between the POPC bilayers that, moreover, accept a non-negligible number of water protons. Taking these HB into account, we find that 75% of the water hydrogen atoms are engaged in hydrogen bonds. Our geometrical criteria attribute a few HB also to the ester oxygens of POPC, but these are probably due to accidental proximity and alignment of atoms in a fluctuating liquid environment. These results are rather insensitive to the choice of the oxygen-oxygen cutoff distance (3.5 Å) but depend somewhat more sensitively on the limiting value for the O–H–O angle. In principle, the NH+3 apex of POPC could also give origin to HB by donating one of its protons. Our united atom model of POPC, however, prevents a detailed analysis of this possibility.

J. Chem. Phys. 142, 124706 (2015)

Despite the low number of ions added to the system, the presence of either [bmim][Cl] or [bmim][PF6] causes characteristic changes in the HB of POPC molecules. To be precise, the coordination of oxygens 1, 3, and 4 increases by only ∼3%, barely significant compared to the 1% error bar. The carbonyl oxygen 2, however, increases substantially the number of hydrogen bonds that it accepts from water, going from the 0.39 ± 0.02 of the pure samples to the 0.48 ± 0.01 of the [bmim][Cl]-added samples and to the 0.51 ± 0.01 of the [bmim][PF6]-added samples. The most likely explanation of this increase is that [bmim]+ enters the bilayer without completely shedding its hydration waters, as already suggested by the results for the intrinsic profile illustrated in Fig. 8. This effect, slightly enhancing the solvent penetration into POPC, increases POPC-water H-bonding, especially at the carbonyl oxygen 2, located on the palmitoyl tail. At the same time, we observe a decrease of the water-water hydrogen bonding from 2.543 ± 0.004 to 2.368 ± 0.004 HB shared among water molecules. The primary reason of this effect, however, is the lower hydration of the RTIL-added samples that has been introduced to keep nearly fixed the system volume at the value of the pure samples. As already stated, the RTIL ions considered in our computations do not form conventional hydrogen bonds. However, the C R carbon of imidazolium (see Fig. 1) is known to be fairly acidic, and in some cases, weak hydrogen bonds centred on this atom have been revealed in the crystal structure of RTIL compounds based on imidazolium.77 Our analysis, instead, reveals virtually no such bonding with the geometric parameters listed in Sec. II. Hydrogen bonds, however, are detected if the geometrical parameters are somewhat relaxed. A cutoff r cut = 3.6 instead of 3.5, for instance, already reveals a small number of long hydrogen bonds connecting C R to water and also to POPC, with a preference for the carbonyl oxygen atoms. The detailed analysis of configurations, however, shows that the proximity of [bmim]+ with POPC oxygens appears to be driven more by purely Coulomb interactions than by bona fide hydrogen bonding that entails a marked directionality of the bonding, and also in this case, the condition on the O–H–O angle being wider than 150◦ is far more selective than the choice of the cutoff distance. Hydrogen bonding between the C R carbon of [bmim]+ and [Cl]− or PF−6 is virtually absent even with the relaxed geometric conditions, while we observe a non-negligible donation of water protons to either [Cl]− or [PF6]−. Distinguishing hydrogen bonding from simple water polarisation by the anions is somewhat ambiguous, but both these effects contribute to the slightly reduced water-water hydrogen bonding. D. Dynamics of molecules and hydrogen bonds

The first dynamical property that we characterise is the diffusion coefficient of lipid and water before and after the addition of RTIL. Because of the anisotropic nature of each sample, diffusion here refers to the motion of molecules parallel to the (x, y) plane, while the motion along z is discussed independently. The diffusion coefficient of molecules of species α, therefore, is estimated from the 2D Einstein relation


124706-13


J. Chem. Phys. 142, 124706 (2015)

TABLE IV. Top half: representative values of the diffusion constant of POPC, water, and RTIL ions in samples #1, #3, and #5. Bottom half: diffusion constant of water and RTIL ions in homogeneous (3D) solutions, without POPC. POPC bilayers in RTIL/water solutions (0.5M) Sample #1 Species D α (cm2/s)

POPC 2.0 ± 0.2 × 10−8

Water 1.20 ± 0.05 × 10−5

Sample #3 Species D α (cm2/s)

POPC 2.8 ± 0.2 × 10−8

Water 0.84 ± 0.04 × 10−5

[bmim]+ 10.1 ± 0.2 × 10−7

[Cl]− 0.66 ± 0.02 × 10−5

Sample #5 Species D α (cm2/s)

POPC 1.5 ± 0.3 × 10−8

Water 0.96 ± 0.04 × 10−5

[bmim]+ 8.1 ± 0.2 × 10−7

[PF6]− 0.19 ± 0.02 × 10−5

Homogeneous RTIL/water solutions (0.5M), no POPC [bmim][Cl]/water D α (cm2/s) [bmim][PF6]/water D α (cm2/s)

Water 2.29 ± 0.03 × 10−5 Water 2.52 ± 0.03 × 10−5

Dα = lim ∆2α (t) = lim t→ ∞

t→ ∞

[bmim]+ 8.2 ± 0.08 × 10−6 [bmim]+ 7.25 ± 0.1 × 10−6

[Cl]− 1.40 ± 0.02 × 10−5 [PF6]− 1.26 ± 0.02 × 10−5

Nα ⟨[X i (t + t 0) − X i (t 0)]2 + [Yi (t + t 0) − Yi (t 0)]2⟩t0 1  . Nα i=1 4t

To average over the initial time t 0, 155 partially overlapping MD segments have been extracted from every 70 ns production run. The starting times of consecutive segments are shifted by 200 ps from each other, and each segment covers 40 ns. Then, the mean square displacement of atoms and molecules is averaged over the 155 time segments before estimating the diffusion coefficient Dα from the linear interpolation of the 30 ≤ t ≤ 40 ns portion of the timedependent ∆2α (t). Representative results for Dα are reported in Table IV. Error bars have been estimated by computing Dα over four subsets of initial times chosen at random among the 155 considered for the global average in Eq. (11). Among the different species present in our samples, the largest relative uncertainty on Dα concerns POPC. The major difficulty is represented by the low value of DPOPC that is nearly three orders of magnitude lower than that of water. In the 40 ns considered for ∆2α (t), POPC molecules move, on average, by less than 10 Å that is larger but still somewhat comparable to the rattling motion of POPC molecules in the local cage of neighbouring molecules. Moreover, the mass, complexity, and self-organisation of POPC into bilayers are responsible for long relaxation times78,79 that are reflected in slight but visible deviations of ∆2α (t) from linearity up to at least 30 ns, as can be seen in Fig. 11. The choice of the interval for the linear fit of ∆2α (t), therefore, affects the estimate for DPOPC, with longer times corresponding to lower estimates of the POPC diffusion constant. The diffusion coefficient of POPC molecules in the pure sample #1 turns out to be DPOPC = 2.0 ± 0.2 × 10−8 cm2/s in the pure sample #1 and DPOPC = 1.80 ± 0.2 × 10−8 cm2/s

(11)

in the pure sample #2. These values are in fair and perhaps partly fortuitous agreement with the results of Ref. 27 (DPOPC = 1.2 ± 0.41 × 10−8 cm2/s) and somewhat disagree with the result of Ref. 12 (DPOPC = 3.9 ± 0.3 × 10−8 cm2/s), despite being computed with virtually the same model potential. The difference might be due partly to the different hydration level and partly to the different time span used for the linear interpolation of the time-dependent mean-square displacement.

FIG. 11. Mean square displacement of POPC molecules as a function of time in pure sample #1, and in two samples with 0.5M RTIL solutions: #3 ([bmim][Cl]) and #5 ([bmim][PF6]). The shaded area covers the time interval considered in the linear interpolation whose linear coefficient is used to compute the diffusion coefficient of POPC.


124706-14


The comparison of the computed DPOPC with experiments is difficult, since even experimental estimates, obtained for instance by quasi-elastic neutron scattering (QENS80), fluorescence recovery after photo-bleaching (FRAP),81 or nuclear magnetic resonance (NMR),82 can vary significantly because of the different length and time scales probed by these techniques.83 These differences between methods that, in most cases, probe time scales much longer than simulation suggest that finite time errors persist over fairly long times, and our results for the diffusion coefficient of lipids, estimated from portions of ∆2POPC(t), might still be affected by small but non-negligible errors due to the residual curvature of ∆2POPC(t) as a function of time. It was observed that fitting the 2 ≤ t ≤ 8 ns portion of ∆2POPC(t) provides results approaching those of FRAP,12 but this observation is very empirical, perhaps too much so to be used as a general basis for comparisons between computations and experiments. As expected, the mean square displacement of lipids as a function of time is highly anisotropic. Although growing slowly, ∆2POPC(t) does not saturate, confirming again the liquidlike character of the bilayer. Its 1D counterpart along z, instead, does saturate at a value of 18 ± 1 Å2, due primarily to slow, wide elongation fluctuations of the POPC tails. An even lower value of 11 ± 1 Å2 is obtained for the mean square displacement of the P atoms, representing the motion along z of POPC polar heads. Compared to lipids, the diffusion coefficient of water is much easier to compute, partly because of the large number of water molecules that eases the averaging implied by Eq. (11) and, more importantly, because D w is relatively high on the scale of our 70 ns simulations, thus reducing the weight of short-time oscillations. The overall 2D diffusion constant of water in sample #1 turns out to be 1.20 ± 0.05 × 10−5 cm2/s. A reference value to compare to is represented by the diffusion constant of bulk water that at T = 300 K is equal to 2.65 × 10−5 cm2/s both in experiments and in SPC/E simulations. A fair comparison, however, has to take into account the mean square displacement in two dimensions only, thus reducing the bulk value to 1.77 × 10−5 cm2/s. Therefore, the diffusion constant of water in the thin films in between POPC layers is roughly 30% less than the 1.77 × 10−5 cm2/s reference value for bulk water. The reduction of D w is not surprising, given the constrained geometry of water in our samples, and the general stabilisation of hydrogen bonds at the lipid/water interface that will be discussed below. Part of the D w reduction, in particular, is due to water molecules hydrogen bonded to POPC,84 whose mobility is much reduced with respect to that of water molecules floating in the middle of the inter-bilayer water film. Since the motion of water molecules is confined to a single inter-bilayer film, the z component of the mean square displacement saturated to a finite value also in the case of water, but the limiting value (∼120 Å2) is higher than that of POPC, since the motion of water molecules explores the entire width (∼2 nm) of the inter-bilayer solvent film. Because of their large mass, and multiple, distributed contacts with water and other lipids, the motion of POPC

J. Chem. Phys. 142, 124706 (2015)

is fairly continuous,85 while water moves in a rather tumble and jump way, especially in close proximity to the bilayer.86 To determine the variation of diffusion properties due to the addition of RTIL ions, we computed the diffusion coefficient of all species in sample #3 and sample #5. The results are collected in Table IV. The features that we highlight for further discussion can be summarised as follows: (i) addition of either [bmim][Cl] or [bmim][PF6] at 0.5M concentration further decreases the diffusion coefficient of water; (ii) the decrease of D w is relatively more important (−30%) in the [bmim][Cl] case than in [bmim][PF6] (−20%); (iii) the POPC diffusion coefficient DPOPC increases in the [bmim][Cl]-added sample and decreases in the [bmim][PF6] case; (iv) the same [bmim]+ ion diffuses faster when added as [bmim][Cl] than as [bmim][PF6]; and (v) [Cl]− diffuses faster − than [PF6] , even taking into account their different masses through a MPF6/MCl = 2.02 scaling factor.87 To distinguish the role of RTIL/water and RTIL/POPC interactions in changing the Dα values, we computed diffusion coefficients in homogeneous solutions of RTIL ions in water at 0.5M concentration, using samples of Nw = 24 040 water molecules, and 224 ion pairs. The results, collected in the second half of Table IV, allow us to rationalise most of the observations (i) to (v). For instance, the slower diffusion of water in sample #3 with respect to sample #5 reflects a similar ordering of D w in homogeneous solutions, i.e., without bilayers, and since the cation is the same, this can be safely attributed to water-anion interactions. Already in homogeneous solutions (no POPC), the same [bmim]+ ion diffuses faster when coupled to [Cl]− than to [PF6]−, possibly because of a slightly different dissociation degree of [bmim][Cl] (more dissociated) and of [bmim][PF6] (less dissociated) in water. The difference, however, is amplified in the presence of the POPC bilayers, pointing to the enhanced absorption of [bmim]+ into POPC in sample #5 as the most likely explanation. In a similar way, the presence of POPC bilayers greatly amplifies the mobility advantage of [Cl]− with respect to [PF6]−, consistently with the enhanced adsorption of [PF6]− at the POPC/water interface (see Fig. 9) with respect to [Cl]−. The only observation that is left unexplained by these comparisons is the opposite change of POPC mobility when [bmim][Cl] (increase) or [bmim][PF6]− (decrease) are added to the solution. To be precise, the variations in DPOPC due to RTIL ions in solution are larger than the statistical error bars, but they might still be affected by the systematic error due to long time relaxations, reflected in the slight curvature of ∆2POPC(t) on the 20 − 40 ns scale. However, our observation qualitatively agrees with the results of Ref. 27 and Ref. 25. We think, therefore, that the opposite effect of [bmim][Cl] and [bmim][PF6] on the lipid diffusion constant is a real effect, due to the different interaction of POPC with the RTIL ions, through a mechanism that we have not identified yet. The dynamics of hydrogen bonds has been analysed by first identifying all HB at an initial time t 0 and then monitoring their intact/broken state at regular time intervals. This analysis has been refined to distinguish the contribution from HB between two water molecules and those connecting water molecules to POPC oxygens. In both cases, occasional large


124706-15


J. Chem. Phys. 142, 124706 (2015)

TABLE V. Lifetime in ps of hydrogen bonds connecting water molecules (τ w ) and donated by water to oxygen 1 (τ 1), 2 (τ 2), 3 (τ 3), and 4 (τ 4), respectively. The numbering of oxygen atoms is specified in Fig. 1. Sample #1 #2 #3

τw

τ1

τ2

τ3

τ4

3.80 ± 0.04 4.54 ± 0.04 4.62 ± 0.04

126 ± 3 148 ± 3 142 ± 3

62 ± 3 72 ± 3 72 ± 3

59 ± 2 72 ± 2 71 ± 2

61 ± 2 71 ± 2 73 ± 2

geometrical (i.e., 2D) surface. We apply this concept to each of the POPC/water interfaces (4 in samples #1, #3, and #5 and 8 in sample #2, #4, and #6), identified with the method described in Sec. II. The instantaneous position of each of these surfaces is Fourier transformed, according to  z(x i , yi ) = Cqei[q x x i +q y yi ], (12) q

amplitude rotations (librations) of water molecules cause the intermittent split of bonds. To remove this flickering effect, we consider intact also the HBs that break but reform within a ps. Analysis of the HB persistence as a function of time displays the same qualitative features seen in previous studies,84,88 with apparent non-exponential stages at very short and long time scales. The results for the average lifetime of different classes of HB in our systems are collected in Table V. As expected, the lifetime τH B of HB donated by water to POPC oxygens is nearly two orders of magnitude longer than the lifetime of HB shared between water molecules. This observation might highlight the microscopic origin underlying the reduced overall mobility and molecular compressibility of water discussed in Subsections III A–III C. Considering here only the pure samples #1 and #2, we remark that our results fully agree with those of Ref. 84. More in detail, in pure samples without RTILs, we find that the τH B of bonds donated to the non-ester oxygens of POPC (3 and 4 in Fig. 1) is of the order of 60 ps. The τH B of bonds donated to the carbonyl oxygen 2 on the palmitoyl tail of POPC is nearly the same, as can be verified in Table V. Bonds donated to the carbonyl oxygen 1 on the oleoyl tail, instead, are significantly more stable, with a lifetime exceeding 120 ps. The non-equivalence of the two carbonyl oxygens is likely to be due again to their different equilibrium depths beneath the lipid/water interface. Addition of RTIL ions, in all cases, stabilises POPC-water hydrogen bonds and increases their τH B by about 20%. The detailed reason for this stabilisation is not clear, being consistent with observations such as the decrease of water diffusion coefficient and increase of the surface compressibility modulus upon RTIL addition, but at odds with other observations such as the general increase of compressibility and decrease of interfacial tension and especially bending rigidity. A broader exploration of lipid/water+RTIL systems might clarify these aspects. The HBs connecting water molecules have a short lifetime, of the order of 3.8 ps, that increases in the same proportion of the water-POPC case to 4.5 ps when RTIL ions are added to the solution. Without the sub-ps recombination provision, the lifetime of water-water hydrogen bonds is of the order of 1 ps,88 and the effect of RTIL ions cannot be detected. E. Bilayer dynamics

The vast difference in the lateral and perpendicular dimension of phospholipid bilayers justifies their representation as elastic membranes, whose shape is described by an ideal

where the q belong to the reciprocal lattice vector of the 2D n lattice q = (qx , q y ) ≡ 2π( Ln xx , L yy ). The {Cq} coefficients depend on time and can be used as generalised coordinates to describe the structure and dynamics of the interface. Such a description is particularly useful at long wavelengths and for single bilayers in excess water, since in that limit, the eigenstates of the interface Hamiltonian can be faithfully approximated by single plane waves. In such a case, general statistical mechanics arguments predict that at equilibrium ⟨|Cq(t)|2⟩ =

1 k BT , 2 A0 γq + k c q4

(13)

where ⟨. . .⟩ indicates time averaging, A0 is the cross sectional area covered by the bilayer, γ is the interfacial tension, and k c is the bending rigidity of the POPC/water interface. Equally general arguments89 state that for a liquid-like bilayer in a fluid environment, the surface tension γ has to vanish at equilibrium. Our simulations are carried out in the NPT ensemble at T = 300 K and P = 1 atm, and since the cell can expand or shrink in the (x y) plane of the bilayer, the sum of tensions from all water/lipid interfaces practically vanishes once averaged over time. In other terms, the interfacial tension defined from inter-molecular forces vanishes by construction in our simulations. However, the presence of multiple independent bilayers separated by a relatively thin water layer can spoil the simple relation in Eq. (13) connecting fluctuations with elastic parameters such as γ and k c , since in this case and even in the long wavelength limit, plane wave displacement patterns along x and y are no longer eigenstates of the membrane Hamiltonian. In our study, the {Cq} computed from configuration selected every 200 ps have been used to carry out the time average indicated in Eq. (13). The long wavelength portion of the logarithmic plot of ⟨|Cq|2⟩ as a function of |q| is shown in Fig. 12 for samples #1, #3, and #5, while the inset shows the results for sample #1 over the full range of |q| accessible to our simulations. The analytical expression of Eq. (13) is used to fit the |q| ≤ 0.25 Å−1 portion of ⟨|Cq|2⟩ data. A few points at very low q (q < 0.08 Å−1) are excluded from the fit because they are apparently affected by the finite size of the sample. Even after these exclusions, the relation between log⟨|Cq|2⟩ and the log q is not linear, unless the range of q is severely restricted. This observation, in turn, shows that both the γq2 and the k c q4 terms contribute significantly to the denominator of Eq. (13). This is confirmed by the quantitative estimate for the water/phospholipid interfacial tension γ and bending rigidity k c in pure samples and in the presence of RTIL ions that are reported in Table VI. We emphasise that this is the estimate of γ and k c based on fluctuations and differs from


124706-16


J. Chem. Phys. 142, 124706 (2015)

FIG. 12. Logarithmic plot of the square amplitude of fluctuation as a function of wave vector |q|. Symbols report simulation data, lines represent their interpolation by Eq. (13). Black empty circles and lines: pure sample #1; red full circles and line: sample #3 with [bmim][Cl]; blue empty squares and line: sample #5 with [bmim][PF6].

the estimate obtained from the mechanical route, based on the evaluation of forces and distances. From this reason, we label the interfacial tension and bending rigidity obtained from the fit of ⟨|Cq|2⟩ with a capital F (γ F , k cF ), to distinguish them from their mechanical counterpart. In all samples, the water/phospholipid interfacial tension (γ F ∼ 15 mJ/m) estimated from fluctuations is about five times lower than at the liquid water surface (72 mJ/m) and also lower than the surface tension of the oil/air surface, estimated at ∼30 mJ/m90 but still significantly above zero, and certainly not as low as in previous simulations.25 As already stated, we think that this deviation of elastic properties computed from the fluctuation and mechanical route is due primarily to the relatively low hydration level that we selected for our simulations, corresponding to a 2 nm inter-layer of water. At this low thickness, the mechanical properties of water are apparently very different from those of the bulk fluid

TABLE VI. Interfacial tension γ and bending rigidity K c of the water/lipid interface. Sample #1 #2 #3 #4 #5 #6

γ (mJ/m) 14.1 ± 0.1 13.8 ± 0.1 11.4 ± 0.15 10.9 ± 0.2 12.0 ± 0.1 11.7 ± 0.2

1020 × K

c

(J)

1.2 ± 0.1 1.2 ± 0.1 0.65 ± 0.1 0.6 ± 0.1 0.6 ± 0.1 0.5 ± 0.1

phase, and, perhaps more importantly, the proximity of other bilayers reduces the amplitude of the interface fluctuations.91 The bending rigidity k cF is comparable to the one computed by simulation in similar systems, and together with the relatively low value of γ F , is responsible for the apparent curvature in the log⟨|Cq|2⟩ function even at moderate |q| values. The same analysis of long range fluctuations discussed in the previous paragraphs for individual lipid/water interfaces can be applied equally well to whole bilayers, defining their instantaneous configuration as the average of their two limiting lipid/water geometric surfaces. It is easy to verify that under mild conditions, the two descriptions are equivalent and give the same estimate of γ F and k cF , provided the reference area A0 in Eq. (13) is replaced by 2A0, accounting for the two interfaces created by each bilayer. Adding RTILs has the effect of decreasing somewhat the surface tension γ F and of halving the bending rigidity k cF . This is consistent with the global picture of reduced stability for the lipid bilayer, and the stronger effect on k cF might be related to the localised disruption of the bilayer at the insertion point of the [bmim]+ cations. The deviation of γ F and k cF from their mechanical counterparts has been attributed to the presence in our simulations of multiple independent bilayers and to the correlation in their motion. This effect can be quantified by looking at the time evolution of the {Cq} coefficients computed at low q for neighbouring surfaces, separated by a lipid bilayer or by a thin water film. To this aim, we write each Cq coefficient as the product of a real amplitude times a complex phase factor, Cq(t) = Aq(t)eiϕ(t).

(14)

The plot of these two geometric parameters as a function of time for adjacent surfaces reveals strong correlations at low |q| (see Fig. 13) and moderate but still apparent correlations at higher |q|. Correlations are stronger for surfaces separated by a water film than by the lipid layer because of the wider thickness and lower bulk modulus of this last with respect to the thinner and stiffer water interlayer. The four bilayer samples allow us to verify that, as expected, correlations among Cq decay quickly with increasing separation between the reference surfaces. Our observation of strong correlations among fluctuations at different interfaces disagrees with the conclusions of Ref. 71(c). The disagreement might be due to the different hydration level of their samples and to the longer wavelength fluctuations analysed in our samples.

IV. SUMMARY AND CONCLUSIVE REMARKS

Molecular dynamics simulations have been carried out to investigate the effect of RTIL solvation on the structure and dynamics of planar stacks of POPC lipid bilayers and of the thin water layer hydrating and stabilising them. The RTIL species considered in our computations consist of [bmim][Cl] and [bmim][PF6], whose interaction with POPC is driven by the incorporation of [bmim]+ into the bilayer. Simulations are carried out by MD in the NPT ensemble, using samples of nlip = 1360 POPC molecules, with a water


124706-17


J. Chem. Phys. 142, 124706 (2015)

FIG. 13. Time dependence of the amplitude ( Aq(t)) and phase (φ q(t)) of the Fourier coefficients {C q}(t) (see Eqs. (12) and (14)) for adjacent interfaces (red dots and blue empty circles, respectively) separated by a lipid or by a water film. Results for sample #1, at q = (2π/L x, 0), i.e., for the longest wavelength surface modulation in our samples. Dashed lines are a guide to the eye.

content 17 ≤ n w ≤ 20 molecules per POPC, corresponding to a basic bilayer-bilayer separation of 6.0 nm, approaching the 6.4 nm measured by experiments on POPC multi-lamellar vesicles at full hydration.6 Two thirds of the sample volume are occupied by POPC, each lipid bilayer is ∼4 nm thick, and thus, each inter-bilayer water film is ∼2 nm thick. RTIL ions are added to the water film at concentration 0.5M. The relatively large number (224) of POPC molecules allows us to consider from two to four distinct bilayers in each sample, opening the way to the investigation of bilayer-bilayer interactions and correlation, both in pure and RTIL-enriched systems. On average, the POPC bilayers and the water films are homogeneous along their nominal plane, and thus, position dependent properties such as the density, concentration, or electrostatic potential profiles are computed and displayed as a function of the coordinate z perpendicular to the bilayer plane. First of all, visualisation of simulation snapshots confirms the qualitative picture proposed in previous studies of the same or similar systems,24–27 in which [bmim]+ cations enter the POPC bilayer tail-first, solvating their butyl chain into the POPC carbohydrate tails. The ionic imidazolium ring of the absorbed cations tends to reside in close proximity of the most polar regions in the phospholipid molecule. These are represented by the two carbonyl groups at the junction with the palmitoyl/oleoyl tails and by the two non-ester oxygens at the phosphate group. On the simulation time scale, the

incorporation of [bmim]+ into POPC is virtually irreversible. The [Cl]− anions remain dispersed in water, while a sizeable portion of the [PF6]− ions is adsorbed at the lipid/water interface. The difference in the behaviour of the two anions is caused primarily by their different radii, with the larger size of [PF6]− decreasing the strength of its water solvation. An additional important role is played also by dispersion forces that are larger for [PF6]− and enhance its affinity for amphiphilic molecules such as POPC. An equivalent description of these same features could be formulated in terms of the different hydrophobicity/hydrophilicity of the [Cl]− and [PF6]− anions. The POPC bilayer retains its basic integrity throughout the RTIL absorption. The compressibility of the whole lipid and water composite samples increases slightly but systematically upon addition of both [bmim][Cl] and [bmim][PF6], pointing to a decrease of the cohesion and stability of the ordered bilayer stack. At the same time, however, we observe the increase of the surface compressibility modulus K A and a smaller but still significant increase of the Young moduli Eα , α = x, y, and z (again of the whole sample), which suggests instead a general stiffening and increase of cohesion following the RTIL addition. As briefly discussed in Subsections III A–III C, the contradictory trends in these diagnostic properties reflect important correlations in the fluctuations of different geometrical parameters that we still do not fully understand. The bilayer cross section per POPC molecule increases slightly upon adding the RTIL. To retain the same POPC


124706-18


molecular volume, the bilayer thickness has to decrease by ∼0.2 Å. This somewhat counterintuitive thinning upon absorption of [bmim]+ into the lipid phase agrees qualitatively with the results of neutron reflectometry measurements22 that reveal an even larger shrinking of up to 1.0 Å. Our observations concerning the surface area and bilayer thickness are also in qualitative agreement with previous experimental determinations of geometrical parameters in similar bilayers in which neutral organic molecules (curcumin, magainin) are absorbed.55,92 The analysis of volumes and of their fluctuations can be specialised to consider the individual compressibilities of the water, POPC, and RTIL components. To this aim, we implemented a Monte Carlo estimate of the volume of Voronoi atomic and molecular cells, whose fluctuations provide the partial compressibilities of interest for our analysis.66 We find that the low compressibility of each sample is due primarily to the stiffness of the inter-bilayer water films, while the compressibility of lipids is comparable to that of similar organic systems at ambient conditions. Addition of RTIL ions increases significantly the compressibility of both water and POPC, but enhanced correlations in their volume fluctuations somewhat reduce the impact on the global sample compressibility to a few percent. An interesting effect revealed by our computations is a slight but nevertheless significant non-additivity in the lipid, water, and RTIL volumes. Although quantitatively small, a negative volume additivity is an indicator of cohesion. We observe that the addition of RTIL ions to the solution decreases the size of this effect, once again suggesting a decreased stability of the ordered lipid stack. In pure samples, the number and charge density profiles of lipids and water across the bilayers agree with the results of previous computations. The same observation applies to the electrostatic potential profile. The variation of these functions upon adding RTIL provides an average but quantitative view of the RTIL interaction with POPC, highlighting once more the [bmim]+ absorption in the lipid phase, and illustrates once again the role of the anion: [Cl]− shows a relatively weak tendency to absorption, while [PF6]− is more likely to be adsorbed at the lipid/water interface. RTIL ions introduce a sizeable amount of free charge into the system, but changes in the orientation of lipids and water molecules effectively screen this perturbation, and the electrostatic signature of each bilayer changes only moderately following the RTIL addition. The diffusion constant of POPC along the (x, y) plane is in the range of 10−8 cm2/s and thus is nearly three orders of magnitude slower than that of water. For this reason, its estimation by molecular dynamics simulations is still challenging. Nevertheless, our computations show that POPC molecules diffuse over the 100 ns time scale, confirming the disordered liquid (L d ) state of the bilayer. Moreover, we found that the absorption of the same [bmim]+ ion has the opposite effect on the POPC mobility in the case of [bmim][Cl] and of [bmim][PF6], apparently for the different coupling of the cation and POPC motion with the [Cl]− or [PF6]− anion dynamics. On the other hand, the diffusion constant of [bmim]+ is faster in the [bmim][Cl] than in the [bmim][PF6] case, probably reflecting the different partition of cations between the lipid and water phases. We also observe that the

J. Chem. Phys. 142, 124706 (2015)

diffusion constant of [Cl]− is faster than that of [PF6]−, even accounting for their mass difference. We think that the slow dynamics of [PF6]− might be due to its closer association with the POPC bilayer. We identify the position of the lipid/water dividing surface implementing a majority rule similar in spirit to previous computational approaches.42,43 This construction allows us to analyse the fluctuations in the surface shape and position, from which we extract an estimate of the interfacial tension γ F and bending rigidity k cF based on equilibrium fluctuations. As expected, the lipid/water interfacial tension of pure samples is much less than that of comparable interfaces, such as the liquid/vapour surface of water and of electrolyte solutions, but still significantly different from zero. This last observation is somewhat surprising, since the NPT simulation automatically ensures zero mechanical tension along the three Cartesian directions. The tension estimated from fluctuations, however, deviates from the mechanical tension, apparently because the low hydration of our samples introduces strong correlations in the motion of adjacent bilayers that remain important even at long wavelengths. Addition of [bmim][Cl] or [bmim][PF6] increases the amplitude of the surface fluctuation and thus reduces the interfacial tension and, even more, the bending rigidity. Both observations point again to a decreased stability of the bilayer phase. The definition of the geometric surface dividing POPC from water allows us also to define an intrinsic profile of all species, in which the origin of the z axis is a function of the x and y coordinates along the plane. The aim of such a construction is to remove the effect of long wavelength, large amplitude capillary oscillations of the interface that could conceal important details in the composition profiles. The resulting intrinsic profiles, indeed, display a richer structure than their Cartesian counterparts and provide a more detailed view of the cation and anion distributions at the interface. Moreover, intrinsic profiles show that the [bmim]+ absorption is accompanied by a small but apparent penetration of water into the head region of POPC, an observation that is confirmed by visual analysis of snapshots. The identification of the lipid/water dividing surface allows us to conclude that such an interface is rough, with roughness defined in the statistical mechanics sense and identified by a logarithmic divergence of the height-height correlation function along the interface. Addition of RTIL increases the interface delocalisation, and thus the roughness. This result could be tested by an (challenging) analysis of neutron reflectometry data. Hydrogen bonding is an essential component of all phospholipid/water systems. In pure samples, the properties of the HB network computed in our simulations agree again with the results of previous studies: water donates its protons preferentially to the non-ester oxygens (3 and 4, see Fig. 1) of POPC. The highly polar carbonyl oxygens (1 and especially 2) are somewhat screened by the POPC head from a more intimate contact with water. Inclusion of RTILs increases slightly the POPC-water HB, presumably because of the enhanced water penetration into the bilayer driven by the [bmim]+ absorption. The lifetime τH B of HB donated by water to POPC oxygens at


124706-19


positions 1, . . . , 4 provide detailed information of the dynamics of POPC and especially of water at the lipid/water interface. Addition of RTIL ions to the solution increases all lifetimes by ∼20%, providing a diagnostic property that could be probed by inelastic neutron scattering. Our extensive simulations and analysis extend our knowledge of the partitioning of RTIL species between the water and lipid phases and of its effect on the bilayer properties. The subject, however, is far from exhausted, since obvious and important questions remain open. First of all, we still have little control of the thermodynamics of [bmim]+ absorption and anion adsorption at the lipid/water interface, when present. Moreover, extending the simulation size and time scales by one order of magnitude (up to 50 × 50 × 50 nm3 and up to the µs range) is both feasible and worthwhile, since it could reveal important changes in the lipid stability and structure caused by the RTIL solvation. Such an extension could be achieved, for instance, by coarse grained approaches and could reveal the structure and in-plane distribution of molecularly thin pores puncturing the phospholipid bilayer. The formation of such pores upon the addition of [bmim][Cl](at high concentration) was inferred in Ref. 19 from the leakage of luminescent molecules trapped inside phospholipid vesicles. Pores, however, have not been seen in the present or previous simulations, probably because of size and especially time limitations. Needless to say, only a few phospholipids and a few imidazolium-based compounds have been considered so far in computational and experimental studies of lipid-RTIL interactions,25–28 while both the lipid and the RTIL classes are large and diverse. Last but not least, all simulation studies of these systems rely on force field models that are known to be rather primitive approximations to their potential energy surfaces. Ab-initio methods are far from being suitable to replace empirical force fields in full blown simulations, but improvements in the force field approach to include polarisation and many-body interactions are possible and desirable. From the application point of view, we speculate that the affinity of the positively charged imidazolium ring for the polar head of phospholipids could provide a novel hook for tethering bilayers to solid surfaces,93 upon tuning the choice of the functional groups at the 1 and 3 positions on imidazolium. The group at position 1 could remain the same alkane chain of standard [Cn mim]+ compounds, whose affinity for the neutral lipid tails provides the driving force for the tethering action. The methyl group at position 3, however, should be replaced by a moderately polar and hydrophilic polymer chain, easily solvated by the water environment in which the bilayer resides.

J. Chem. Phys. 142, 124706 (2015)

APPENDIX: STATISTICAL VORONOI VOLUMES

The Voronoi construction represents a suitable procedure to attribute a volume to atoms and thus also to molecules in extended and relatively homogeneous systems. Unfortunately, the construction is cumbersome in large, complex samples such as those considered in our case. For this reason, we adopt a simpler approach, giving the volume of the Voronoi cell for each atom, without explicitly computing its sides, edges, and vertices. Since, by definition, the Voronoi cell of each atom i is the geometrical locus of points closest to i than to any other atom, we uniformly distribute M points over the simulation cell, and we count the number ni of points having atom i as their nearest neighbour. The Voronoi volume attributed to each atom is estimated as vi = V ni /M, where V is the instantaneous system volume. The accurate determination of vi for each atom requires M to be at least 5000 times larger than the number of atoms, and possibly more, because of the characteristically slow √ 1/ ni decay of the relative error in purely random processes. However, the identification of the closest atom can be made very efficient by using a neighbouring list or subdividing the orthorhombic sample into partially overlapping sub-cells. Moreover, the estimation of the average volume of sizeable subsets of atoms, such as lipids, water, and ionic liquids, is far less expensive (see Fig. 14), since in this case, the error is averaged out over thousands of atoms. In our computations, convergence is further enhanced by averaging over different snapshots from each simulation. Averaging, however, can be used only less extensively in quantifying the mean square fluctuations of atomic and molecular volumes, since low accuracy data would add spurious fluctuations from the sampling, resulting into an overestimation of both fluctuations and partial compressibilities. For this reason, the computation of atomic and molecular compressibility is somewhat more expensive than the evaluation of average volumes but remains a task suitable for a laptop. The volume and compressibility data discussed in the paper, for instance, have been computed upon distributing

ACKNOWLEDGMENTS

We thank G. Cottone for a careful reading of the manuscript. A.B. acknowledges support from the European Community under the Marie Curie Intra-European Fellowship for Career Development (IEF) Grant HYDRA No. 301463; computations have been carried out on the Stokes and Fionn supercomputer of ICHEC, whose access has been made possible by a Class B Grant No. ndphy051b. P.B. acknowledges support from the Istituto Italiano di Tecnologia (IIT) under the SEED project SIMBEDD Grant No. 259.

FIG. 14. Convergence of the estimated volume per atom of POPC and of water as a function of the number of sampling points per atom (total: lipids + water) used in the statistical Voronoi construction described in the Appendix.


124706-20


107 points on each simulation frame, requiring a few minutes per frame on a single core. For consistency with the united atom representation of POPC implied by our force field, in the analysis of Voronoi volumes, we consider non-hydrogen atoms only also in the case of water and of RTIL ions. 1P.

Ball, Chem. Rev. 108, 74 (2008). M. Berg, J. L. Tymoczko, and L. Stryer, Biochemistry, 7th ed. (W.H. Freeman, New York, 2012). 3H. Khesbak, O. Savchuk, S. Tsushima, and K. Fahmy, J. Am. Chem. Soc. 133, 5834 (2011). 4T. Welton, Chem. Rev. 99, 2071 (1999). 5N. H. Tattrie, J. R. Bennett, and R. Cyr, Can. J. Biochem. 46, 819 (1968). 6N. Kuˇ cerka, S. Tristram-Nagle, and J. F. Nagle, J. Membr. Biol. 208, 193 (2005). 7J. F. Nagle and S. Tristram-Nagle, Biochim. Biophys. Acta 1469, 159 (2000). 8S. Leekumjorn and A. K. Sum, J. Phys. Chem. B 111, 6026 (2007); M. Ceccarelli and M. Marchi, Biochimie 80, 415 (1998); A. W. Mauk, E. L. Chaikof, and P. J. Ludovice, Langmuir 14, 5255 (1998); H. Heller, M. Schaefer, and K. Schulten, J. Phys. Chem. 97, 8343 (1993). 9T. D. Lazzara, A. Janshoff, and C. Steinem, in Handbook of Nanomaterials Properties, edited by B. Bhushan, D. Luo, S. R. Schricker, W. Sigmund, and S. Zauscher (Springer, Berlin, 2014). 10M. L. Berkowitz and R. Vácha, Acc. Chem. Res. 45, 74 (2012). 11R. Zimmermann, D. Küttner, L. Renner, M. Kaufmann, and C. Werner, J. Phys. Chem. A 116, 6519 (2012). 12R. A. Böckmann, A. Hac, T. Heimburg, and H. Grubmüller, Biophys. J. 85, 1647 (2003). 13S. A. Tatulian, in Phospholipids Handbook, edited by G. Cevc (Marcel Dekker, New York, 1993), pp. 511-552; P. T. Vernier, M. J. Ziegler, and R. Dimova, Langmuir 25, 1020 (2009). 14Ionic Liquids: From Knowledge to Applications, ACS Symposium Series Vol. 1030, edited by N. V. Plechkova, R. D. Rogers, and K. R. Seddon (American Chemical Society, Washington, DC, 2009). 15M. Petkovic, K. R. Seddon, L. P. N. Rebelo, and C. S. Pereira, Chem. Soc. Rev. 40, 1383 (2011); J. Ranke, A. Müller, U. Bottin-Weber, F. Stock, S. Stolte, J. Arning, R. Störmann, and B. Jastorff, Ecotoxicol. Environ. Saf. 67, 430 (2007); D. J. Couling, R. J. Bernot, K. M. Docherty, J. K. Dixon, and E. J. Maginn, Green Chem. 8, 82 (2006). 16B. K. Paul, A. Ganguly, and N. Guchhait, J. Photochem. Photobiol., B 133, 99 (2014); V. F. Curto, S. Scheuermann, R. M. Owens, R. Vijayaraghavan, D. R. MacFarlane, F. Benito-Lopez, and D. Diamond, Phys. Chem. Chem. Phys. 16, 1841 (2014); N. Byrne and C. A. Angell, Chem. Commun. 2009, 1046. 17G. Portella, M. W. Germann, N. V. Hud, and M. Orozco, J. Am. Chem. Soc. 136, 3075 (2014); M. Nakano, H. Tateishi-Karimata, S. Tanaka, and N. Sugimoto, J. Phys. Chem. B 118, 379 (2014); A. Chandran, D. Ghoshdastidar, and S. Senapati, J. Am. Chem. Soc. 134, 20330 (2012); R. Vijayaraghavan, A. Izgorodin, V. Ganesh, M. Surianarayanan, and D. R. MacFarlane, Angew. Chem., Int. Ed. 49, 1631 (2010). 18S. Jeong, S. H. Ha, S.-H. Han, M.-C. Lim, S. M. Kim, Y.-R. Kim, Y.-M. Koo, J. S. So, and T.-J. Jeon, Soft Matter 8, 5501 (2012). 19K. O. Evans, Colloids Surf., A 274, 11 (2006). 20K. O. Evans, J. Phys. Chem. B 112, 8558 (2008). 21K. O. Evans, Int. J. Mol. Sci. 9, 498 (2008). 22A. Benedetto, F. Heinrich, M. A. Gonzalez, G. Fragneto, E. Watkins, and P. Ballone, J. Phys. Chem. B 118, 12192 (2014). 23M. Galluzzi, S. Zhang, S. Mohamadi, A. Vakurov, A. Podestà, and A. Nelson, Langmuir 29, 6573 (2013). 24S. R. T. Cromie, M. G. Del Pópolo, and P. Ballone, J. Phys. Chem. 113, 11642 (2009). 25R. J. Bingham and P. Ballone, J. Phys. Chem. B 116, 11205 (2012). 26G. S. Lim, J. Zidar, D. W. Cheong, S. Jaenicke, and M. Klähn, J. Phys. Chem. B 118, 10444 (2014). 27B. Yoo, J. K. Shah, Y. Zhu, and E. J. Maginn, Soft Matter 10, 8641 (2014). 28H. Lee and T.-J. Jeon, Phys. Chem. Chem. Phys. 17, 5725 (2015). 29M. Klähn and M. Zacharias, Phys. Chem. Chem. Phys. 15, 14427 (2013). 30J. Katsaras, Biophys. J. 75, 2157 (1998). 31G. Fragneto, Eur. Phys. J.: Spec. Top. 213, 327 (2012). 32H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987). 33C. Oostenbrink, A. Villa, A. E. Mark, and W. F. van Gunsteren, J. Comput. Chem. 25, 1656 (2004). 2J.

J. Chem. Phys. 142, 124706 (2015) 34O.

Berger, O. Edholm, and F. Jähnig, Biophys. J. 72, 2002 (1997). Chiu, M. Clark, V. Balaji, S. Subramanian, H. J. Scott, and E. Jacobsson, Biophys. J. 69, 1230 (1995). 36J. N. Canongia Lopes, J. Deschamps, and A. A. H. Padua, J. Phys. Chem. B 108, 2038 (2004); J. N. Canongia Lopes and A. A. H. Padua, J. Phys. Chem. B 108, 16893 (2004); J. N. Canongia Lopes, J. Deschamps, and A. A. H. Padua, J. Phys. Chem. B 108, 11250 (2004). 37B. L. Bhargava and S. Balasubramanian, J. Chem. Phys. 127, 114510 (2007); T. Cremer, C. Kolbeck, K. R. J. Lovelock, N. Paape, R. Woelfel, P. S. Schulz, P. Wasserscheid, H. Weber, J. Thar, B. Kirchner, F. Maier, and H. Steinrueck, Chem. - Eur. J. 16, 9018 (2010). 38J. P. M. Jämbeck and A. P. Lyubartsev, J. Chem. Theory Comput. 8, 2938 (2012); N. Bhatnagar, G. Kamath, and J. J. Potoff, J. Phys. Chem. B 117, 9910 (2013); A short discussion of the state of the art in lipid modelling is in: M. Manna, T. Róg, and I. Vattulainen, Biochim. Biophys. Acta 1841, 1130 (2014). 39I. V. Voroshylova and V. V. Chaban, J. Phys. Chem. B 118, 10716 (2014); O. Borodin, ibid. 113, 11463 (2009). 40H. Berendsen, D. van der Spoel, and R. van Drunen, Comput. Phys. Commun. 91, 43 (1995); E. Lindahl, B. Hess, and D. van der Spoel, J. Mol. Mod. 7, 306 (2001). 41B. Hess, H. Bekker, H. J. C. Berendsen, and J. G. E. M. Fraaije, J. Comput. Chem. 18, 1463 (1997). 42E. Chácon and P. Tarazona, Phys. Rev. Lett. 91, 166103 (2003); F. Bresme, E. Chácon, P. Tarazona, and K. Tay, ibid. 101, 056102 (2008). 43W. J. Allen, J. A. Lemkul, and D. R. Bevan, J. Comput. Chem. 30, 1952-1958 (2009). 44A. Tardieu, V. Luzzati, and F. C. Reman, J. Mol. Biol. 75, 711 (1973). 45P. D. Tieleman, J. L. MacCallum, W. L. Ash, C. Kandt, Z. Xu, and L. Monticelli, J. Phys.: Condens. Matter 18, S1221 (2006). See also the web site: http://www.ucalgary.ca/tieleman/. The original configuration corresponded to POPE. It was changed to POPC, and the bilayer in water has been re-equilibrated with the POPC potential. 46D. M. Leneveu, R. P. Rand, and V. A. Parsegian, Nature 259, 601 (1976); T. J. McIntosh and S. A Simon, Annu. Rev. Biophys. Biomol. Struct. 23, 27 (1994). 47S. J. Marrink, M. Berkowitz, and H. J. C. Berendsen, Langmuir 9, 3122 (1993). 48K. Gawrisch, H. C. Gaede, M. Mihailescu, and S. H. White, Eur. Biophys. J. 36, 281 (2007). 49R. N. A. H. Lewis, B. D. Sykes, and R. N. McElhaney, Biochemistry 27, 880 (1988). 50K. R. Seddon, A. Stark, and M.-J. Torres, Pure Appl. Chem. 72, 2275 (2000). 51J. L. Antony, E. J. Maginn, and J. F. Brennecke, J. Phys. Chem. B 105, 10942 (2001). 52P. A. Hyslop, M. Morel, and R. D. Sauerheber, Biochemistry 29, 1025 (1990). 53G. Pabst, M. Rappolt, H. Amenitsch, and P. Laggner, Phys. Rev. E 62, 4000 (2000). 54J. F. Nagle and S. Tristram-Nagle, Curr. Opin. Struct. Biol. 10, 474 (2000). 55S. Ludtke, K. He, and H. Huang, Biochemistry 34, 16764 (1995). 56R. Kubo, Rep. Prog. Phys. 29, 255 (1966). 57W. Rawicz, K. C. Olbrich, T. McIntosh, D. Needham, and E. Evans, Biophys. J. 79, 328 (2000). 58Q. Waheed and O. Edholm, Biophys. J. 97, 2754 (2009). 59L. D. Landau and E. M. Lifshitz, Statistical Physics (Addison-Wesley, Reading, MA, 1958). 60N. M. Philip, Proc. Indian Acad. Sci. - Sect. A 9, 109 (1939). 61H. L. Pi, J. L. Aragones, C. Vega, E. G. Noya, J. L. F. Abascal, M. A. Gonzalez, and C. McBride, Mol. Phys. 107, 365 (2009). 62Z.-Y. Gu and J. F. Brennecke, J. Chem. Eng. Data 47, 339 (2002). 63W. Schrader, H. Ebel, P. Grabitz, E. Hanke, T. Heimburg, M. Hoeckel, M. Kahle, F. Wente, and U. Kaatze, J. Phys. Chem. B 106, 6581 (2002). 64CRC Handbook of Chemistry and Physics, 95th ed., edited by W. M. Haynes (Francis and Taylor, Boca Raton, 2014). 65P. W. Bridgman, Proc. Am. Acad. Arts Sci. 47, 347 (1911); O. V. Nagornov and V. E. Chizhov, Zhur. Prik. Mek. Tekh. Fiz. 3, 41 (1990). 66M. Marchi, J. Phys. Chem. B 107, 6598 (2003). 67We considered only non-hydrogen atoms of water and RTIL molecules for consistency with the united atom representation of POPC provided by the Gromos force field. 68C. A. Croxton, Statistical Mechanics of the Liquid Surface (John Wiley & Sons, New York, 1980). 35S.-W.


124706-21 69R.


Vácha, P. Jurkiewicz, M. Petrov, M. L. Berkowitz, R. A. Böckmann, J. Barucha-Kraszewska, M. Hof, and P. Jungwirth, J. Phys. Chem. B 114, 9504 (2010). 70T. J. McIntosh, Curr. Opin. Struct. Biol. 10, 481 (2000). 71(a) S. W. I. Siu, R. Vácha, P. Jungwirth, and R. A. Böckmann, J. Chem. Phys. 128, 125103 (2008); (b) W. Shinoda, M. Shimizu, and S. Okazaki, J. Phys. Chem. B 102, 6647 (1998); (c) F. Zhou and K. Schulten, J. Phys. Chem 99, 2194 (1995); (d) M. A. Wilson and A. Pohorille, J. Am. Chem. Soc. 116, 1490 (1994). 72L. S. Vermeer, B. L. de Groot, V. Réat, A. Milon, and J. Czaplicki, Eur. Biophys. J. 36, 919 (2007). 73 n = 128, n = 64 in Samples #1, #3 and #5; n = 64, n = 64 in x y x y Samples #2, #4 and #6. 74A. Zangwill, Physics at Surfaces (Cambridge University Press, Cambridge, 1988). 75R. Kumar, J. R. Schmidt, and J. L. Skinner, J. Chem. Phys. 126, 204107 (2007). 76D. Xenides, B. R. Randolf, and B. M. Rode, J. Mol. Liq. 123, 61 (2006). 77J. S. Wilkes and M. J. Zaworotko, J. Chem. Soc., Chem. Commun. 1992, 965; A. K. Abdul-Sada, A. M. Greenway, P. B. Hitchcock, T. J. Mohammed, K. R. Seddon, and J. A. Zora, ibid. 1986, 1753; C. J. Dymek, Jr., D. A. Grossie, A. V. Fratini, and W. W. Adams, J. Mol. Struct. 213, 25 (1989). 78J. Yang, C. Calero, and J. Marti, J. Chem. Phys. 140, 104901 (2014). 79E. Flenner, J. Das, M. Rheinstädter, and I. Kosztin, Phys. Rev. E 79, 011907 (2009). 80W. Pfeiffer, G. Schlossbauer, W. Knoll, B. Farago, A. Steyer, and E. Sackmann, J. Phys. (France) 49, 1077 (1988); J. Tabony and B. Perly, Biochim. Biophys. Acta 1063, 67 (1990).

J. Chem. Phys. 142, 124706 (2015) 81W.

L. C. Vaz, R. M. Clegg, and D. Hallmann, Biochemistry 24, 781 (1985). Wennerstrom and G. Lindblom, Q. Rev. Biophys. 10, 67 (1977); A. L. Kuo and C. G. Wade, Biochemistry 18, 2300 (1979). 83W. L. C. Vaz and P. F. Almeida, Biophys. J. 60, 1553 (1991). 84S. Y. Bhide and M. L. Berkowitz, J. Chem. Phys. 123, 224702 (2005). 85E. Falk, T. Róg, M. Karttunen, and I. Vattulainen, J. Am. Chem. Soc. 130, 44 (2008). 86Y. von Hansen, S. Gekle, and R. R. Netz, Phys. Rev. Lett. 111, 118103 (2013). 87J. P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic Press, London, 1986). 88A. Luzar and D. Chandler, Phys. Rev. Lett. 76, 928 (1996); A. Luzar and D. Chandler, Nature 379, 55 (1996). 89F. Brochard, P. G. de Gennes, and P. Pfeuty, J. Phys. 37, 1099 (1976). 90F. J. Nagle, Ann. Rev. Phys. Chem. 31, 157 (1980); F. Jahnig, Biophys. J. 46, 687 (1984). 91F. Jähnig, Biophys. J. 71, 1348 (1996). 92W.-C. Hung, F.-Y. Chen, C.-C. Lee, Y. Sun, M.-T. Lee, and H. W. Huang, Biophys. J. 94, 4331 (2008). 93R. Neumann, S. M. Schiller, F. Gless, B. Grohe, K. B. Hartman, I. Kärcher, I. Köper, J. Lübben, K. Vasilev, and W. Knoll, Langmuir 19, 5435 (2010). 94See supplementary material at http://dx.doi.org/10.1063/1.4915918 for Table of the Young’s moduli E α , α = x, y and z of the simulated samples; dipole orientation of water molecules along the direction perpendicular to the POPC bilayers (FigS1.eps); charge density and electrostatic potential profiles for sample #1 (FigS2.eps); intrinsic profiles for sample #3 (FigS3.eps); acyl-chain order parameters SC D of all simulated systems (FigS4.eps); snapshots of typical hydration shells around closely associated POPC/[bmim]+ complexes (hydration.xyz file). 82H.


Structure and dynamics of room temperature ionic liquids with bromide anion: results from 81Br NMR spectroscopy.

Dissolution of cellulose in room temperature ionic liquids: anion dependence.

Ion transport and structural dynamics in homologous ammonium and phosphonium-based room temperature ionic liquids.

water solution: a neutron reflectometry study.

Dynamics of electrical double layer formation in room-temperature ionic liquids under constant-current charging conditions.

Ultrafast solvation dynamics and charge transfer reactions in room temperature ionic liquids.

Phosphonium chloromercurate room temperature ionic liquids of variable composition.

The joint effects of room temperature ionic liquids and ordered media on fluorescence characteristics of estrogens in water and methanol.

Orientational dynamics in a lyotropic room temperature ionic liquid.

Standard and absolute pKa scales of substituted benzoic acids in room temperature ionic liquids.

Ionic liquids at nonane-water interfaces: molecular dynamics studies.

Application of room-temperature aprotic and protic ionic liquids for oxidative folding of cysteine-rich peptides.

Water Solution: A Neutron Reflectometry Study".

Terahertz and Infrared Spectroscopy of Room-Temperature Imidazolium-Based Ionic Liquids.

Self-assembly of imidazolium-based surfactants in magnetic room-temperature ionic liquids: binary mixtures.

What type of nanoscopic environment does a cationic fluorophore experience in room temperature ionic liquids?

Dynamic propensity as an indicator of heterogeneity in room-temperature ionic liquids.

Properties of water confined in ionic liquids.

Stability of the Liquid State of Imidazolium-Based Ionic Liquids under High Pressure at Room Temperature.

Room-temperature ionic liquids meet bio-membranes: the state-of-the-art.

Physicochemical properties of tri-n-butylalkylphosphonium cation-based room-temperature ionic liquids.

Transparent and colourless room temperature ionic liquids having high refractive index over 1.60.

Room-temperature ionic liquids enhanced green synthesis of β-glycosyl 1-ester.

SEM observation of hydrous superabsorbent polymer pretreated with room-temperature ionic liquids.