INTRODUCTION

TOP ABSTRACT INTRODUCTION REFERENCES

Polyunsaturated fatty acid chains are an essential component of biomembranes. For instance, the retinal rod outer segment disk membrane is exceptionally rich in polyunsaturated fatty acids, with ~50% of docosahexaenoic fatty acid (DHA). Contents of DHA close to native levels are needed for the proper function of the visual receptor rhodopsin (Brown, 1994; Litman and Mitchell, 1996), the first three-dimensional crystal structure of which has recently been resolved at 2.8-Å resolution (Palczewski et al., 2000). In addition, polyunsaturated fatty acid chains are found in high concentrations in cerebral gray matter and synaptic plasma membranes and, from nutritional studies, it is well known that polyunsaturated lipids are important in the development of human brain function (Bloom, 1998).

Polyunsaturated phospholipid bilayers are characterized by low temperatures for the main (gel-to-liquid-crystalline) phase transition, T_m (Lipowsky and Sackmann, 1995). This allows fluidity of model membranes under physiological conditions in contrast to lipids with saturated chains of similar lengths (Small, 1986; Barry et al., 1991). In biological membranes, this fluidity is regulated by varying its composition (lipid and cholesterol content, for instance) since cells are usually constrained to an environment where temperature and pressure are fixed. For polyunsaturated lipids or even phospholipids with mixed saturated/polyunsaturated chains, the fluid lamellar phase can be achieved at room temperature for a unicomponent lipid bilayer. In the case of the 1-stearoyl-2-docosahexaenoyl-sn-glycero-3-phosphocholine (SDPC, 18:0/22:6 PC) lipid, in particular, the main order-disorder (L_'L) transition temperature in multilamellar dispersions containing 50 wt % H₂O was measured by NMR spectroscopy and found to be T_m = 7.7°C ± 0.4 (5.3°C ± 0.7) on decreasing (increasing) the temperature (Barry et al., 1991).

Using computer modeling two linear conformations for DHA were predicted (Applegate and Glomset, 1986): helical and an angle-iron-shaped form. Although the lower energy conformations correspond to a hairpin shaped molecule, the former were thought to be more suitable for biological interests, in which DHA can be paired with stearoyl in typical biological membranes. That study suggested that this kind of mixed-chain lipids with a saturated chain in the sn1 position and the DHA in the sn2 position may form chain arrays with relatively tight packing in certain conditions. The importance of polyunsaturated lipids in the local properties of membranes has recently been studied experimentally and specific lipid-protein and lipid-lipid interactions have been interpreted in terms of models with lateral segregation and the formation of domains (Litman and Mitchell, 1996). The aim of the present work is thus to study at a microscopic (atomistic) level these type of complex systems under physiological conditions. In this way, we are able to probe for the first time the molecular origin of the peculiarities that the system confers to the membranes of which it is a component.

In this paper we report the results of a molecular dynamics (MD) study of an SDPC lipid bilayer in the liquid crystalline phase at room temperature and ambient pressure. The MD simulations include structural properties of the lipid bilayer as well as conformational properties of the lipid molecules. To our knowledge, earlier studies of water-lipid phosphatidylcholine (PC) systems in the fluid lamellar phase were restricted to disaturated lipids (Merz and Roux, 1996) or lipids with a low degree of unsaturation. Those systems were constituted mainly by dimonounsaturated lipids, such as dioleoyl phosphatidylcholine (DOPC, 18:1/18:1) (Feller et al., 1997), or mixed saturated/unsaturated chains with one double bond, such as palmitoyl oleoyl phosphatidylcholine (POPC, 16:0/18:1) (Chiu et al., 1999), or two double bonds, such as palmitoyl linoleoyl phosphatidylcholine (PLPC, 16:0/18:2) (Hyvönen et al., 1997). In the case of systems with unsaturations, the maximum number of double bonds in a chain was two, and in the case of mixed chains, the simulations were not performed using a flexible cell geometry, which has been proved to be essential, at least during equilibration (Tobias et al., 1997; Venable et al., 2000). In all cases, except for simulations of dimyristoyl phosphatidylcholine (DMPC) and DOPC, the simulated temperatures were >50°C.

Computer simulation details

The initial configuration for the computer simulation was constituted by a monolayer of 32 lipid molecules arranged in the xy plane, which transformed into a bilayer after a rotation of 90° about the z axis and of 90° about one of the axis in the xy plane at the center of the bilayer. The bilayer system was then solvated by a water slab and the water molecules overlapped with the headgroups and those located deeper than the carbonyl region were removed. The initial cell dimensions were chosen to give an area per lipid, A, of A  69 Å²/lipid, and a lamellar spacing, d, of d  66 Å. These values were chosen to agree with the NMR and x-ray experimental data at room temperature (Koenig et al., 1997). The initial lipid conformation was the following: 1) the headgroups were pointing away from the acyl chain region toward the water-rich zone, 2) the saturated chain was built as all-trans, the minimum energy conformation at zero temperature, 3) the unsaturated chain consisted of an angle-iron structure for the double-bond region with the double bonds adopting a cis conformation, which is one of the conformations predicted by computer modeling for DHA in membrane environments (Applegate and Glomset, 1986), and 4) the molecular axes with the smallest inertia moments of the lipids were oriented randomly in the xy plane. The system with the previous characteristics was constituted by 64 lipid molecules and 27.5 water molecules per lipid (n_w), which corresponds to a fully hydrated lipid bilayer.

