Article Outline

Article Information

PubMed

Related Articles

  • …more

Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 92, Issue 4, 1114-1124, 15 February 2007

doi:10.1529/biophysj.106.086272

Biophysical Theory and Modeling

Atomic-Scale Structure and Electrostatics of Anionic Palmitoyloleoylphosphatidylglycerol Lipid Bilayers with Na+ Counterions

Wei Zhao*, Tomasz Róg*, ‡, Andrey A. Gurtovenko§, Ilpo Vattulainen†, , ‖ and Mikko Karttunen**, Go To Corresponding Author,  

* Biophysics and Statistical Mechanics Group, Laboratory of Computational Engineering, Helsinki University of Technology, Espoo, Finland
† Laboratory of Physics and Helsinki Institute of Physics, Helsinki University of Technology, Espoo, Finland
‡ Department of Biophysics, Jagiellonian University, Kraków, Poland
§ Computational Laboratory, Institute of Pharmaceutical Innovation, University of Bradford, Bradford, West Yorkshire, United Kingdom
Institute of Physics, Tampere University of Technology, Tampere, Finland
‖ MEMPHYS–Center for Biomembrane Physics, University of Southern Denmark, Odense, Denmark
** Department of Applied Mathematics, The University of Western Ontario, London, Ontario, Canada

Address reprint requests to Mikko Karttunen.

Abstract

Anionic palmitoyloleoylphosphatidylglycerol (POPG) is one of the most abundant lipids in nature, yet its atomic-scale properties have not received significant attention. Here we report extensive 150-ns molecular dynamics simulations of a pure POPG lipid membrane with sodium counterions. It turns out that the average area per lipid of the POPG bilayer under physiological conditions is ∼19% smaller than that of a bilayer built from its zwitterionic phosphatidylcholine analog, palmitoyloleoylphosphatidylcholine. This suggests that there are strong attractive interactions between anionic POPG lipids, which overcome the electrostatic repulsion between negative charges of PG headgroups. We demonstrate that interlipid counterion bridges and strong intra- and intermolecular hydrogen bonding play a key role in this seemingly counterintuitive behavior. In particular, the substantial strength and stability of ion-mediated binding between anionic lipid headgroups leads to complexation of PG molecules and ions and formation of large PG-ion clusters that act in a concerted manner. The ion-mediated binding seems to provide a possible molecular-level explanation for the low permeability of PG-containing bacterial membranes to organic solvents: highly polar interactions at the water/membrane interface are able to create a high free energy barrier for hydrophobic molecules such as benzene.

Introduction

Phosphatidylglycerols (PGs) are among the most abundant lipids in nature. In higher plants, they constitute ∼20%–30% of lipids 1. They are also one of the major constituents of bacterial membranes 2; mixed PG/phosphatidylethanolamine (PE) lipid bilayers have been used as models for bacterial membranes 3,4,5,6. In contrast, the concentration of PG lipids in animal cell membranes is low. In red blood cells PGs constitute ∼2% of all phospholipids 7. Though, whereas being less common in animal membranes, they are present in mitochondria and some tissues. From the physical point of view, PGs, together with phosphatidylserine (POPS) and cardiolipin, are anionic, the former two carrying a unit positive charge and the latter being divalent under physiological conditions.

The nonzero net charge of PGs makes them different from commonly studied zwitterionic (neutral) lipids. Since direct electrostatic interactions between charged species are rather strong, PGs, as well as other charged lipids, are considered membrane stabilizers and destabilizers 8 and are believed to play an important role in controlling membrane-peptide/protein interactions. For this reason, PGs have been suggested as one of the possible components in immunoliposomes, generic drug delivery vehicles 9. It is also noteworthy that PGs constitute up to 80% of membrane lipids of Staphylococcus aureus, a common human pathogen with its methicillin-resistant strain (MRSA)—the so-called hospital bacteria—being a serious problem in hospitals around the world. Hence, it has been suggested that electrostatics may play an important role in membrane selectivity of bacteria. Atomic-scale details of these electrostatic interactions would, therefore, be very helpful and important for a better understanding of how such a selectivity can be controlled to make drugs more efficient.

Furthermore, studies of toxic effects of organic solvents on bacteria have revealed that relative concentrations of PG and PE lipids in their membranes could change depending on the chemical character of solvent molecules 10,11. Such an alteration in headgroup composition seems to be a means to change permeability of a membrane and to preserve its stability. Precise mechanisms responsible for solvent-induced effects are not known, due to the lack of detailed microscopic information.

Molecular dynamics (MD) simulations provide a detailed microscopic picture of atomic-scale interactions and processes, which are not easily accessible by experimental methods. Systems having an abundance of charged components, such as anionic or cationic membranes, and systems with salt ions are computationally demanding—first reliable MD studies in the field were reported only recently: anionic POPS lipid bilayers 12,13, PG monolayers 14, mixed cationic/zwitterionic lipid membranes 15,16,17, phosphatidylcholine lipid membranes with salt 18,19,20,21, and ion leakage through water pores in PC membranes 22.

Given the abundance of PGs in nature, it is somewhat surprising that, to the best of our knowledge, there has been only one atomistic computational study of pure PG lipid bilayers 23 focusing on the hydrogen bonding between the lipids. All related simulation studies to date have dealt with mixtures of PGs with other lipids 3,6 or with PG monolayers at the air/water interface 14. On the experimental side, pure POPG bilayers, or POPG mixed with PC at high concentrations, have been studied (to our knowledge) only by Cowley et al. 24 in the 1970s and more recently by Pozo-Navas et al. 25,26. The phase behavior 27,28 and various aspects of the structure-function relationship of POPGs in PE bilayers 8,29,30,31,32 have been experimentally characterized, though.

