Article pubs.acs.org/JPCB
Driving Force for Crystallization of Anionic Lipid Membranes Revealed by Atomistic Simulations Bao Fu Qiao† and Monica Olvera de la Cruz*,†,‡ †
Department of Materials Science and Engineering, and ‡Department of Chemistry, Northwestern University, Evanston, Illinois 60208, United States S Supporting Information *
ABSTRACT: Crystalline vesicles are promising nanomaterials due to their mechanical stability in various environments. To control their fabrication, it is essential to understand the effects of different experimental conditions on crystallization. Here we perform atomistic molecular dynamics simulations of anionic lipid membranes of 1,2-dilauroyl-sn-glycero-3-phosphol-L-serine. In the presence of Na+ monovalent counterions, we access the phase transition from the liquid-like disordered liquid-crystalline phase to the ordered gel phase by lowering the temperature of the system. The phase transition is conclusively evidenced by the scattering structure factor. Quantitative calculations show that the enhancement of the intertail van der Waals interaction (about −6 kBT) plays a dominant role in driving the phase transition rather than the increase in the cohesive interaction (−0.5 kBT) between lipids and counterions. Meanwhile, in the presence of multivalent counterions of Zn2+ or La3+ the gel phase is found throughout the temperature range investigated. Moreover, the van der Waals interaction per hydrocarbon group is ∼20% stronger in the gel phase (∼ −1.8 kBT regardless of the counterions) than in the liquid-crystalline phase (−1.5 kBT).
■
INTRODUCTION ell membranes are of broad interest due to their critical roles in various biological functions, such as protecting the interior of a cell from the surrounding media and allowing metabolic flow to occur. Apart from the typical spherical structure, crystalline cellular membranes can take the shape of regular and irregular polyhedra including icosahedral buckled viral capsids,1 the square envelope of hallophilic organism,2 and the triangular membrane of archaea and various organelles.3 Micrometer size hollow icosahedral vesicles were artificially fabricated by self-assembling oppositely charged monovalent surfactants in cationic−anionic (catanionic for short) mixtures.4 Recently, buckled vesicles consisting of trivalent cationic and monovalent anionic amphiphiles were synthesized experimentally5 and obtained via analytical methods and Monte Carlo simulations.6 One would also expect that faceted vesicles could be obtained by decreasing the temperature, evidenced by the coarse-grained molecular dynamics (MD) simulations on dipalmitoyl-phosphatidylcholine (DPPC) vesicles.7 Closed polyhedral vesicles are primarily composed of rigid crystalline parts in addition to the minor soft curved parts.7−10 For the purpose of fabricating faceted vesicles with the potential application of use as nanocontainers, it is thus critical to investigate what factors affect the crystallization and how strong their influences are. Even though zwitterionic lipid membranes have been extensively investigated in the last decades (see recent reviews11−13), the polymorphism14 is only occasionally ex-
plored by virtue of simulation approaches. (See refs 15 and 16 and references therein.) The effect of temperature on the liquid-crystalline phase (also known as fluid phase) to the gel phase transition has been investigated using atomistic MD simulations17−19 and coarse-grained MD simulations20 by heating or cooling the systems investigated. Stepniewski et al. compared the liquid-crystalline phase of an unsaturated phosphatidylcholine (PC) bilayer with the gel phase of a saturated PC bilayer at the same temperature.21 Pimthon et al. compared the gel phases of a zwitterionic phosphatidylethanolamine (PE) bilayer with that of a charged phosphatidylglycerol (PG) bilayer.22 The phase transition was shown to also be affected by the concentration of cholesterol in the DPPC or palmitoylsphingomyelin bilayers23 and by the hydration level of glycerol-1-monopalmitin membranes.16 In contrast with the abundant literature on zwitterionic lipid bilayers, ionic lipid membranes are less explored by experiments and simulations13,22,24−30 due to the instability that results from the electrostatic repulsion between the charged headgroups of lipid molecules. Nevertheless, in the presence of multivalent ions, crystalline structures could arise in ionic lipid membranes, owing to the enhanced ionic correlations. Very recently, the crystallization of anionic phosphatidylserine (PS) membranes induced by divalent ions was investigated by means of X-ray
C
© 2013 American Chemical Society
Received: February 19, 2013 Revised: April 3, 2013 Published: April 7, 2013 5073
dx.doi.org/10.1021/jp401767c | J. Phys. Chem. B 2013, 117, 5073−5080
The Journal of Physical Chemistry B
Article
scattering.31 Upon addition of zinc (or calcium), the PS membranes, which are in the liquid-crystalline phase in the presence of monovalent Na+ ions, evolve into the gel phase and eventually the crystalline phase. By virtue of computer simulation approaches, the effects of multivalent ions have been investigated on the structures of charged PG25,30 and phosphatidic acid28 lipid bilayers. However, to the best of our knowledge, the liquid-crystalline phase to gel phase transition of ionic lipid membranes in the presence of multivalent ions has not been analyzed by computer simulations. The aim of the present work is to investigate the phase behavior of anionic lipid membranes by varying the valence of the counterions or the temperature. We aim to quantitatively determine the driving force that leads to the crystallization of anionic lipid membranes. The classical MD simulation at atomistic resolution is employed here, which is capable of providing the packing structures of hydrocarbon tails in lipid bilayers.
each leaflet) is prepared with the proper amount of counterions to neutralize the system (Table 1). See Figure 1b for an example, and Figure S1 in the Supporting Information (SI) for the snapshots of all of the initial structures at 298 K. The initial area per lipid molecule is set to be ∼0.56 nm2 in all of the systems, which is close to the equilibrium value in the liquidcrystalline phase obtained in Figure 2. Such initial bilayer
■
SIMULATION METHODOLOGY Systems under Study. In the present work, bilayers of negatively charged phospholipid of 1,2-dilauroyl-sn-glycero-3phosphol-L-serine (DLPS, Figure 1a) are investigated. These
Figure 2. Calculated area per lipid in the small systems of DLPS-Na+ (317 K), DLPS-Na+ (310 K), DLPS-Zn2+ (310 K), and DLPS-La3+ (310 K). See Table 1 for the compositions. It is indicated that all of the systems become equilibrated after a simulation time of ∼150 ns.
structures shorten the simulation durations for equilibrium structures in atomistic simulations. No additional salt ions are added in the present work. (Note that extra Ca2+ ions have been reported to further tighten palmitoyloleoyl phosphatidylglycerol (POPG) bilayers, whereas extra Na+ ions do not.30) All water molecules are initially distributed outside the hydrophobic interior of the bilayer structures. The ratio of the number of water to lipid molecules is ∼42 (water concentration: ∼54% w/w), which is comparable to that in a previous simulation study on 1,2-dimyristoyl-sn-glycero-3phosphol-L-serine (DMPS).25 Simulation Methodology. Classical MD simulations are performed using the package GROMACS (version 4.5.5).34 The force-field parameters of DLPS lipid, Na+, and Zn2+ are taken from the GROMOS96 54A7 force field.35 The SPC water model, under which the GROMOS96 54A7 force field was parametrized, is employed using the SETTLE algorithm36 to constrain the structure. We use the same set of Lennard-Jones force-field parameters of trivalent ions of La3+ (σ = 0.375 nm, ε = 0.25104 kJ/mol) that were employed to investigate the charge reversal of anionic liposomes,37 obtained by van Veggel and Reinhoudt.38 We use the steepest descent algorithm to equilibrate the initial conformations. Three-dimensional periodic boundary conditions are applied. The electrostatic interactions are calculated using the smooth particle mesh Ewald (PME) method39,40 with a direct space cutoff of 1.2 nm. Lennard-Jones interactions and forces are cut at 1.2 nm, applying the longrange dispersion corrections to the potential energy and the pressure. To speed up the simulations, all covalent bond lengths are constrained using the LINCS algorithm,41 which supports the stable MD integration of a time step of 2.5 fs.42−44 The isobaric−isothermal ensemble (NTP, that is, constant number of particles, temperature, and pressure) is employed, where the temperatures of DLPS, counterions, and water molecules are separately scaled via the Nosé−Hoover thermostat (characteristic time τ = 0.5 ps) in combination with the
Figure 1. (a) Schematic representation of the anionic lipid of DLPS and (b) initial structure of the DLPS-Na+ system, where the Na+ ions are highlighted in blue.
systems are chosen because of the available experimental evidence of their crystallization.31,32 The valence of the counterions is varied starting with monovalent Na+, then divalent Zn2+, and finally trivalent La3+ ions. The temperatures investigated in the simulations range from 298 to 333 K (Table 1). Three initial systems (at 298 K) are built independently using three different random seeds under the package Packmol,33 where a bilayer structure (64 DLPS molecules on Table 1. Compositions of the Small Systems (128 DLPS Molecules) and the Simulation Temperatures Investigated number of components DLPS
Na+
DLPS-Na+
128
128
DLPS-Zn2+ DLPS-La3+
128 128
2
Zn2+
La3+
64 42
water
temperature (K)
5360
298, 310, 317, 323, 333 298, 310, 317, 333 298, 310, 317, 323, 333
5360 5360
5074
dx.doi.org/10.1021/jp401767c | J. Phys. Chem. B 2013, 117, 5073−5080
The Journal of Physical Chemistry B
Article
Table 2. Compositions of the Large Systems (512 or 1152 DLPS molecules) and Some Simulation Data and Results number of components +
DLPS-Na DLPS-Na+ DLPS-Zn2+ DLPS-La3+
DLPS
Na+
512 1152 1152 1152
512 1152
Zn2+
La3+
water
T (K)
duration (ns)
378
21440 48240 48240 48240
317 310 310 310
50 100 100 100
576 18
anisotropic Parrinello−Rahman barostat algorithm (characteristic time τ = 4 ps, reference pressure 1 bar, compressibility 4.5 × 10−5 bar−1). Each system is first equilibrated at 298 K for a duration of 200 ns. The obtained area per lipid values with respect to the simulation time are provided in Figure S2 in the SI, where a transition from a less ordered state (area/lipid > 0.55 nm2) to an ordered state (area/lipid of ∼0.45 nm2) can be inferred. On the basis of the final structures obtained at 298 K, the production simulations at 298 K are performed by extending the simulations for another 200 ns each, and the production simulations at the elevated temperatures (i.e., 310, 317, 323, and 333 K) are performed by increasing the temperature to the desired value and running for a duration of 200 ns each. (One control simulation performed by cooling the system of DLPS-Na+ down from 333 to 298 K will be discussed later.) For the data collection and analysis of area per lipid and tail thickness, only those frames in the last 50 ns are taken into consideration. It is noteworthy that simulation longer than 250 ns has been claimed to be necessary for precise calculation of the area per lipid for POPG lipid bilayers.30 In the present work, the simulation for each system is first equilibrated at 298 K for a duration of 200 ns and then heated to the desired temperature for another 200 ns. To further validate the simulations, we extended four simulations of the systems of DLPS-Na+ (317 K), DLPS-Na+ (310 K), DLPS-Zn2+ (310 K), and DLPS-La3+ (310 K) to a simulation duration of 400 ns. The calculated area per lipid values with respect to the simulation time (Figure 2) indicates that a simulation duration of 150 ns is long enough for the systems to equilibrate. Scattering Structure Factor. To better calculate the scattering structure factor, we enlarged the four systems, that is, DLPS-Na+ (317 K), DLPS-Na+ (310 K), DLPS-Zn2+ (310 K), and DLPS-La3+ (310 K), by stacking 2 × 2 pre-equilibrated structures at 400 ns side by side in the X × Y dimensions for the DLPS-Na+ (317 K) system (or 3 × 3 for the other systems). See Table 2 for the compositions of the large systems. The simulations on these large systems are then performed for the durations indicated in Table 2. The same simulation parameters are employed for these large systems. The obtained area per lipid as a function of the simulation time (Figure S3 in the SI) supports that the systems are well-equilibrated. The calculated area per lipid values and the hydrophobic tail thicknesses, presented in Table 2, are also in good agreement with the corresponding results in Figure 3 in small systems. (See also Table S1 in the SI.) The consistency indicates that the finite size effect is negligible in the simulations presented here.
area/lipid (nm2) 0.562 0.454 0.440 0.434
± ± ± ±
0.008 0.001 0.001 0.001
tail thickness (nm) 2.30 2.80 2.86 2.88
± ± ± ±
0.03 0.01 0.01 0.01
Figure 3. Calculated area per lipid as a function of the simulation temperature in the three systems with different counterions of Na+, Zn2+, or La3+. The dotted lines are to guide the eye. The standard deviations are provided as error bars. The inset shows the corresponding tail thicknesses of the bilayers.
temperature range 310−317 K. Moreover, such a phase transition in DLPS-Na+ is supported by the obtained bilayer tail thickness, as illustrated in the inset of Figure 3. The tail thickness is calculated based on the carbon atoms in the ester groups in the opposite leaflets. In the following section, we demonstrate that the DLPS-Na+ system is in the tilted gel phase (Lβ′) at low temperatures (≤310 K) and in the liquid-crystalline phase (Lα) at higher temperatures (≥317 K). The experimental main phase transition temperature (from the liquid-crystalline phase to gel phase) of DLPS with the counterion of sodium was reported to be 287 K using differential scanning calorimetry45 or 291 K using X-ray scattering.32 Therefore, for the DLPS-Na+ system, the main phase transition temperature obtained in our simulations is about 20−30 K higher than the experimental data. Similar discrepancies have been observed in previous simulation studies,17−20 which are probably attributed to the force field employed15,35 and the hydration degree16,46 of the lipid membranes. In contrast, no phase transition is observed for the DLPS-Zn2+ or DLPS-La3+ membranes in the temperature range of 298−333 K. The wide-angle X-ray scattering (WAXS) experimental area per lipid of DLPS (with the counterion of Na+) was reported to be ∼0.603 nm2 at 293 K in the Lα phase, with a slight increase to 0.606 nm2 as the temperature is elevated to 333 K.32 Consistently, the area per lipid at 333 K in the present work is calculated to be 0.59 nm2, deviating −3% from the experimental value at the same temperature. The dependence of the area per lipid on temperature is also reproduced here, even though the slope of ∼0.18 Å2/K calculated from our simulations is higher than the experimental value of ∼0.08 Å2/K.32 In comparison, the influence of temperature on the area per lipid of the systems in the gel phase is negligible irrespective of the valence of the counterions. Meanwhile, using the WAXS experiments,31 the gel phase of DLPS with the counterions of Zn2+ was reported to possess an area per lipid of 0.4176 nm2 (0.3834
■
RESULTS AND DISCUSSION Phase Transition: Area Per Lipid. The lateral area per lipid, which is frequently employed to validate computer simulations and can also be used to elucidate the phase behavior, is calculated and shown in Figure 3. A sharp change in the area per lipid is observed in the DLPS-Na+ system in the 5075
dx.doi.org/10.1021/jp401767c | J. Phys. Chem. B 2013, 117, 5073−5080
The Journal of Physical Chemistry B
Article
nm2 for the crystalline phase, LC) at room temperature. Our simulation on the DLPS-Zn2+ system at 298 K results in an area per lipid of 0.432 ± 0.002 nm2, which deviates ∼3% from the experimental datum in the gel phase. To validate the area per lipid values in Figure 3 resulting from the increasing temperature method in our simulations, we performed a control simulation, where the DLPS-Na+ system is cooled from 333 to 298 K and then equilibrated for another 200 ns. In the control simulation, we find an area per lipid of 0.459 ± 0.003 nm2, which is in good agreement with the value of 0.448 ± 0.004 nm2 in Figure 3. Such consistency suggests the reproducibility of our data in Figure 3 and indicates that the experimental area per lipid in the gel phase of 0.599 nm2 at 283 K32 in the presence of Na+ counterions is unconvincing. Moreover, as shown in the inset of Figure 3, the tail thickness of the DLPS-Na+ system is ∼2.8 nm in the gel phase and 2.2 to 2.3 nm in the liquid-crystalline phase, which is in good agreement with the experimental values32 of 3.0 nm (gel phase) and 2.4 to 2.6 nm (liquid-crystalline phase), respectively. The last frames from the 400 ns simulations of the systems DLPS-Na+ (317 K), DLPS-Na+ (310 K), DLPS-Zn2+ (310 K), and DLPS-La3+ (310 K) are shown in Figure 4. The liquidcrystalline phase of the hydrocarbon tails is evidently illustrated in the system of DLPS-Na+ (317 K) (Figure 4a) and the ordered gel phase in the other systems (Figure 4b−d). It is necessary to note that a mismatch exists between the
orientation of the crystalline structure in the upper leaflet and that in the lower leaflet. A movie where the DLPS-Na+ membrane at 310 K is rotated around the normal is provided in the SI, which undoubtedly shows the coexistence of the crystalline structures in the upper and lower leaflets in the gel phase. To our knowledge, the physics of such a mismatch between the crystalline orientations in the upper and lower leaflets, which may originate from the correlation between the opposite leaflets, has never been discussed. The correlation between the opposite leaflets has been indicated by the transbilayer complementarity phenomenon (i.e., a preference for the pairing of a long phospholipid with a short phospholipid across the bilayer) evidenced by virtue of nearest-neighbor recognition experiments48−51 and simulation work52 in the coarse-grained resolution. Figure 4 also indicates that while some monovalent sodium ions (∼10% from the density distribution in Figure S4 in the SI) are distributed in the aqueous phase, approximately zero ions are in the aqueous phase in the Zn2+ and La3+ cases due to the competition between the hydration ability and the electrostatic interactions. That is, a higher valence of the counterions enhances the electrostatic attraction between the oppositely charged groups on DLPS and the counterions and decreases the hydration degree of the charged groups on DLPS (Figure 5). Similar findings on the distributions of Na+ and Ca2+ in the vicinity of DMPS25 or POPG30 bilayers have been reported.
Figure 5. Radial distribution function, g(r), of water molecules around the COO− groups of DLPS molecules. The decrease in the dehydration degree of the COO− groups is induced by either the transition from the liquid-crystalline phase (317 K) to the gel phase (310 K) for the Na+ counterions or the presence of a higher valence of the counterions.
Phase Transition: Scattering Structure Factor. Strictly speaking, the calculation of the area per lipid can only suggest the occurrence of a phase transition. Obtaining conclusive evidence of the liquid-crystalline phase to gel phase transition requires the calculation of the scattering structure factor, which is absent in previous simulation studies,17−22 due to the necessity of simulating a large, computationally expensive system. The small q range of the scattering structure factor stands for the long-range correlations, such as multilamellar structure, which is beyond the scope of the present work. The structure factors in the large-q range, which will be discussed here, reflect the packing behavior at molecular resolution. The liquid-like disordered packing of the hydrocarbon tails in liquidcrystalline phase (Lα) is indicated by a broad peak; the 2D
Figure 4. Snapshots of the last simulation frames of (a) DLPS-Na+ (317 K), (b) DLPS-Na+ (310 K), (c) DLPS-Zn2+ (310 K), and (d) DLPS-La3+ (310 K). Order in the orientation of the hydrocarbon tails is evidently observed in panels b−d but not in panel a. The counterions are highlighted using the VDW drawing method under VMD,47 with waters omitted for clarity. 5076
dx.doi.org/10.1021/jp401767c | J. Phys. Chem. B 2013, 117, 5073−5080
The Journal of Physical Chemistry B
Article
scattering peaks are in good agreement with the experimental datum of q(0,1) = 14.78 nm−1 in the Lβ phase for the DLPS-Zn2+ system.31 Moreover, there exist two additional peaks located at q ≈ 25.3 and 29.3 nm−1, which correspond to the second (√3q) and third (2q) peaks of a hexagonal packing structure, respectively.14 The order parameter15,16 (Figure 8) and the tilt angle (Figure 9) of the hydrocarbon tails are calculated to further
hexagonal packing in the gel phase (untilted Lβ, or tilted Lβ′) leads to one sharp principal peak at q, followed by some weak peaks at √3q and 2q; the crystalline phase (LC) is supported by multiple principal peaks in the large-q range of the scattering structure factor.14 In the calculation of the scattering structure factor, a higher resolution in the reciprocal space requires a larger simulation cell. Four systems are enlarged to that end: DLPS-Na+ (317 K), DLPS-Na+ (310 K), DLPS-Zn2+ (310 K), and DLPS-La3+ (310 K) (Table 2). Hereafter, only the results from these large systems are discussed. On the basis of the large systems, the scattering structure factors of the hydrocarbon atoms, which are obtained via the Fourier transformation of the corresponding radial distribution function (RDF, g(r)) of the hydrocarbon atoms (Figure 6), are calculated and plotted in Figure 7.
Figure 8. Order parameter, SCD, of the hydrocarbon tails as a function of the hydrocarbon number in the four large systems. The obtained order parameters in the liquid-crystalline and the gel phases are consistent with the corresponding ones reported in ref 16.
Figure 6. Radial distribution function, g(r), of the hydrocarbon atoms in the hydrophobic bilayer core region in the four large systems. The fluctuations of g(r) in the last three systems support the existence of ordered packing structures of the hydrocarbon tails.
Figure 9. Tilt angle of the hydrocarbon chains in the large systems, which is the angle between the bilayer normal and the vector defined by the second and the last hydrocarbon atoms from the corresponding ester group. The tilt angles of ∼10° suggest the tilted gel phase (Lβ′) in these systems.
characterize the packing structures in the gel phase. The order parameters of the methylene groups in the gel phase are evidently higher than those in the liquid-crystalline phase, depicting the highly ordered orientation of the hydrocarbon tails in the gel phase. The obtained order parameters are consistent with the corresponding values in the liquidcrystalline and the gel phases in the glycerol-1-monopalmitin membranes reported in ref 16. Moreover, the consistency between the order parameters in the upper and lower leaflets also suggests that both the upper and lower leaflets are in the same phase. (See Figure S5 in the SI for an example of the DLPS-Na+ system at 310 K.) The tilt angle of the hydrocarbon tails is calculated to be ∼10°, indicating the tilted packing in the gel phase (Lβ′). Note that even though the untilted packing structure (Lβ) was claimed for the DLPS-Zn2+ system,31 the tilted packing structure (Lβ′ ) was not experimentally excluded.53
Figure 7. Log−log plot of the scattering structure factors of the hydrocarbon atoms on the DLPS tails. The peaks of q ≈ 14.6, 25.3, and 29.3 nm−1 in the last three systems indicate the hexagonal packing (see inset) of the hydrocarbon tails in the ordered gel phase, whereas the liquid-crystalline phase (Lα) in DLPS-Na+ (317 K) is supported by the broad principal peak and the lack of the corresponding second and third peaks.14
Because of the limited length of the simulation cell in the present atomistic simulations (