To be able to perform a quantitative analysis of the internal structure of the acyl chains, we have considered the following definitions for the different conformations associated with the dihedral angles (_k and _k) related to the torsional motions about the different bonds. For the saturated chain, the trans conformation is defined by the range 120°  _k < 240°, whereas the gauche conformations are defined by the ranges 0°  _k < 120° (gauche⁺) and 240°  _k < 360° (gauche). For the polyunsaturated chain, the skew⁺ (skew) conformations of the dihedral angles are given by the ranges 60°  _k < 180° (180°  _k < 300°) and the cis conformation corresponds to dihedral angles in the interval 60°  _k < 60°. For the latter chain type, we have considered the usual definitions for the dihedral angles between consecutive cis double bonds, which correspond to the structures with the following different shapes (Applegate and Glomset, 1986): helical (==skew^±skew^±== skew^±skew^±==; where "==" indicates the position of the double bonds), angle-iron (==skew^±skew^±==skewskew ==), and hairpin (==skew^±skew==skewskew^±== ).

The molecular and potential model used for the lipid molecules was the recent version of the all-atom CHARMM force field (CHARMM27) for lipids (Schlenkrich et al., 1996; Feller and MacKerell, 2000), which has been shown to give excellent results for disaturated lipids and for lipids with a low degree of unsaturation (Feller et al., 1997; Hyvönen et al., 1997). All the motions involving hydrogen atoms were frozen since those degrees of freedom are not expected to be relevant for the properties analyzed and this allowed us to use a longer time step in the integration of the classical equations of motion. We used a rigid TIP3P model for water (Jorgensen et al., 1983), which is consistent with the force field chosen for the lipid. The intermolecular parts of the force fields are pairwise additive functions, which consist of simple Lennard-Jones plus Coulomb terms.

The MD simulations were performed at constant temperature, T = 303 K, and pressure, p = 1 atm, and consisted of an equilibration period of 1.8 ns and an equilibrium run of 1 ns. During the equilibration period, the system was initially simulated at constant volume for 500 ps. Once the system was equilibrated at constant volume, it needed ~1300 ps to get a good convergence of the energy and the cell dimensions at constant pressure under NPT conditions using the Nosé-Hoover thermostat chain extended system isothermal-isobaric dynamics method with an orthorhombic simulation cell, as implemented in the PINYD computational package (Tuckerman et al., 2000). The reversible multiple time step algorithm (Martyna et al., 1996) permitted the use of a longer time step (5 fs) while the smallest time step was 1 fs. After equilibration, different properties were evaluated over a production run of 1 ns.

We used periodic boundary conditions and the constraints were handled by means of the SHAKE/ROLL and RATTLE/ROLL methods (Martyna et al., 1996). The short-range forces were computed using a cutoff of ~10 Å, and the minimum image convention, and the long-range forces were taken into account by means of the particle mesh Ewald summation technique (Frenkel and Smit, 1996).

Convergence and membrane dimensions

Due to the complexity (high degree of unsaturation) of the studied lipid, a long equilibration period (1.8 ns) was needed to achieve a convergence of the dimensions of the system. After the equilibrium was reached, the dimensions of the lipid bilayer evaluated over an additional 1 ns were A = 61.4 ± 0.5 Å² and d = 70.0 ± 0.5 Å, for the area per lipid and lamellar spacing, respectively. These data correspond to deviations of 10% (area per lipid) and 5% (lamellar spacing) from the initial set-up values.

In Fig. 1, we show the configuration of the system at the beginning of the equilibrium period. Although the initial values were chosen to correspond to the experimental data in excess water at the same temperature (Koenig et al., 1997), the combined NMR and x-ray experiments were performed on liposomes. In this kind of systems, previous studies on PC lipids in the L phase suggest that only the first 15 water molecules per lipid are incorporated homogeneously. For higher water contents, its addition causes the formation regions with water pockets apart from the well defined multibilayer lattice (Koenig et al., 1997). Therefore, one expects that, for n_w > 15, the reported experimental value for the lamellar spacing constitutes only a lower bound for our system. In MD simulations, however, upon the addition of water, this is incorporated in the interbilayer region homogeneously while the L phase is stable. Concerning the area per lipid, as discussed in Koenig et al. (1997), the value obtained by x-ray experiments (A^Xray) is predicted to be probably too large for n_w > 15 and, consequently, the area per lipid obtained by NMR measurements, A^NMR, may also be too large because A^NMR was found to be always > A^Xray in the range explored. The experimental results for the area per lipid are thus reliable only for low water content (Koenig et al., 1997), constituting an upper bound for our system. The evolution of the lipid bilayer dimensions in our simulation is, in this sense, in good agreement with the above discussion.

View larger version (138K): [in this window] [in a new window]   FIGURE 1   Configuration of the model bilayer system after 2 ns at constant pressure (NPT ensemble). Only the molecules in the simulation cell are shown. Molecules are depicted by their atomic van der Waals radii, and, for the sake of clarity, the hydrogen atoms of the lipids are not displayed. The color code is as follows: blue, nitrogen; yellow, phosphorus; red, oxygen; white, hydrogen; gray, carbon atoms; and green, C==C atomic groups.

Concerning the number of water molecules per lipid, n_w, corresponding to the fully hydrated SDPC bilayer, it is expected to be higher than the corresponding one for DMPC (n_w = 25.7) because the SDPC membrane absorbs more water than the DMPC at medium water content (Koenig et al., 1997).

The area compressibility modulus of the tensionless model membrane was estimated from the area fluctuations through the relation:
where k_B is the Boltzmann constant and A is the area per lipid, which gives K_A = 86.7 ± 0.7 dyn/cm. This is in quite good agreement with the experimental value of 121 dyn/cm (Koenig et al., 1997).

Orientational order

The orientational order has been studied through the orientational order parameter, S_CD, which can be directly measured experimentally by deuterium substitution NMR spectroscopy and is given by

(1)