In this work we employ atomistic MD simulations to gain a thorough understanding of the properties of PG-rich membranes through a systematic comparison with corresponding PC bilayers. This is needed to establish the fundamental molecular and atomistic level physical origin of the observed phenomena. This knowledge should be useful also in interpreting experimental findings, and it offers a possible explanation as to why nature has chosen to use charged lipids structures such as bacterial membranes.

We focus on the structural and electrostatic properties of single-component palmitoyloleoylphosphatidylglycerol (POPG) lipid membranes. Zwitterionic palmitoyloleoylphosphatidylcholine (POPC) was used as a control as its properties are well known and, importantly, since PCs are the most abundant lipids in biological membranes. As they have a dipole in their headgroup, their behavior is expected to be quite different from that of anionic PGs. Indeed, when compared to a POPC bilayer, the average area per lipid in POPG bilayer simulations turns out to be much smaller. That seems to be characteristic for charged bilayer as similar compression has also been seen in other anionic bilayers (see, e.g., Petrache et al. 33). Our simulations reveal the physical mechanism responsible for that, and the key results of this work are the following: there are strong attractive interactions between anionic POPG lipids, which overcome the electrostatic repulsion between negative charges of PG headgroups. We show that interlipid counterion bridges and strong intra- and intermolecular hydrogen bonding lead to this seemingly counterintuitive behavior. In particular, the strength and stability of ion-mediated binding between anionic lipid headgroups leads to complexation of PG molecules and ions and formation of large PG-ion clusters that act in a concerted manner. The ion-mediated binding may provide a possible molecular-level mechanism for the low permeability of PG-containing bacterial membranes to organic solvents: highly polar interactions at the water/membrane interface are able to create a high free energy barrier for hydrophobic molecules such as benzene.

The rest of this work is organized as follows: in the next section we present the model and simulation details. Then, we first show the results for the basic structural quantities such as the area per lipid, density profiles, and order parameters. After that we come to the main results: hydrogen, ion bonding, and their consequences such as clustering. Finally, we finish with conclusions and a discussion.


Model and simulation details

We have studied a lipid bilayer system consisting of 128 anionic POPGs and 3527 water molecules. A total of 128 Na+ counterions was included to ensure electroneutrality. The force field parameters for POPG (Fig. 1) were taken from the united atom force field of Berger et al. 34. They are based on the GROMOS force field for lipid headgroups, the Ryckaert-Bellemans potential functions 35,36 for hydrocarbon chains, and the optimized parameters for liquid simulations parameters 34 for the Lennard-Jones interactions between united CHn groups of nonpolar lipid chains. Partial charges for the glycol group (H16-O16-C13-C12-O15-H15 in the headgroup of POPG; see Fig. 1) were taken from 1,2-propanediol (D. Warren, unpublished). The force field parameters for POPG are available online at http://www.softsimu.org/downloads.shtml. The SPC model 38 was used for water, and standard GROMACS force field parameters were employed for sodium ions. For subtleties related to modeling of Na+, please see Patra and Karttunen 39.

Display large version of this figure
Figure 1
Chemical structures and the numbering of atoms of POPG and POPC molecules used in this study.

All Lennard-Jones interactions were cut off at 1.0nm without shift or switch functions. Lipid bonds were constrained using the LINCS algorithm 40, and the SETTLE algorithm 41 was used for water molecules. The electrostatic interactions were handled by the particle mesh Ewald (PME) method 42,43, which has been shown to perform well in membrane simulations 44,45,46.

The simulations were performed in the NpT ensemble with a time step of 2fs, using the GROMACS simulation package 47. The temperature was kept constant at T=310K using the Berendsen thermostat 48 with a coupling time constant of 0.1ps. Lipid molecules and solvent (water molecules and sodium ions) were separately coupled to a heat bath. The Berendsen barostat 48 with a coupling constant of 1.0ps was employed to keep pressure constant (1bar). The pressure coupling was applied semiisotropically: the box size in the direction of the bilayer normal (z axis) and the cross sectional area (xy-plane) were allowed to vary independently.

The starting structure of a POPG bilayer was prepared using an equilibrated POPC bilayer 49: the choline moieties of POPC lipids were replaced with the glycerol groups such that equal numbers of D-POPG and L-POPG lipids were created, leading to a lipid bilayer system with neutral chirality (racemic). The two cases differ with respect to the chirality of the C12 site 50, and hence this choice guarantees that our lipid bilayer system is characterized by neutral chirality. After initial energy minimization, a short 10-ps run in the NVT ensemble was performed to remove unphysical voids due to structural modifications. The resulting POPG system was simulated for 150ns. By analyzing the time evolution of the area per lipid and the coordination of Na+ counterions with lipids and water, we concluded that the system had equilibrated within 70ns (see Fig. 2). The time range from 70 to 150ns was therefore used for analysis. In total, the MD simulations took ∼12,000 CPU hours on AMD Opteron 2.2GHz computers (AMD, Sunnyvale, CA).

Display large version of this figure
Figure 2
Time evolution of the area per lipid for a POPG bilayer. Dashed line shows the average value in equilibrium, 〈A〉.

Below, the properties of POPG are compared with those of POPC. Hence, for the purpose of comparison, a complementary simulation of a POPC lipid bilayer containing 128 POPC and 3655 water molecules was also conducted at a temperature of 310K for 30ns. Our choices for the force field and the initial (and already equilibrated) configuration are the same as in our previous simulations 49. For the given temperature and hydration, both studied systems, POPG and POPC, are in the fluid (liquid-disordered) phase. The experimenatlly determined main transition temperatures are below 273K for both POPG (≈ 273K) 26,27 and POPC (268K) 51. Our data (in the following sections) are also consistent with the liquid-disordered phase.


Results

Area per lipid

The average area per lipid, 〈A〉, is perhaps the most widely used quantity to characterize lipid bilayer systems, since it affects a variety of physical processes such as lateral diffusion, membrane elastic properties, and permeation. The time evolution of the area per lipid is shown in Fig. 2. The long equilibration time visible in the figure is mainly due to the slow kinetics associated with the binding of counterions with the charged POPG lipids and in part due to hydrogen-bond formation among the lipids. Similar findings have also been made in simulations of anionic POPS bilayers 13 and in bilayer simulations under the influence of salt 17,18,20,21.