where _n is the angle between the orientation of the vector along the CH bond of the nth carbon atom of the saturated (sn1) and/or the unsaturated (sn2) chains and the bilayer normal, and the brackets indicate averages over time and lipid molecules. The values adopted by S_CD are as follows: S_CD = 1 when the reference vectors are parallel, S_CD = 0 when the orientations are random, and S_CD = 0.5 when the reference vectors are perpendicular. In previous simulation studies, it was found already that there is a strong effect of unsaturation on the orientational order of the unsaturated chains (see for instance, Feller et al. (1997); Hyvönen et al. (1997)). Experimentally, it was observed that the polyunsaturated chains affect the order of the saturated chains (Salmon et al., 1987; Rajamoorthi and Brown, 1991; Holte et al., 1995). Specifically, the latter exhibits an increase in disorder at the bottom half of the chain as the degree of unsaturation is increased in PCs with mixed chains (Holte et al., 1995), whereas the effect of increasing temperature on the chain order consisted of a decrease in order in the first two-thirds of the saturated chain.

In Fig. 2, we plot the results obtained for S_CD as a function of the position, n, of the carbon atoms along the saturated and the polyunsaturated chains. As a general trend, the different molecular structure and dynamics of the polyunsaturated chain leads to values of the order parameters significantly lower than those of the saturated chain, in good agreement with experiments (Safley et al., 2000). The orientational order parameters for the DHA chain have been found to be strongly affected by unsaturation. This effect is correlated with the unsaturation position, especially at the first and second double bonds and in less degree at the third one. This correlation with unsaturation position suggests a dependence of the order parameters on the specific conformation of the polyunsaturated chain (Saiz and Klein, 2001). A decrease in order occurs at the position of the first double bond (C4==C5) and S_CD takes very small values between the first two bonds. The effect is just the opposite at the second cis double bond (C7==C8) where the order is increased, although S_CD takes smaller values than at the beginning of the chain. After the third double bond (C10==C11), S_CD decays to .05 for positions deeper into the bilayer interior but the behavior is not visibly as strongly correlated with the cis double bond position as in the case of the first two double bonds. This can be related with the fact that from computer modeling it was shown how the molecular shape (and packing) of diacyl glycerols depends on the position of the unsaturations (Applegate and Glomset, 1991a,b). For instance, stearoyl acid when paired with DHA in the sn2 position, the situation considered in this study, was classified in the same group as arachidonyl and eicosatrienoyl acids, which have at least three double bonds starting at the beginning of the chain. Diacylglycerols with two or three unsaturations starting at the middle of the chains pertained to the same type as distearoyl, and molecules with only one cis double bond at the middle of the chain showed a different behavior. Our results suggest that the addition of more cis double bonds after the first three when the unsaturation begins close to the headgroup does not have a strong effect on the order of the polyunsaturated chain.

View larger version (26K): [in this window] [in a new window]   FIGURE 2   Orientational order parameter as a function of the position of the carbon atom along the chains, SCD(n). The results of the saturated (sn1) and polyunsaturated (sn2) chains are plotted separately. Error bars are given as the standard deviations.

A similar correlation of the chain order with the position of a characteristic atomic group in the lipid chains was previously found (Husslein et al., 1998). There, the authors reported an MD simulation of diphytanol phosphatidylcholine (DPhPC) where the behavior of S_CD(n) suggested a stepwise decrease in order, correlated with the positions of the methyl group substitutions.

Density profiles

The electron density profiles (EDPs) have been computed assuming a gaussian distribution located at the atomic positions with variance, , equal to  = 2^3/2_LJ, where _LJ is the range of the Lennard-Jones potential. Hence, _LJ corresponds to the width of the gaussian distribution at half-height. The prefactors for each gaussian correspond to the atomic numbers, n_i. The EDPs computed in this way are proportional to the density profiles measured along the bilayer normal obtained by x-ray scattering experiments at low angles (White and Wiener, 1992, 1996) and are given by

(2)

where V is the volume of the slab between z and z + dz.

The EDP has been calculated separately for the different atomic groups of the lipid molecules. Maintaining the decomposition of White and Wiener (1992, 1996), the SDPC molecule is constituted by the following groups: choline group ([CH₃]₃NCH₂CH₂, denoted by CHOL.), phosphate group (PO₄), glycerol (C₃H₅, denoted by GLYC.), carbonyl groups (COO), alkene (HC==CH, denoted by C==C), methylene (CH₂) and methyl (CH₃) groups. It is worth noting that the C==C groups are located only in the unsaturated chain, whereas methyl and methylene groups are located in both chains. The component due to the carbonyl groups has been evaluated for the saturated (C⁽¹⁾OO) and the unsaturated (C⁽²⁾OO) chains. In Fig. 3 a, we show the results obtained for the (total) EDP, and the contributions arising from the water and lipid molecules separately. The components due to the acyl chains are depicted in Fig. 3 a as well. The contributions of the lipid headgroups, glycerol and the carbonyl groups to the EDP are decomposed in Fig. 3 b and compared to that of water.

View larger version (27K): [in this window] [in a new window]   FIGURE 3   Electron density profiles. (a) Total, components arising from the water and lipid molecules, and subcomponents coming from the hydrocarbon chains, namely, methyl, methylene and alkene groups; (b) components due to the lipid headgroups (phosphate, PO4, and choline, CHOL., groups), glycerol (GLYC.), and carbonyl groups (C(1)OO and C(2)OO for the saturated and polyunsaturated chains, respectively) and water molecules. The bilayer center is located at z = 0 Å.