For the average area per lipid in a POPG bilayer we obtain 〈APOPG〉=0.530±0.006nm2. This is a remarkably small number compared to typical values in corresponding zwitterionic lipid systems. For example, for a POPC bilayer studied in this work we found 〈APOPC〉=0.658±0.009nm2 at the same temperature. Sphingomyelin bilayers in turn have yielded average areas per lipid of the order of 0.52nm2, and the lipid hydrocarbon chains in these systems are known to be highly ordered. Similar enhanced ordering is therefore expected for POPG too (see below).

Comparison of our results to experimental data is difficult since, to our knowledge, there are no experimental estimates for the average area per lipid in single-component POPG bilayers. However, comparison to previous related simulations is possible. Murzyn et al. 6 used the Voronoi tessellation technique to extract the average area per lipid from simulations of a POPE/POPG mixture. They reported a value of 0.628±0.003nm2 for the average area per POPG, which is considerably larger than the value found in our work. One has to note, however, that the result by Murzyn et al. cannot be directly compared to our findings. This stems from the fact that they applied the Voronoi tessellation technique to a binary system, and then there is no unique way to decompose the free area (volume) among the different lipid components. Hence the areas of two-dimensional Wigner-Seiz cells in binary systems do not correspond to the area per lipid measured in single-component lipid bilayers. This issue has been extensively discussed in computational studies of cholesterol in phospholipid bilayers 52,53,54,55 and also in experiments on mixtures of nonsterol natural lipids 56.

A better, but also indirect comparison can be provided by MD simulations of pure POPS bilayers 13. Although headgroups of POPG and POPS lipids have different chemical structures, both of these lipids are anionic under physiological conditions. Mukhopadhyay et al. obtained 0.55±0.01nm2 for the average area per POPS under simulation conditions similar to those employed in this study—close to 0.53nm2 found here for a POPG bilayer. These findings are also in agreement with Cowley et al. 24, who studied PC/POPG bilayers at different concentrations under varying hydration conditions.

We also computed the autocorrelation function to measure the characteristic timescales for area fluctuations of POPGs (data not shown). As expected on the basis of the above discussion, they turned out to be much longer than those observed for POPC bilayers 49. The above observation indicates that an anionic POPG bilayer is less fluid than one would naïvely expect 27: its densely packed structure, which seems to be counterintuitive at first sight, is caused by tight binding of anionic lipids through relatively stable ion bridges (see the section Interactions of Na+ ions with water/membrane interface).


Density profiles

Next, we consider the density profiles of various bilayer constituents (see Fig. 3).

Display large version of this figure
Figure 3
Component-wise density profiles across the bilayer as a function of distance from the membrane center (z=0). Gly: the glycol group (H16, O16, C13, C12, O15, H15); Phos: the phosphate group (O11, O12, O13, O13, P); and Carb: the carbonyl ester group (C21, O21, O22, C31, O31, O32).

The maxima of the (mass or electron) density profiles are often used to estimate the thickness of the bilayer 6,44,57. As in previous studies, we considered the distance between the average positions of phosphate atoms in the two leaflets of a lipid bilayer. That distance computed from the electron density profile yielded 4.39±0.02nm for the phosphate-phosphate distance dP-P. The result is in line with the previous simulation study of a POPG bilayer at 310K 23.

The maximum of the mass density profile of POPG lipids turns out to be located at the average position of the ester groups instead of the phosphate groups as has been found for phosphatidylcholine bilayers 49. The most remarkable feature is the distribution of Na+ ions characterized by two peaks in the interface region (see Fig. 4). The first (main) peak is located at ≈1.9nm from the bilayer center and overlaps with the peak due to the POPG carbonyl ester groups. Elmore has recently made a similar observation 23. The second, much smaller peak is at ≈2.65nm. Further analysis shows that this peak corresponds to Na+ ions mainly located in the water phase, close to the membrane but beyond the phosphate groups, although the distribution of Na+ is strongly affected by the surface potential of the negatively charged surface of a POPG bilayer. Hence the positively charged sodium ions tend to avoid being in close contact with the negatively charged phosphate groups and instead prefer interactions with the ester groups.

Display large version of this figure
Figure 4
Scaled number densities for various constituents of the POPG system. The number density of lipids is not scaled, whereas the densities of water, phosphorus, O22, O32, and Na+ ions are scaled by dividing them by 70, 10, 10, 10, and 7, respectively. The bilayer center is at z=0.

The peaks of the density profiles of phosphodiester oxygens (O13, O14), hydroxyl oxygens O16 and O15 (see Fig. 1) are located at ≈2.2nm, 2.05nm, and 2.0nm, respectively (data not shown). When these numbers are compared to the distribution of Na+ ions, they suggest that sodium ions interact more favorably with ester groups rather than phosphodiester moieties. Meanwhile, sodium ions corresponding to the second, smaller peak of the distribution (≈2.65nm) are not bound to lipids (see the section Interactions of Na+ ions with water/membrane interface for a detailed discussion). Due to strong electrostatic interactions, all the Na+ ions stay near lipid headgroups and interact with various oxygen atoms. Yet the mutual electrostatic repulsion between the Na+ ions together with their thermal motion and energy compensation from hydration competes with the tendency of sodium ions to reside close to the headgroups. The two effects together lead to the existence of the broad second peak in the distribution of the Na+ ions.


Headgroup orientation

Due to the glycerol group, it is not straightforward to characterize the orientation of PG headgroups as compared to that of phosphatidylcholine headgroups defined by their P-N dipoles. Hence, for POPG we considered distributions of the P-C12, P-O15, and P-O16 vectors (see labeling in Fig. 1) as shown in Fig. 5.