The total EDP has the expected characteristics: 1) a higher density at the position of the lipid-water interface, corresponding to the headgroup, glycerol, and carbonyl groups of the lipids and water distributions; 2) a lower density in the bulk water and hydrocarbon region; and 3) a (slight) depletion at the center of the bilayer. The differences found between the two monolayers give an idea of the errors present in the simulations, whereas the asymmetry, splitting of the headgroup distributions shown in Fig. 3 b, is a consequence of shape deformations, fluctuations in curvature, of the model membrane. The decomposition of the total EDP into its components allow us to ascribe the distinct features of the curves to the different quasimolecular groups. Therefore, the first shoulder in the total distribution close to the bilayer center is due to the carbonyl groups and their covalently bonded ester oxygens, and the maximum corresponds to the phosphate and choline groups. It is interesting to note that the water distribution is constrained between the glycerol distributions of opposite monolayers, which is in excellent accord with the diffraction results for a DOPC lipid bilayer at low hydration (White and Wiener, 1992). Furthermore, we found a remarkable overlap between the water distribution and the double-bond region of the lipid chains, and between C==C distributions of lipids at the two opposite leaflets of the bilayer. These two features suggested that the double bonds might ferry water across the bilayer explaining, in this way, the enhanced permeability of membranes upon increase of unsaturation (White and Wiener, 1992). This mechanism was already confronted with simulations of DOPC (Feller et al., 1997) and PLPC (Hyvönen et al., 1997) bilayers. In the present study, however, due to the higher hydration and higher degree of unsaturation of the system, the overlap of water and C==C curves is visibly enhanced and the C==C distributions continuously expand the membrane interior. The favorable interactions between water molecules and the hydrocarbon region of the lipid will be considered in detail below.

The localization of the different quasiatomic groups at both sides of the membrane provides an approximate measure of the bilayer thickness. Nevertheless, it is difficult to give a unique definition for this quantity due to the complexity of the lipid bilayer, which can be divided into lipid-water interface and hydrocarbon region. The headgroup region lies at the interface and constitutes the most polar zone of the system with effective charges of +1e and 1e for the choline and phosphate groups, respectively. In this broad region, water is polarized and its effects counteract those of the lipid headgroup dipoles (Saiz and Klein, 2000). So far, several definitions have been proposed for the bilayer thickness that are generally accepted, namely, the Luzzati thickness, d_l, the Small thickness, d_L, and the distance between the maximums in the EDP, d_pp. In the case of the SDPC lipid bilayer, we obtained d_l = d_Chol.Chol.  45.6 Å, d_L = d_GLYC.GLYC.37.6 Å, and d_pp  42.7 Å, (basically, d_PO4PO4  43.3 Å). The distance between the choline and the glycerol groups, actually (d_L  d_l)/2, gives an estimate of the width of the water-lipid interface, which in our case amounts for d_interface  4 Å. A minimum of 10 Å of bulk water is located between the two interfaces of the SDPC lipid bilayer. Hence, our system is constituted by a low polar region in the core of the bilayer of 37.6 Å, two interfaces of 4 Å, each, and a bulk water slab of 24.4 Å.

The fact that the bilayer thickness is (slightly) smaller than the dimensions of the simulation cell in the bilayer plane indicates that the shape fluctuations of the membrane we observed represent bending deformations (fluctuations in curvature) since shape fluctuations represent bending modes of the model membrane for characteristic lengths somewhat larger than the bilayer system (Goetz et al., 1999). The bending rigidity of the system, , can be estimated from the membrane thickness, l_me, and the area compressibility modulus, K_A, through the relationship:  = K_Al/48, obtained by classical elasticity theory for thin solidlike films (Goetz et al., 1999). Although this relationship give excellent results for disaturated lipids, in the case of lipids with unsaturations it seems to predict a value for the /K_A ratio somewhat higher than the experiments (Rawicz et al., 2000). For l_me = d_L  37.6 Å, which corresponds to the deformable hydrocarbon region, we obtained  = (0.25 ± .01) × 10¹⁹ J which represents ~6 times the thermal energy, k_BT  4.18 × 10²¹ J. This result is consistent with experimental data for similar systems, even though, due to the dimensions of the system, we are just at the onset of the bending modes (collective undulations). In a recent simulation, however, an extensive study and a spectral decomposition of the mesoscopic undulations and thickness fluctuations modes was performed by simplifying the model and extending temporal and spatial scales (Lindahl and Edholm, 2000).

Hydrocarbon region and water-hydrocarbon chain interaction

To get more insight into the interactions of water with the hydrocarbon region of the lipid, we computed the number density profile along the bilayer normal for the acyl-chain carbon atoms as a function of their position in the chain. In Fig. 4, we plot the results obtained for the carbon atoms of the saturated and polyunsaturated chains and for the water molecules. There, the water number density was augmented conveniently for the sake of clarity. Our results show that carbon atoms at the end of the chains and, especially, those participating in the double bonds (C4==C5, C7==C8, C10==C11, C13==C14, C16==C17, and C19==C20) can reach the lipid-water interface. Furthermore, it is remarkable the disorder (accessible space and mobility) of the unsaturated chain atoms, which present quite broad distributions, especially when saturated and polyunsaturated chains are compared. It is worth realizing that the saturated-chain distributions (Fig. 4, top) are not specific of the studied system but display the common features of saturated chains in lipid bilayers. Note also that both chains show carbon atoms located at distances deeper than the bilayer center, manifesting an interpenetration of the hydrocarbon region of the monolayers. The broad distributions of the polyunsaturated chain carbon atoms, especially after the first cis double bond, and the accessibility of the interface by the C==C atomic groups is in good agreement with Cantor's theoretical (lattice statistical thermodynamics) calculations on the effect of unsaturations on the lateral pressure profiles in lipid bilayers (Cantor, 1999). Cantor found that the addition of unsaturated chains redistributes pressures from a broad region close to the center of the lipid bilayer to a region centered at ~45 Å below the lipid-water interface. For polyunsaturated lipids, this effect was found to be more pronounced. Cantor reported also a correlation with the position of the first cis double bond and, thus, the effect on the pressure profiles was more marked when the polyunsaturation began closer to the headgroup and continued deeper into the bilayer. The situation manifested in Fig. 4 (top) and (bottom) indicates as well the molecular origin of the low area compressibility modulus (high projected area fluctuations) of SDPC membranes (Koenig et al., 1997; Rawicz et al., 2000) when compared with disaturated lipid bilayers.

View larger version (58K): [in this window] [in a new window]   FIGURE 4   Number density of the acyl-chain carbon atoms along the bilayer normal for the saturated (top) and polyunsaturated (bottom) chains and of the water molecules (bottom). The distribution corresponding to water has been augmented for the sake of clarity. Bold curves stand for atoms located at the beginning and at the end of the chains and the bilayer center is located at z = 0 Å.