Display large version of this figure
Figure 5
The probability distribution P(α) of the angle α between PG headgroup vectors and the outward bilayer normal.

The average orientation angles between the PG headgroup vectors (P-C12, P-O15, and P-O16) and the outward bilayer normal are found to be 85±2°, 98±3°, and 89±6°. Consequently, PG headgroups are on average almost parallel to the membrane surface. This is also confirmed by the density profile shown in Fig. 3, where the maxima of phosphate and glycol groups overlap to a large extent. As for the P-O15 vector, the average inward orientation possibly results from the intramolecular hydrogen bond O15-H15-O12.

The orientation of the PG headgroup with respect to the membrane surface has been examined earlier by small angle x-ray 58 and neutron scattering techniques 59. These experiments were performed in the ordered phase of DPPG or Escherichia coli PG, respectively, at very low hydration (nwater∼2.6 per lipid). An approximate tilt of 30° was derived from these data for the PG headgroup, though an alternative interdigitated structure could not be strictly excluded 59. The above experimental results are not comparable to our simulations as they were performed in the gel phase, and the former one even under varying pH. In recent simulations in the fluid phase, Elmore 23 found an average angle of 102° with a broad variation.

It has been observed that the PG headgroup undergoes changes of orientation as the hydration level is changed 60. This orientation flexibility, as suggested by Kurze et al., is probably related to the breakdown of the continuous primary hydration shell of the PG headgroup, accompanied by the formation of an interheadgroup hydrogen-bonding network that replaces the headgroup-water interaction. Unfortunately, experimental data for the orientation of PG headgroup are not available at the hydration level we applied in our bilayer system.


Ordering of hydrocarbon tails

The ordering of lipid acyl chains is typically characterized by the order parameter tensor

(1)
where α, β=x, y, z, and θα is the angle between the αth molecular axis and the bilayer normal (z axis). The order parameter is then computed separately for all carbons along the acyl chains. The above is related to the experimentally measured deuterium order parameter
(2)

We therefore present our results in terms of |SCD|. Since we used a united atom model and therefore had no explicit hydrogen positions, they had to be reconstructed assuming a perfect tetrahedral arrangement.

The order parameters for saturated sn-1 chains and unsaturated sn-2 chains are presented separately in Fig. 6. For comparison, the order parameters from a reference simulation for a POPC bilayer are presented as well. Fig. 6 shows that POPG lipids are more ordered than lipids in a POPC bilayer and that the qualitative features of the order parameter profiles are quite different. Unlike for POPC, the order parameter of POPG acyl chains increases up to the eighth carbon along the chain. Eventually it decreases for the last few carbons but, nevertheless, the POPG chains are highly ordered even at their ends. The higher ordering is in line with what one can expect from the behavior of the area per lipid (see Fig. 2).

Display large version of this figure
Figure 6
Deuterium order parameters SCD of saturated sn-1 (top) and unsaturated sn-2 (bottom) chains in POPG and POPC bilayers.

Hydrogen bonding

Experimental observations have suggested that the hydroxyl group of a PG lipid has the potential to form both intra- and intermolecular hydrogen bonds 61. To determine whether these hydrogen bonds exist, we follow the generally accepted approach and define two geometrical parameters: the distance r between the hydrogen and the acceptor, and the donor-hydrogen-acceptor angle (zero extended) α. The conditions r≤0.25nm and α≤60° were then used as criteria for a hydrogen bond (H-bond) to exist.

On the average, we found that O16 forms ∼0.15 intramolecular and ∼0.47 intermolecular hydrogen bonds with various lipid oxygen atoms per POPG molecule. The intermolecular hydrogen bonds between the O16-H16 groups and carbonyl groups are the most abundant hydrogen bonds found in our study. We will discuss this issue in more detail below. Table 1 summarizes the average numbers of intra- and intermolecular hydrogen bonds.

Table 1 Average number of intermolecular and intramolecular hydrogen bonds in a POPG lipid molecule. The errors given are standard error estimates. See Fig. 1 for atom labeling.
H-donorH-acceptorIntermolecularIntramolecular
O15O12–1.00±0.01
O16O110.06±0.02–
O16O120.01±0.010.02±0.01
O16O13, O140.07±0.02–
O16O150.05±0.020.05±0.02
O16O160.02±0.01–
O16O210.07±0.010.04±0.01
O16O220.14±0.020.02±0.01
O16O310.01±0.010.01±0.01
O16O320.05±0.020.01±0.01

Interactions of Na+ ions with water/membrane interface

Radial distribution functions and coordination numbers

Fig. 7 shows the radial distribution functions (RDFs) of various oxygen atoms (including intermolecular and intramolecular parts) with the hydrogen atom H16 in the POPG headgroup. The results for H15 are not shown since it forms intramolecular hydrogen bonds with the phosphoether oxygens. In Fig. 7, the first minima in the RDF for all the carbonyl and phosphodiester oxygen atoms are located at ∼0.24nm, indicating hydrogen bonding with H16 hydrogen atoms. This confirms intermolecular hydrogen bonding in the POPG bilayer.

Display large version of this figure
Figure 7
RDFs of various oxygen atoms with the hydroxyl hydrogen H16; see Fig. 1 for atom labeling.

The distribution of Na+ in Fig. 8 suggests that a large amount of sodium ions are located in the ester region of the water/membrane interface, rather than in the phosphate region. Clearly, there are prominent peaks with oxygen atoms O22 and O32 associated with the ester bonds, followed by O15 in the glycerol part of the headgroup. The remaining peaks related to O16 (also part of the glycerol group) and oxygen atoms O13 and O14 in the phosphate group are considerably less significant.

Display large version of this figure
Figure 8
RDFs of various oxygen atoms with Na+ ions.

The structure and population of oxygen atoms around the Na+ ions can also be characterized by coordination numbers 39 computed by integration of the corresponding RDF over the first coordination shell,

(3)
where VNa+ is the van der Waals volume of a sodium ion and rmin is the position of the first minimum of the RDF. This coordination number describes the number of atoms of interest in the first coordination shell of a Na+ ion.

From Eq. (3) one finds, on average, 3.36, 0.91, 0.60, 0.40, 0.19, and 0.04 oxygen atoms from water, O22, O32, O15, O16, and from phosphodiester oxygen atoms, respectively, around the Na+ ion.

To get a more detailed view on how Na+ ions interact with different charged groups, in Fig. 9 we show the coordination numbers of Na+ ions with various oxygen atoms as a function of the distance from the bilayer center.

Display large version of this figure
Figure 9
Average coordination numbers of Na+ ions with various lipid oxygens, water oxygens, as well as the total coordination numbers (including contributions from all oxygen atoms in POPG) as a function of the distance from the bilayer center.

The total coordination number shows a small increase, from ∼5.0 deep inside the lipid bilayer to ∼5.7 in the aqueous phase. The coordination numbers of Na+ ions with lipid oxygen atoms show a maximum in the ester region of the water/membrane interface, indicating strong interactions between Na+ ions and carbonyl groups. Due to binding of sodium ions to the membrane, water molecules are squeezed out from the first coordination shell.

In the section Density profiles we found that the density profile of sodium ions was characterized by two peaks in the membrane-water interfacial region. Fig. 10 elucidates the origin of this behavior: the first peak in the density profile of sodium ions around z≈1.9nm is due to ions bound to lipid oxygens and especially those bound to the oxygens in the ester groups. The second peak at z∼2.65nm, in turn, arises from those sodium ions that are not bound to any lipid oxygens but are rather confined to water molecules in the vicinity of the membrane-water interface.

Display large version of this figure
Figure 10
Distributions of Na+ bound and unbound to lipid oxygens as functions of distance from the bilayer center.

Ion bridges and ion-lipid clusters

The term “ion bridge” refers to the case in which two lipids bind to each other via a Na+ ion. To characterize the existence of ion bridges, we define an ion-lipid bond to occur when the carbonyl oxygen is within a distance of 0.33nm from a Na+ ion, which is the radius of the first shell of Na+ ions.

Here, we only define lipid carbonyl oxygens (O22, O32) instead of all lipid oxygens as binding sites available for ion bridges 19. This is justified since the bonding between Na+ ions and carbonyl oxygens is much stronger compared to interactions between Na+ ions and other lipid oxygens. These ion bridges constitute the majority of all interactions between Na+ ions and lipid, as seen above.

When analyzing the trajectory we found, on average, 193±5 ion-lipid bonds (in the model of 128 PGs and 128 sodium ions), involving 63% of Na+ ions and 97% of lipids. The percentages for Na+ ions binding to one, two, three, and four lipids are, respectively, 9%, 24%, 24%, and 6%. Fig. 11 illustrates some of the cases. As for lipids, the percentages binding to one, two, and three Na+ ions are, respectively, 45%, 51%, and 1%.

Display large version of this figure
Figure 11
Snapshots of lipids binding to Na+ ions. The cases shown here correspond to the binding of (a) three lipids and two waters with a Na+ ion, and (b) four lipids with a Na+ ion.

As we can see, the majority of Na+ ions bind to two or three lipids, whereas most of the lipids bind to one or two Na+ ions with almost equal probabilities. It is observed that, although each lipid has two binding sites for ion bridging, the cases with an individual Na+ ion binding to two carbonyl oxygens within the same lipid are rare.

As seen above, an ion often binds to more than one lipid. That makes it possible for the lipids to form larger connected clusters in which the lipids are linked to each other by ion bonds. We define such a continuous and connected network as an ion-lipid cluster. Figure 12a shows a typical configuration of an ion-lipid cluster. We also studied the distribution of the cluster sizes as shown in the inset in Fig. 13. Although smaller clusters appear frequently, most of the lipids belong to larger clusters. As seen in the inset in Fig. 13, some very large clusters also appear. An example of such a cluster is shown in Figure 12b.

Display large version of this figure
Figure 12
A snapshot of a typical configuration (one leaflet) demonstrating the existence of clusters. Each lipid (red) is represented by the average position of its two carbonyl oxygens. An ion-lipid bond (light blue dotted line; an ion can be bound to more than one lipid) appears between a lipid and a Na+ ion (blue). b) A snapshot of a giant cluster consisting of 51 lipids, 32 Na+ ions, and 86 ion-lipid bonds.
Display large version of this figure
Figure 13
Time evolution of the number of clusters during the whole simulation. Inset: Size distribution of clusters. The cluster size is defined by the number of lipids it contains, e.g., a cluster of size 1 might or might not contain a Na+ ion, whereas clusters of size more than 1 always contain Na+ ions.

It can be easily imagined that the existence of clusters has an effect on the dynamics of the whole bilayer since the robust ion bridges force the lipids to move collectively. This suggestion is supported by the strong correlation between the evolution of the area per lipid and the number of clusters (Figure 2 and Figure 13). The slow process of the formation and breaking of ion-lipid bonds is likely to be the key factor that makes the equilibration of the POPG bilayer slow.



Water orientation

Ordering of water molecules in the vicinity of the water/membrane interface can be characterized by the time-averaged projection of the water dipole unit vector onto the bilayer normal ,

(4)
where z is the z-component of the center of mass of a water molecule, and ϕ is the angle which the water dipole vector makes with the bilayer normal.

The orientation of water dipoles as a function of the distance from the bilayer center is presented in Fig. 14. Due to the negative surface charge of the POPG bilayer, interfacial water is highly polarized and distinctly different from the POPC case. Whereas in POPC there is a single minimum in the interfacial region, the ordering of water within the POPG bilayer is characterized by a rather deep minimum around 2.4nm and a broad maximum around 1.9nm. The water orientation profile reflects the total charge density profile of the POPG/water/ion system shown in Fig. 15.