The water permeability of lipid membranes is known to be affected by the introduction of cis double bonds in the lipids. For instance, the apparent coefficient for water permeability at 21°C measured recently by micropipette aspiration (Olbrich et al., 2000) varies from ~30 to 40 × 10⁶ nm/ns for mono- and dimonounsaturated PCs, whereas, with two or more cis double bonds in the chain, the apparent permeability rises to ~50 × 10⁶ nm/ns for C18:0/2, to ~90 × 10⁶ nm/ns for diC18:2, and to ~150 × 10⁶ nm/ns for diC18:3. Those results are in qualitative agreement with previous measurements (see for instance, Huster et al. (1997), in which the authors reported a value of 239 ± 67 × 10⁶ nm/ns for the water permeability coefficient of SDPC at 25°C). The marked overlap between the water distribution and the double bond region of the lipid chains observed in the EDPs suggested a mechanism for the enhanced permeability of membranes upon increase of unsaturation (White and Wiener, 1992). However, it is difficult to predict favorable interactions between water and the C==C groups only from this fact and also to quantitatively compare the behavior observed on computer simulations of different lipid bilayers (Feller et al., 1997), which usually were performed at different conditions (temperature, hydration, membrane dimensions, etc.).

The use of a mixed chain (saturated/polyunsaturated) lipid allows us to probe the interactions of the two different chains with water under the same circumstances. Hence, we have evaluated the three-dimensional radial distribution functions (RDFs) of the oxygen atoms of the water molecules around the carbon atoms of the methyl (C^(0.27e)  H; with the effective charges in units of e, the electron charge, indicated in parenthesis), methylene (C^(0.18e)  H) and alkene (C^(0.15e) = H^(+0.15e)) groups of the polyunsaturated chain and for the carbons of the methylene groups of the saturated chain located at the same position as those considered for the polyunsaturated chain. The results obtained are plotted in Fig. 5. Due to the different length of the chains, the interactions of some atoms at the end of the polyunsaturated chain were not included in the calculations. The structure found in the g_XOw(r) (X = CH₃, CH₂, CH) functions for the polyunsaturated chain is an indication of preferential interactions of carbon atoms and the oxygen atoms of water. The first maximum (3.8 Å) corresponds to the distance between the carbon atom interacting with the water oxygen atom, whereas the second maximum (4.9 Å) corresponds to that between the carbon atom covalently bonded to a carbon atom interacting with the water O_w atom. In general, the probability of finding a water molecule around a carbon atom in the tails is clearly higher for the polyunsaturated chain than for the saturated one. For X = CH and CH₃, the curves corresponding to the two chains are qualitatively different. For the former, the peak in the g(r) for the polyunsaturated chain is indicative of favorable interactions with water, whereas for the saturated chain is structureless (when compared with that for X = CH₂ in the same chain). For the latter, although becomes nonzero at the same distance, indicating that some end atoms are reaching the water-lipid interface, there is not any affinity, whereas there is definitely some for the polyunsaturated chain. The similar qualitative results found for X = CH₂ are probably due to the major contribution of the atoms at the beginning of the chains.

View larger version (25K): [in this window] [in a new window]   FIGURE 5   Three-dimensional partial radial distribution functions, gXOw, of the oxygen atoms of the water molecules (Ow) around the carbon atoms of the X = CH, X = CH2, and X = CH3 groups for the polyunsaturated chain (sn2). The results are compared with those of the saturated chain (sn1) for those carbon atoms at the same position along the chain as the X = CH and X = CH2 groups for the polyunsaturated chain, and for the X = CH3 group.

Headgroup-headgroup and headgroup-water interactions

Lipid headgroups and water molecules are strongly organized at the interface (Saiz and Klein, 2000), which constitutes the most polar region of the membrane. This organization is not only found along the bilayer normal, where the water molecules are polarized to counteract the effect of the headgroup dipoles. At the plane of the interface, molecular dipoles and charges are also arranged forming a network of molecules interacting through electrostatic forces. In the present work, we have studied this organization at the interface of the SDPC membrane by evaluating the three-dimensional partial RDFs of the different pairs of polar components.

The choline and phosphate groups of the SDPC lipid are well hydrated as shown in Fig. 6 by the RDFs. The g_POw(r) function presents a sharp maximum at r = 3.75 Å, and the first minimum at r = 4.5 Å. The integration of the curve for the first hydration shell gives a coordination number of 6.2 water molecules around each phosphate group. The first maximum for the g_NOw(r) function (at r = 4.2 Å) is broader, and we obtained a mean number of 16 water molecules around the nitrogen atoms for distances r  5.75 Å. The RDFs (data not shown) for the non-ester oxygens, O_p, indicate that these atoms strongly interact with water. The g_OpOw(r) function has the first maximum at 2.55 Å, and the mean coordination number of water oxygens around O_p is 2.42 for r  3.25 Å. The presence of a small peak at similar distances (2.8 Å) for the ester oxygens of the headgroup gives a mean number of water oxygens of 0.6 for r  3.25 Å. These values (2 × 2.42 + 2 × 0.6) correspond to the 6 water oxygens around the phosphorus atom. Regarding the carbonyl oxygens (O_c), the g_OcOw(r) functions present a sharp first maximum at 2.75 Å, and a mean number of 0.85 to 0.9 water oxygens for r  3.25 Å, and of 0.175 for the ester oxygen atoms in the chains. This corresponds to ~1 water molecule around each carbonyl carbon. In summary, about 16, 6, and 2 water molecules are located around the choline, phosphate and carbonyl groups of the lipid, respectively. The analysis of the coordination number of the distinct oxygens around the choline group indicates that, although there are only 16 water oxygens around the nitrogen atom (for r  5.75 Å), its number increases up to 25.2 (16 + 9.2) when one considers all the oxygen atoms around a nitrogen atom including oxygen atoms in the same molecule (4.1) and neighboring lipid molecules (5.1). This corresponds to 6.3 oxygen atoms per methyl group around the tetramethylammonium group for the first minimum of the g_NOw(r) or 4.6 for the first minimum of the g_NOp(r).

View larger version (18K): [in this window] [in a new window]   FIGURE 6   Three-dimensional radial distribution functions of the oxygen atoms of the water molecules around the phosphorus (gPOw(r)) and nitrogen (gNOw(r)) atoms of the lipid headgroups.

The headgroup-headgroup interactions are evidenced by the structure of the RDFs for the nitrogen (N) and phosphorus (P) atoms shown in Fig. 7. There are very strong interactions between P and N atoms with a mean number of 1.2 N (P) atoms around each P (N) atom for distances r  6.2 Å, being 4.5 Å, the most probable distance for such a charge pair. It is worth remembering that the phosphate and choline groups have an effective charge of 1e and +1e, respectively. The previous result is in excellent agreement with former simulations of a DMPC lipid bilayer, in which charge pairs were also found between the lipid headgroups (Pasenkiewicz-Gierula et al., 1999). Furthermore, our findings are thus compatible with a picture in which the interface plane is constituted by a network of chains, rings, or pairs of lipid molecules connected via PN···P interactions. This representation is supported by the results obtained for the g_PP(r) and g_NN(r) functions. These two RDFs present similar characteristics, namely, a pronounced first and second maximums located at similar positions. Note that the second maximum is as high as the first one for the N-N interactions. In the case of g_PP(r) (g_NN(r)), the positions of the first and second maximums correspond to r  6.4 Å, (6.45 Å) and r  8.7 Å, (8.4 Å), respectively. The integration of the P-P (N-N) curve gives a coordination number of 1.7 (1.35) for r  7.4 Å, which is identical to the P-N coordination number at the same position, and 4.1 (4.1) for r  10.1 Å. The second maximum is compatible with the expected distances between P-P (N-N) pairs when there is a PN···P charge pair between the molecules and the angle between the intramolecular PN vector and the bilayer normal is close to the most probable value for the SDPC membrane (Saiz and Klein, 2000). For these interactions, similar distances are expected for P and N atoms and the fact that the g_NN(r) function is smoother than the g_PP(r) one can be attributed to the higher mobility of the choline group. The first peak in the g_PP(r) function was identified previously (Pasenkiewicz-Gierula et al., 1999) with water bridges. It is worth noting the similarities between intermolecular and intramolecular RDFs for the most probable distances between P and N atoms (Fig. 7, inset). This finding gives further support, in this case from computer simulations at the molecular level, to the fact that similar properties are found for free and nonfree ions in lipid bilayers. Examples of this phenomenon can be found, for instance, in the screening of charges in DNA complexes or planar surfaces of PCs, which leads to forces that does not depend critically on whether the phosphate and counterion are bonded, as in PCs, or not, as in DNA/tetramethylammonium systems (Parsegian and Rand, 1995), or in the equal screening of DNA charges by neutral and polar lipids in DNA/charged-lipid-bilayer complexes (Bandyopadhyay et al., 2000).

View larger version (29K): [in this window] [in a new window]   FIGURE 7   Intermolecular three-dimensional radial distribution functions of phosphorus and nitrogen atoms of the lipid headgroups. The inset corresponds to the distribution of intramolecular distances of phosphorus (nitrogen) atoms around nitrogen (phosphorus) atoms.

Interestingly, these PN···P interactions can actually take place directly or via interactions of O_p oxygens and the choline group through water molecules simultaneously hydrogen bonded to both groups. This explains the double first peak in the g_PN(r) function. The shortest distances can be ascribed to direct interactions whereas the largest ones are due to interactions through water molecules. The latter is clearly shown in the three-dimensional intermolecular distribution of the oxygen atoms of the water molecules and nitrogen atoms of the lipids around the phosphate groups depicted in Fig. 8 by the water rings around the O_p atoms. Because similar distances are preferentially adopted by those atoms in the same molecule (Fig. 7, inset), intramolecular charge pairs are present, as well, and interact in a similar fashion, directly or/and through water molecules.

View larger version (53K): [in this window] [in a new window]   FIGURE 8   Three-dimensional average intermolecular density isosurfaces of the oxygen atoms of the water molecules and the nitrogen atoms of the lipid molecules around the phosphate group of the SDPC molecule.

Further indication of the different nature of the double peak in the g_PN(r) functions is confirmed by the analysis of the orientational correlations present at the interface between the headgroup dipole moments (basically, the PN vector). The orientational correlations in molecular liquids can be evaluated through the functions (Bohm et al., 1983)

(3)

where P_l is the lth Legendre polynomial and (r) is, in this case, the angle between the headgroup dipole moments of two molecules the  and  atoms of which are located at a distance r. We have calculated G_l for l = 1 and l = 2 and the results obtained for ,  = P, N are shown in Fig. 9. The most important features are, on the one hand, the indication of antiparallel headgroups dipoles for nearest neighbors (all G₁(r) exhibit negative first peaks) and, on the other hand, the lack of dipole-dipole correlation for longer separations (G₂(r) decays to zero). In addition, the headgroup dipole moments of neighboring molecules are antiparallel for r _PN 5 °A whereas for 5 Å  r_PN  6 Å, which corresponds to the position of the second distribution forming the first peak, the headgroup dipole moments are oriented parallelly, however, the probability of an antiparallel or parallel orientation is quite similar, giving in this way a positive G₂(r) curve and zero values for G₁(r). For the region with some structure after the first minimum in the g_PN(r) function, the correlations are reversed and molecules tend to orient their headgroup dipoles parallelly. The third peak in the g_PN(r) corresponds to the P···N distribution for a pair of molecules whose other P···N atoms form a charge pair, whereas the second peak is not associated with such chains (or pairs or rings) of charge pairs.