Display large version of this figure
Figure 14
The average cosine of the angle ϕ between the water dipoles and the outward bilayer normal.
Display large version of this figure
Figure 15
Charge density profiles due to POPG with Na+ ions, water, and the total charge density profile. The profiles are shown as functions of the distance from the bilayer center.

As Fig. 15 shows, ordering of water compensates the excess local charge carried by lipid and Na+ ions. This is in agreement with observations for an anionic POPS lipid bilayer 62.


Electrostatic potential

To a large extent, the electrostatic potential determines the permeability of ionic solutes through the lipid bilayer 63. Hence, to calculate the electrostatic potential across the bilayer, the average charge density profile was first computed, the bilayer center being set to z=0 for each simulation frame. Then, the electrostatic potential was determined by integrating the charge density twice with the initial condition V(z=0)=0 (see, e.g., Patra et al. 44).

Fig. 16 shows that the electrostatic potentials for POPG and POPC bilayers differ especially in the inner side of the interface region. The potential barrier in POPG is much lower, which means that it is easier for a positively charged ion to migrate across the interface zone. This difference is expected to arise from a number of factors related to the charge density profiles of a POPG bilayer (see Fig. 15). First for a POPG bilayer the net negative charge carried by the phosphate group, along with the positive charge carried by the nonuniformly distributed Na+ ions, resulting in a strong surface potential, play an important role. The dipole potential due to the orientation of ester carbonyls is another potential source. Furthermore, the oriented water also makes an important contribution.

Display large version of this figure
Figure 16
Electrostatic potential across the bilayer as a function of distance from the bilayer center.


Discussion

We have presented results from an extensive MD study of a negatively charged POPG lipid bilayer. A zwitterionic POPC bilayer has been used as a reference. Thus far, to our knowledge, previous computational studies of PG systems 6,14,23 have been rather limited in scope and focused on hydrogen bonding 23, mixtures of PG with POPE 6, or PG molecules in monolayers 14. Here, we have presented a systematic study of the structural properties of a POPG bilayer with special attention to the role of electrostatics in PG bilayers.

Our studies revealed that PG headgroups have a rather weak intramolecular hydrogen bonding ability. We observed hydrogen bonds mainly between the hydroxyl-1 O16-H16 groups and carbonyl. This is in agreement with infrared spectroscopy analysis 61. Existence of this hydrogen bonding was suggested to prevent the interactions between hydroxyl-1 O16-H16 groups and Na+ ions (see the section Interactions of Na+ ions with water/membrane interface). The fact that the PG headgroup is not strongly hydrogen bonding is emphasized by the observation that in a mixture with POPE 6, the number of intramolecular bonds is three times smaller than what we have observed here. That is mainly due to competition for hydrogen bonding sites with the strongly hydrogen bonding PE group 6.

It seems that weakening hydrogen bonding of the PG headgroup is due to the strong direct electrostatic interactions caused by the negative charge within the POPG headgroup. A similar conclusion has been suggested for anionic POPS bilayers 13. On the other hand, the hydroxyl-2 (O15-H15) group is involved with high affinity in intramolecular hydrogen bonds with the phosphate oxygen O12. The existence of the intramolecular hydrogen bond O15-H15-O12, however, has not been reported in experimental studies thus far (to our knowledge). Infrared spectroscopy studies would be useful to investigate the existence and stability of this hydrogen bond.

The main finding in our study is that although the intermolecular hydrogen bonds are rare, the strong and frequent interaction between ions and PG carbonyl oxygens seems to be the most important factor determining the order of hydrocarbon chains and membrane surface area. About 60% of Na+ ions are involved in interlipid interactions in terms of ion bridges between the lipids. The resulting interaction network where lipids are connected through ion bridges results in ion-lipid clusters whose sizes range from a few to ∼100 particles. These interlipid interactions are stable, some of them persisting over the whole simulation run. Consequently, the ion-lipid clusters are also rather stable, implying that the lipids involved in a cluster form entities that act in a concerted fashion. It is likely that this plays a role in the structural as well as the dynamical properties of the bilayer, including, e.g., the movement of lipids in the plane of the bilayer and the possible formation of small-scale domains of PG-ion rich regions in many component membranes. Further, the strength and stability of positively charged ions and negatively charged headgroups provide a possible molecular level explanation for the action of PG in bacterial membranes under the influence of organic solvents. These highly polar interactions create a high free energy barrier for hydrophobic molecules such as benzene, thus possibly giving one reason nature uses charged lipids in these structures that have to adjust quickly to external changes.

PG headgroups also strongly modify the orientation of water near the bilayer. Two peaks of opposite sign in the interfacial region of the POPG bilayer related to the local water dipole order were observed (Fig. 14). Correlation of the position of these peaks with the position of the two peaks of opposite sign on the charge density (Fig. 15) along the bilayer normal suggests that the ordering of water molecules originates from the shielding of local charges by water dipoles. The snapshot of a Na+ ion coordinated with PG oxygen atoms and water molecules, which is additionally hydrogen bonded with a phosphate oxygen (Figure 11a), illustrates well the origin of the first peak. The second peak is due to water molecules hydrogen bonded with phosphate and glycerol oxygens from the side of the water phase. In the case of PC, water is ordered by the phosphate groups in a manner similar to PG, but the positively charged choline group does not act the same way as PG (which has no positive charge around) and Na+ ions around the choline group create a clathrates-like structure, thus the lack of ordering relative to the bilayer normal. The long-range polarization of the water molecules throughout the simulation box suggests that the amount of interlayer water in our model is too low for comparisons with unilamellar vesicles. That means that there is an additional repulsive interaction through the water phase between headgroups of opposite layers which potentially can influence the structure and dynamics of simulated bilayers. This is in agreement with experimental studies of charged lipids, which show continuous swelling with water added between lamellae up to a certain threshold at which two phases are formed, a fully hydrated unilamellar vesicle with bulk water 64. However, the exact amount of interlamellar water for POPG is not known. Despite that, we can exclude the possibility that the increase of the order and decrease of the surface area observed in our studies are additionally facilitated by a too small number of interbilayer water molecules. A similar increase for the order of PG acyl chains was observed in mixed PE-PG bilayers 6 and PG monolayers 14, and thus it is unlikely that it can influence our conclusions.