View larger version (21K): [in this window] [in a new window]   FIGURE 9   Orientational correlation functions Gl(r), for l = 1, 2, of the headgroup dipoles for molecules whose phosphorus (solid line) and nitrogen (dashed line) atoms are separated by a distance r and, similarly, for two molecules whose P(N) and N(P) atoms are separated by r (dotted line).

Dihedral angles of the saturated and polyunsaturated chains

Up to now, we have dealt with general properties of the membrane or with those properties which arise from lipid-lipid and/or water-lipid interactions. In this section and the next one, we will address the structure of the individual acyl chains of the lipid. To study the conformation of the lipid chains, we have computed the dihedral angle distributions for the saturated and the polyunsaturated chains and the results are plotted in Fig. 10. In order to get meaningful information about the structure of the lipid molecules from the simulations, an equilibrium state should have been reached. The symmetry of the curves depicted in Fig. 10 shows that such an equilibrium has been attained even for the polyunsaturated chain. The distributions for the dihedral angles close to the beginning of the chain, the mobility of which is rather low compared to the rest of the atoms, are also reasonably symmetrical. Moreover, this is especially remarkable because we started from a single configuration of the lipid.

View larger version (26K): [in this window] [in a new window]   FIGURE 10   Distribution of dihedral angles for torsional motions about the covalent bonds along the chains for both the saturated (sn1) and polyunsaturated (sn2) chains.

The distributions for the dihedral angles in the saturated (sn1) chain show the typical behavior expected for saturated chains in lipid bilayers. The curves exhibit the usual maximums at ±60° for the gauche^± conformations and at 180° for the trans conformation. The presence of some gauche defects in these chains is an indication of the disordered state of the saturated chains. In the case of the polyunsaturated (sn2) chain, the distributions of the dihedral angles corresponding to the positions of the double bonds have only one maximum at 0°, which corresponds to a cis conformation. For the dihedral angles located between two (cis) double bonds, the distributions present the two expected symmetric maximums at ±120°, which correspond to the skew^± conformations.

Polyunsaturated chain conformation

The starting configuration of the lipid was an angle-iron structure for the double-bond region of the polyunsaturated chain. This structure consists of the successive disposition of the  ==skew^±skew^±==skewskew== conformations. This initial conformation was adopted because it was one of the two (angle-iron and helical) predicted by a molecular modeling approach (Applegate and Glomset, 1986; Albrand et al., 1994) to be relevant for biological membranes, where DHA could be paired with a saturated chain. These two structures have almost straight chain axes, and Applegate and Glomset (1986, 1991b) showed that rather tight intermolecular packing arrangements were possible, especially for the angle-iron conformation case, for DHA chains alone or when DHA was paired with stearoyl acid in diacylglycerols.

In the present study, after the equilibrium was reached, we obtained that most of the pairs of dihedral angles (81%) between two consecutive cis double bonds adopted a ==skew^±skew^±== conformation, where the two consecutive double bonds are almost parallel, in contrast to the ==skew^± skew== one (19%), where the two consecutive double bonds are almost perpendicular. The percentage of pairs of dihedrals between two consecutive cis double bonds was computed also as a function of its position along the chain. The results obtained are summarized in Table 1. We found only a small variation of the two populations as a function of their position, which indicates that the results are reliable and that a meaningful equilibrium has been reached. For instance, the percentage of ==skew^±skew^±== angles decreases from an 85% at the beginning of the chain to a 78% at the chain end.

View this table: [in this window] [in a new window]   TABLE 1   Population of the ==skew± skew±== and == skew± skew ==conformations

The population of the helical (==skew^±skew^±==skew^± skew^±==), angle-iron, and hairpin (==skew^±skew== skew^±==) conformations were studied (among others) for the four dihedral angles located between three consecutive cis double bonds. Those populations are reported in Table 2 as a function of its position along the chain of the double-bond segments and the mean values are included, as well. The helical conformation was the other structure predicted by Applegate and Glomset (1986) to allow a relatively tight packing of the chains. The hairpin structure, however, is the conformation with the lowest energy in the gas-phase in which the chain is curled toward itself to maximize intramolecular interactions. We found that helical and angle-iron conformations are quite stable for the present thermodynamic state. On average, helical and angle-iron structures represent 37 and 29% of the groups of dihedral angles, respectively, which together constitute a 66% of the total population. Nevertheless, we found a significant fraction of dihedrals (34%) with conformations in which a tight packing of the chains is not possible. In particular, 2% of the dihedral angles adopt a hairpin conformation and 32% adopt other types of hairpin-like structures. It is also worth evaluating the fluctuations of these populations. In Fig. 11, we show the evolution of the previous quantities with time. The curves present oscillations around the mean values without great differences depending on the position of the groups of dihedral angles along the chain. Variations in the different populations indicate transitions between skew^±  skew conformations.

View larger version (52K): [in this window] [in a new window]   FIGURE 11   Evolution of the populations of the different conformations (helical, angle-iron and hairpin) with time as a function of the position of the three consecutive double bonds along the chain (n = 1 (upper/headgroup), 2, 3, and 4 (lower/chain-end)).

The study of the temporal evolution of the dihedral angles between two cis double bonds gives more insight into the stability of the conformations and the connection between chain conformation and the physical properties of the model membrane. Two main features can be inferred from the dependence of the dihedral angles with time. On the one hand, transitions between skew^± and skew take place in a nanosecond time scale. On the other hand, it is remarkable the correlation found between (==skew^±skew^±==)  (== skewskew==) transitions when the two dihedral angles between two consecutive double bonds have the same sign in contrast to the uncorrelated transitions of dihedral angles with opposite signs. Fig. 12 illustrates this situation for a pair of consecutive dihedral angles. This correlation is an indication of the stability of the dihedrals of the same sign (helical and angle-iron conformations) between double bonds, which is in excellent agreement with the predictions from molecular modeling (Applegate and Glomset, 1986, 1991a). It shows that when two consecutive double bonds are aligned, they continue to be so. When they are perpendicular, however, they are less correlated. This behavior is easily explained by the fact that transitions from dihedrals of the same sign to dihedrals of different sign require a structural reorganization with relatively large transversal (local) fluctuations compared to the situation where those transitions that take place from dihedrals of the same sign to dihedrals of the same sign.