Acknowledgments

We thank the Finnish IT Center for Science (CSC) and the HorseShoe (DCSC) supercluster computing facility at the University of Southern Denmark for computer resources.

This work has been supported by the Emil Aaltonen foundation, the Academy of Finland through its Center of Excellence Program (I.V.), and the Academy of Finland projects.

References

1. Mead, J.F., Alfin-Slater, R.B., Howton, D.R., and Popják, G. (1986). Lipids—Chemistry, Biochemistry and Nutrition. (New York: Plenum Press). PubMed

2. Dowhan, W. (1996). Molecular basis for membrane phospholipid diversity: why are there so many lipids?. Annu. Rev. Biochem. 66, 199–232. CrossRef | PubMed

3. Pasenkiewicz-Gierula, M., Murzyn, K., Róg, T., and Czaplewski, C. (2000). Molecular dynamics simulation studies of lipid bilayer systems. Acta Biochim. Pol. 47, 601–611. PubMed

4. Hallock, K.J., Lee, D.-K., Omnaas, J., Mosberg, H.I., and Ramamoorthy, A. (2002). Membrane composition determines pardaxin’s mechanism of lipid bilayer disruption. Biophys. J. 83, 1004–1013. Abstract | Full Text | PDF (176 kb) | PubMed

5. Lee, S.-Y., and MacKinnon, R. (2004). A membrane-access mechanism of ion channel inhibition by voltage sensor toxins from spider venom. Nature 430, 232–235. CrossRef | PubMed

6. Murzyn, K., Róg, T., and Pasenkiewicz-Gierula, M. (2005). Phosphatidylethanolamine-phosphatidylglycerol bilayer as a model of the inner bacterial membrane. Biophys. J. 88, 1091–1103. Abstract | Full Text | PDF (398 kb) | CrossRef | PubMed

7. Uran, S., Larsen, A., Jacobsen, P.B., and Skotland, T. (2001). Analysis of phospholipid species in human blood using normal-phase liquid chromatography coupled with electrospray ionization ion-trap tandem mass spectrometry. J. Chrom. B 758, 265–275. PubMed

8. Tari, A., and Huang, L. (1989). Structure and function relationship of phosphatidylglycerol in the stabilization of phosphatidylethanolamine bilayer. Biochemistry 28, 7708–7712. PubMed

9. Maruyama, K. (2002). Peg-immunoliposomes. Biosci. Rep. 22, 251–266. CrossRef | PubMed

10. Isken, S., and de Bont, J.A.M. (1998). Bacteria tolerant to organic solvents. Extremophiles 2, 229–238. CrossRef | PubMed

11. Weber, F.J., and de Bont, J.A.M. (1996). Adaptation mechanisms of microorganisms to the toxic effects of organic solvents on membranes. Biochim. Biophys. Acta 1286, 225–245. PubMed

12. Pandit, S.A., and Berkowitz, M.L. (2002). Molecular dynamics simulation of dipalmitoylphosphatidylserine bilayer with Na+ counterions. Biophys. J. 82, 1818–1827. Abstract | Full Text | PDF (779 kb) | PubMed

13. Mukhopadhyay, P., Monticelli, L., and Tieleman, D.P. (2004). Molecular dynamics simulation of a palmitoyl-oleoyl phosphatidylserine bilayer with Na+ counterions and NaCl. Biophys. J. 86, 1601–1609. Abstract | Full Text | PDF (316 kb) | PubMed

14. Kaznessis, Y.N., Kim, S., and Larson, R.G. (2002). Simulations of zwitterionic and anionic phospholipid monolayers. Biophys. J. 82, 1731–1742. Abstract | Full Text | PDF (732 kb) | PubMed

15. Bandyopadhyay, S., Tarek, M., and Klein, M.L. (1999). Molecular dynamics study of a lipid-DNA complex. J. Phys. Chem. B 103, 10075–10080. PubMed

16. Gurtovenko, A.A., Patra, M., Karttunen, M., and Vattulainen, I. (2004). Cationic DMPC/DMTAP lipid bilayers: molecular dynamics study. Biophys. J. 86, 3461–3472. Abstract | Full Text | PDF (256 kb) | CrossRef | PubMed

17. Gurtovenko, A.A., Miettinen, M., Karttunen, M., and Vattulainen, I. (2005). Effect of monovalent salt on cationic lipid membranes as revealed by molecular dynamics simulations. J. Phys. Chem. B 109, 21126–21134. PubMed

18. Böckmann, R.A., Hac, A., Heimburg, T., and Grubmüller, H. (2003). Effect of sodium chloride on a lipid bilayer. Biophys. J. 85, 1647–1655. Abstract | Full Text | PDF (616 kb) | PubMed

19. Pandit, S.A., Bostick, D., and Berkowitz, M.L. (2003). Molecular dynamics simulation of a dipalmitoylphosphatidylcholine bilayer with NaCl. Biophys. J. 84, 3743–3750. Abstract | Full Text | PDF (186 kb) | PubMed

20. Böckmann, R.A., and Grubmüller, H. (2004). Multistep binding of divalent cations to phospholipid bilayers: a molecular dynamics study. Angew. Chem. Int. Edn. 43, 1021–1024. PubMed

21. Gurtovenko, A.A. (2005). Asymmetry of lipid bilayers induced by monovalent salt: atomistic molecular dynamics study. J. Chem. Phys. 122, 244902. CrossRef | PubMed