View larger version (36K): [in this window] [in a new window]   FIGURE 12   Temporal evolution of a pair of dihedral angles between two consecutive cis double bonds.

The results obtained for the conformations and intramolecular dynamics of the polyunsaturated chains involve a broad distribution of projected area per chain, and quite large local fluctuations when a transition between a linear and a nonlinear conformation takes place. This is directly connected with the small area compressibility modulus, i.e., with large fluctuations of the projected area, of polyunsaturated lipids found theoretically (Cantor, 1999) and experimentally (Koenig et al., 1997). The opposite picture is likely, as well. The large fluctuations of membrane embedded proteins (bundles of transmembrane peptides, for example), such as the Metarhodopsin I  Metarhodopsin II transition of the visual protein rhodopsin, the light-sensitive photopigment of the rod cells of the vertebrate retina, may be rather easily accommodated by the membrane by inducing a reorganization of the lipid chains. Our results suggest that this exceptional system, the mixed polyunsaturated-saturated membrane, can be especially capable of adjusting such a changes in volume, i.e., area in the membrane plane and bilayer thickness along the membrane normal (Litman and Mitchell, 1996; Brown, 1994; Cantor, 1999). Actually, recent NMR experiments on SDPC model membranes found differences in the order parameters along the DHA chain in presence and in absence of the protein rhodopsin (Safley et al., 2000). These variations suggested a change in the average conformation of the polyunsaturated chain.

Concluding remarks

The structure of a fully hydrated mixed (saturated/polyunsaturated) chain phosphatidylcholine lipid bilayer in the biologically relevant liquid crystalline phase has been examined by performing a molecular dynamics study at constant (room) temperature and (ambient) pressure. In general, we have found a reasonable agreement between MD findings and the available experimental structural information under similar conditions.

The order parameters obtained for the docosahexaenoic fatty acid chain were lower than those for the saturated chain, in good agreement with experiment (Safley et al., 2000). The orientational order parameter, as a function of the carbon position in the polyunsaturated chain, has been found to be affected by unsaturation, as expected, especially at the first three cis double bonds. In this region, a sudden loss and gain of order occurs, which is correlated with the unsaturation position. This correlation indicates a dependence of the order parameters on the specific polyunsaturated chain conformation (Saiz and Klein, 2001). Our results suggest, as well, that the addition of more cis double bonds after the first three, when the unsaturations begin close to the headgroup, does not have an additional effect on the order of the polyunsaturated chain. This is compatible with previous studies on the effects of the degree and unsaturation position on the molecular shape of lipids (Applegate and Glomset, 1991a).

We found remarkably overlapping atomic distributions for end-of-chain carbons and, especially, for the double bond region of the polyunsaturated chain and the water molecules at the interface, in agreement with previous computer simulations and experiments (White and Wiener, 1992; Feller et al., 1997; Hyvönen et al., 1997). In the case of the SDPC/water system, we found that this overlapping of distributions is more pronounced. The present study of a mixed (saturated/polyunsaturated) chain lipid allowed us to analyze the differences between the water-hydrocarbon chain interactions for the polyunsaturated and saturated chains at the same conditions of temperature, pressure, hydration, membrane dimensions, etc.. Thus, we observed that the water molecules interact mainly with the polyunsaturated chains, which agrees with the enhancement of membrane permeability to water and small organic solvents with increasing unsaturation (Olbrich et al., 2000). The atomic distributions of the saturated chains are similar to those obtained previously for other saturated lipids. On the contrary, we found broad distributions for the polyunsaturated chain carbon atoms, especially after the first cis double bond. These two features, together with the accessibility of the interface by the C==C atomic groups can be connected with the effects of unsaturation on the lateral pressure profile in lipid bilayers (Cantor, 1999).

Lipid headgroups and water molecules have been observed to be strongly organized at the broad lipid-water interface. In this region, molecular dipoles and charges are arranged forming a network of molecules interacting through electrostatic forces. For instance, P···N interactions among headgroups (charge pairs), both directly and through water molecules are evident, with mean distances for the P···N charge pairs similar to those corresponding to the intramolecular interactions. Lipid headgroups thus form chains, rings, or pairs of strongly interacting molecules. The study of the orientational correlations between headgroup dipoles indicates antiparallel dipoles for nearest neighbors and a lack of dipole-dipole correlations for longer separations.

The helical and angle-iron conformations of the region of the polyunsaturated chains comprised between three consecutive cis double bonds are shown to be quite stable for the studied thermodynamic state. These conformations permit a relatively tight packing of the chains since consecutive cis double bonds are parallelly oriented. Nevertheless, we have found a significant fraction of molecules with conformations (hairpin and hairpin-like) in which such a tight packing is not possible. This leads to a high degree of inhomogeneity in the system. The results obtained for the conformations and intramolecular dynamics of the polyunsaturated chains involve a broad distribution of projected area per chain, and quite large local fluctuations when transitions between conformations where double bonds are parallel and those with perpendicular double bonds take place. There is a remarkable correlation found between (skew^±skew^±  skewskew) transitions when the two dihedrals have the same sign, in contrast to the transitions of dihedral angles with opposite signs. The small area compressibility modulus of polyunsaturated membranes (Cantor, 1999; Koenig et al., 1997) and capability of polyunsaturated rich membranes to accommodate the quite large changes in volume occurring in membrane embedded proteins (Litman and Mitchell, 1996; Brown, 1994; Cantor, 1999), such as the Metarhodopsin I  Metarhodopsin II transition of the visual protein rhodopsin, are thus explained in terms of lipid reorganization at a microscopic (molecular) level.