22. Gurtovenko, A.A., and Vattulainen, I. (2005). Pore formation coupled to ion transport through lipid membranes as induced by transmembrane ionic charge imbalance: atomistic molecular dynamics study. J. Am. Chem. Soc. 127, 17570–17571. CrossRef | PubMed

23. Elmore, D.E. (2006). Molecular dynamics simulation of a phosphatidylglycerol membrane. FEBS Lett. 580, 144–148. CrossRef | PubMed

24. Cowley, A.C., Fuller, N.L., Rand, R.P., and Parsegian, V.A. (1978). Measurement of repulsive forces between charged phospholipid bilayers. Biochemistry 17, 3163–3168. PubMed

25. Pozo-Navas, B., Raghunathan, V.A., Katsaras, J., Rappolt, M., Lohner, K., and Pabst, G. (2003). Discontinuous unbinding of lipid multibilayers. Phys. Rev. Lett. 91, 028101. CrossRef | PubMed

26. Pozo-Navas, B., Lohner, K., Deutsch, G., Sevcsik, E., Riske, K.A., Dimova, R., Garidel, P., and Pabst, G. (2004). Composition dependence of vesicle morphology and mixing properties in a bacterial model membrane system. Biochim. Biophys. Acta 1716, 40–48. PubMed

27. Wiedmann, T., Salmon, A., and Wong, V. (1993). Phase behavior of mixtures of DPPC and POPG. Biochim. Biophys. Acta 1167, 14–20. PubMed

28. Koynova, R. (1997). Liquid crystalline phase metastability of phosphatidylglycerols. Chem. Phys. Lipids 89, 67–73. CrossRef | PubMed

29. Wohlgemuth, R., Waespe-Sarcevic, N., and Seelig, J. (1980). Bilayers of phosphatidylglycerol. A deuterium and phosphorus nuclear magnetic resonance study of the head-group region. Biochemistry 19, 3315–3321. PubMed

30. Pascher, I., Lundmark, M., Nyholm, P.-G., and Sundell, S. (1992). Crystal structures of membrane lipids. Biochim. Biophys. Acta 1113, 339–373. PubMed

31. Garidel, P., and Blume, A. (2000). Miscibility of phosphatidylethanolamine-phosphatidylglycerol mixtures as a function of pH and acyl chain length. Eur. Biophys. J. 28, 629–638. CrossRef | PubMed

32. Lohner, K., Latal, A., Degovics, G., and Garidel, P. (2001). Packing characteristics of a model system mimicking cytoplasmic bacterial membranes. Chem. Phys. Lipids 111, 177–192. CrossRef | PubMed

33. Petrache, H.I., Tristram-Nagle, S., Gawrisch, K., Harries, D., Parsegian, V.A., and Nagle, J.F. (2004). Structure and fluctuations of charged phosphatidylserine bilayers in the absence of salt. Biophys. J. 86, 1574–1586. Abstract | Full Text | PDF (245 kb) | PubMed

34. Berger, O., Edholm, O., and Jahnig, F. (1997). Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature. Biophys. J. 72, 2002–2013. Abstract | | PubMed

35. Ryckaert, J.P., and Bellemans, A. (1975). Molecular dynamics of liquid n-butane near its boiling point. Chem. Phys. Lett. 30, 123–125. PubMed

36. Ryckaert, J.P., and Bellemans, A. (1975). Molecular dynamics of liquid alkanes. Faraday Discuss Chem Soc. 66, 95–106. PubMed

37. Reference deleted in proof..

38. Berendsen, H.J.C., Postma, J.P.M., van Gunsteren, W.F., and Hermans, J. (1981). Interaction models for water in relation to protein hydration. In Intermolecular Forces. Pullman, B., ed. (Dordrecht: Reidel), pp. 331–342. PubMed

39. Patra, M., and Karttunen, M. (2004). Systematic comparison of force fields for microscopic simulations of NaCl in aqueous solutions: diffusion, free energy of hydration and structural properties. J. Comp. Chem. 25, 678–689. PubMed

40. Hess, B., Bekker, H., Berendsen, H.J.C., and Fraaije, J.G.E.M. (1997). LINCS: a linear constraint solver for molecular simulations. J. Comp. Chem. 18, 1463–1472. PubMed

41. Miyamoto, S., and Kollman, P.A. (1992). SETTLE: an analytical version of the SHAKE and RATTLE algorithms for rigid water models. J. Comput. Chem. 13, 952–962. CrossRef | PubMed

42. Darden, T., York, D., and Pedersen, L. (1993). Particle mesh Ewald: an N log(N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089–10092. CrossRef | PubMed

43. Essman, U., Perela, L., Berkowitz, M.L., Darden, H.L.T., and Pedersen, L.G. (1995). A smooth particle mesh Ewald method. J. Chem. Phys. 103, 8577–8592. CrossRef | PubMed

44. Patra, M., Karttunen, M., Hyvönen, M.T., Lindqvist, P., Falck, E., and Vattulainen, I. (2003). Molecular dynamics simulations of lipid bilayers: major artifacts due to truncating electrostatic interaction. Biophys. J. 84, 3636–3645. Abstract | Full Text | PDF (167 kb) | PubMed

45. Anézo, C., de Vries, A.H., Höltje, H.-D., Tieleman, D.P., and Marrink, S.-J. (2003). Methodological issues in lipid bilayer simulations. J. Phys. Chem. B 107, 9424–9433. PubMed

46. Patra, M., Karttunen, M., Hyvönen, M.T., Falck, E., and Vattulainen, I. (2004). Lipid bilayers driven to a wrong lane in molecular dynamics simulations by truncation of long-range electrostatic interactions. J. Phys. Chem. B 108, 4485–4494. PubMed

47. Lindahl, E., Hess, B., and van der Spoel, D. (2001). GROMACS 3.0: a package for molecular simulation and trajectory analysis. J. Mol. Model. (Online) 7, 306–317.