INTRODUCTION

TOP ABSTRACT INTRODUCTION MODEL AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

A biological membrane is an assembly of a diverse set of molecules such as lipids, sterols, and proteins. Among the different lipids found in biomembranes, phospholipids play an important role. Phospholipids can be charged, as in the case of phosphatidylserine (PS) molecules, or neutral but zwitterionic-like phosphatidylcholine (PC) and phosphatidylethanolamine (PE). The molecules in a membrane assemble into a bilayer that is nonsymmetrical with respect to charge distribution. Thus, for example, in the mammalian plasma membrane, the negatively charged PS molecules are mostly found in the cytosolic leaflet of the membrane where they interact with a variety of proteins involved in the signal transduction pathway (Langner and Kubica, 1999).

To simplify the complex problem presented by natural membranes, one studies model membranes, although even in model membranes the situation remains complicated. Thus, it was observed that, in a model membrane containing a mixture of PS and PC, the mixing is nonideal (Huang et al., 1993). In other model membranes containing PS, interactions between PS and proteins produce domain formation (Denisov et al., 1998).

It is often assumed that interactions between charged membranes and charged peptides are governed by electrostatic forces. Theoretical studies of such interactions frequently use continuum description of the system (Ben-Tal et al., 1996). Computer simulations may provide a more detailed molecular description of these interactions. Such simulations can also provide justification for the use of continuum models and provide the reason for its success or failure. Computer simulations on model membranes were initially performed on homogeneous membranes containing only one type of phospholipid molecules, usually of the zwitterionic type (Tobias et al., 1997). More recently, simulations on membranes that contain mixtures of phospholipids with sterols (Smondyrev and Berkowitz, 1999a) or phospholipids with surfactants (Schneider and Feller, 2001; Bandyopadhyay et al., 2001) appeared in the literature. We are interested in a simulation study of the electrostatic properties of a membrane that contains a mixture of charged and uncharged phospholipid molecules, for example, a membrane containing a mixture of PS and PC molecules. Before performing a simulation on such a mixture, we decided to perform a simulation on a pure PS bilayer. This paper presents results obtained from our simulation and comparison of these results with the results from a previous simulation study performed on a similar bilayer (López Cascales et al., 1996).

It was observed that many biologically important properties of dipalmitoylphosphatidylserine (DPPS) bilayers are dependent on the pH of the environment, nature, and concentration of counterions. Among these properties are the transition temperature (MacDonald et al., 1976), packing, headgroup hydration (Ermakov et al., 2001; Papahadjopoulos, 1968), bilayer fusion, and phase separation in mixtures containing DPPS (Jacobson and Papahadjopoulos, 1975). A comparison of the effect of monovalent counterions and bivalent counterions is a subject of our future work. Here, we concentrate on a study of a PS bilayer, because PS molecules play an important role in supporting homeostasis in the brain by actively participating in such processes as conduction of nerve impulses, and accumulation, storage, and release of neurotransmitters.

	MODEL AND METHODS

TOP ABSTRACT INTRODUCTION MODEL AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Force field parameters

Previous simulations on phospholipid membranes performed in our group were mostly done on phosphatidylcholine molecules. We used a united atom (UA) force field in those simulations, and we would like to continue using a UA force field in simulations that also involve phosphatidylserine molecules. Nevertheless, to treat the NH group of the PS headgroup, we decided to include explicit hydrogens into our description. This turned out to be useful for the analysis of the intermolecular hydrogen bonding. Because the structure of the PS molecule is similar to that of PC molecule (see Fig. 1), we, wherever possible, used the force field developed by Smondyrev and Berkowitz (1999b) for the simulation of PC bilayers. The force field parameters that are different from those used in the description of PC were derived in a way to be consistent with the parameters adopted for the PC.

View larger version (17K): [in this window] [in a new window]   FIGURE 1   The structure of PS and PC molecules.

To find the missing parameters for the PS molecule, we started with a calculation of the partial charges on the atoms of the headgroup of PS molecule. We assigned no charges to the UA of the tails and, therefore, excluded these hydrocarbon tails from our quantum calculations. The calculations were done using the quantum chemistry software, Gaussian 98 (Frisch et al., 1998) at the HF level with the 6-31+G* basis set. The initial input for the PS molecule was prepared using SYBYL 6.7 and minimized using the built-in energy minimization routine with the Tripos Force Field (SYBYL, St. Louis, MO). The geometry of the molecule was optimized with Gaussian 98 before charge computations. Finally, the charges were calculated using the CHELPG method for the 15 different conformations generated using Free Software Ghemical (Hassinen et al., 2000). Figure 2 shows the final partial charges. The reported charges were averaged over 15 different conformations. We observed that partial charges on the atoms common to both PC and PS are similar in value to each other (Smondyrev and Berkowitz, 1999b).

View larger version (23K): [in this window] [in a new window]   FIGURE 2   The partial charges used in the simulation (Hydrogen charges are summed into the heavy atoms except in the case of NH).

Because all the bond, angle, nonbonded parameters, and most of dihedral angle parameters were taken from the work of Smondyrev and Berkowitz (1999b), we needed to determine parameters for only three dihedrals: ₀  C2-CH-N3-H, ₁  C2-CH-C-O₂ and ₂  N3-CH-C2-OS. We decided to adopt the Amber parameter for ₀ and determine the parameters for the other two by performing calculations in the spirit of the work by Smondyrev and Berkowitz (1999b). Therefore, we chose small model compounds for which we calculated the torsional energy profiles. The following procedure was used in the computation of torsional parameters:

1.	We calculated the ab initio energy profile of a compound as a function of the dihedral angle of interest. For each dihedral angle, we optimized the respective model compound at the HF level with the 6-31+G* basis set.
2.	We took the optimized configurations of the model compounds from the ab initio calculations and computed single point energies using the Sander module of AMBER 6.0 (Case et al., 1999). In these calculations, we omitted the contribution from the term corresponding to the dihedral angle of interest (parameters for the other dihedrals were taken from the work of Smondyrev and Berkowitz (1999b)). We verified that the equilibrium angles and bond lengths were in agreement with the OPLS/AMBER (Jorgensen and Tirado-Rives, 1988) bond lengths and angles. For Sander calculations, we chose the value of the dielectric constant  = 1, cutoff of 10 Å and 1-4 screening parameters as 8.0 and 2.0 for van der Waals and electrostatic interactions, respectively.
3.	The difference between the two energy profiles from the ab initio and Sander calculations was fitted to the function so that the total energy profile from the force field calculation produced a reasonable fit to the profile obtained from ab initio calculations. Figure 3 shows the model compounds used to compute the torsional energy for 1, 2. The first compound we chose to compute potential parameters associated with the dihedral 1 was Serine. The partial charges used in these computations were the same as those obtained for the DPPS molecule (See Fig. 2). The charge on the OH group was adjusted so that the total charge on the compound was zero. Following the procedure described above, we obtained the parameters for 1 (See Table 1). Figure 4 a shows the match of fitted energy with the ab initio energy.

View larger version (14K): [in this window] [in a new window]   FIGURE 3   The model compounds used in the computation of force field parameters for dihedral angles.

View this table: [in this window] [in a new window]   TABLE 1   Parameters used in the simulation with PS that are different from the ones used in simulations with PC (Smondyrev and Berkowitz, 1999b)

View larger version (5K): [in this window] [in a new window]   FIGURE 4   (a) The energy profile for 1. (b)The energy profile for 2.

Considering that there is a possibility of hydrogen bonding between NH and PO in PS, we chose a larger model compound, serine methyl phosphate ion, for the computation of parameters for ₂. Again, the charge on methyl group was adjusted so that the total charge on the compound was 1 e.s.u. Parameters for torsional motion around ₁ obtained from calculations on serine were used in these calculations. Figure 4 b shows the force field energy and the ab initio energy. As the figure shows, the fit in this case is not as good as in the case of serine. It even indicates the wrong location of the global minimum. Nevertheless, quantum calculations show that the proton transfer from the ammonium to carbonyl group occurs around the global minimum at 80^°. This means that our dihedral potential should restrict the variation of the dihedral angle to the region around the other minimum located at 30^°. This was accomplished by the fitted potential. The distribution of the dihedral angle ₂ obtained from the simulation (not shown) confirms this. We list the parameters associated with dihedral angles ₀, ₁, and ₂ in Table 1.

The molecular dynamics simulation

The DPPS bilayer in our simulation consisted of 128 DPPS molecules, 64 molecules in each monolayer. The monolayers were formed by randomly rotating and copying a DPPS molecule 64 times. The overlap cutoff for the random placement was 2.1 Å. We placed the two monolayers on top of each other so that the phosporus-phosphorus distance was 50 Å with the headgroups pointing outside. Two slabs of TIP3P water, each with 1312 water molecules, were added above and below this bilayer. The Na⁺ counterions were randomly added to both slabs of water. The initial distance between the ions and the plane of the bilayer (the initial plane of the bilayer was considered to be located 3 Å away from the plane of phosphorus atoms) was between 5 and 10 Å. The distance between monolayers was gradually reduced to the phosphorus-phosporus distance of 38 Å, performing short molecular dynamics runs, each of 2-ps duration at 350 K. After every run the distance was reduced by 0.5 Å. The phosphorus atoms were kept frozen during this process. After adjusting the z dimension of the simulation cell to 63 Å, we performed a 20-ps simulation with free phosphorus atoms. To ensure the disorder in hydrocarbon tails, we raised the temperature to 430 K and then brought it down to 350 K by performing a set of 20-ps molecular dynamics runs, each at a temperature lowered by 20 K. Finally, we equilibrated the system for 50 ps at 350 K.

After equilibration, we performed a 4-ns simulation at constant pressure p = 1 atm and temperature 350 K. The main transition temperature for the DPPS bilayer with Na⁺ counterions is 329 K (MacDonald et al., 1976). Hence, the simulations were performed at 350 K to ensure that the bilayer was in the liquid crystalline state. Moreover, because the previous simulations (López Cascales et al., 1996) were also performed at 350 K, it was easier to perform comparisons between our simulation and previous ones. Periodic boundary conditions in all three dimensions were applied. The thermostat and barostat relaxation times were 0.2 and 0.5 ps, respectively. Smooth particle mesh Ewald (SPME) (Essmann et al., 1995) was used with 10⁴ tolerance for computation of long range electrostatic contributions. All the bond lengths in the simulation were held fixed using the SHAKE algorithm with tolerance 10⁴. The calculations were done at the North Carolina Supercomputing Center using the DL_POLY (Smith and Forester, 1999) simulation package version 2.12. We monitored the area per headgroup and observed that it was stable during the last 2 ns of the run (See Fig. 5). Therefore, we performed our analysis using the data obtained from the last 2 ns. The size of the simulation cell during the last 2 ns was ~58.7, ~58.7, and ~63.0 Å along x, y, and z directions, respectively. To make sure that the ions were in a stationary state during the run, we monitored the z coordinate of the distance between the ions center of mass and the center of the bilayer and observed that it was stable during the 2-ns run used for the analysis. We also plotted the distribution of ions as a function of z coordinate obtained from the first and second nanosecond of simulation (see Fig. 6). The figure shows the similarity in the distribution.

View larger version (31K): [in this window] [in a new window]   FIGURE 5   Area per headgroup versus time.

View larger version (10K): [in this window] [in a new window]   FIGURE 6   Distribution of counterions as function of the z coordinate during the first half and the second half of the simulation.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION MODEL AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Area per headgroup and atomic distribution

From our simulations, we observed that the average area per headgroup for DPPS molecules is 53.75 ± 0.10 Å². A similar value of 54 Å² per DPPS headgroup was obtained in the simulation performed by López Cascales et al. (1996). The values inferred from the experiments are in the region between 45 and 55 Å² (Cevc et al., 1981; Demel et al., 1987). Note that the average area per headgroup for DPPC bilayer obtained from simulations and measurements is ~62 Ų. One would expect that the area per headgroup in case of DPPS molecules will be larger than the area per DPPC molecule, because the headgroups of DPPS are charged and therefore repel each other. Contrary to this expectation, the area per headgroup in DPPS is ~13% smaller than that of DPPC. López Cascales et al. (1996) proposed that this reduction in the area per headgroup is due to a strong intermolecular coordination between DPPS molecules. We will return to this issue shortly.

To appreciate the dimensions of the DPPS membrane and sizes involved in the problem, we show, in Table 2, the average distances (from the bilayer center) of some headgroup atoms, and, in Fig. 7, the distributions of various atoms in the bilayer as a function the z coordinate. Absence of water in the region between 10 and 10 Å clearly indicates the hydrophobic nature of the hydrocarbon tails. Figure 7 b also shows that Na⁺ ions are found mostly next to the serine group, therefore substantially compensating the charge on the headgroup. The Na⁺ ion distribution peak is located around ~21 Å and it falls-off rapidly. This is contrary to the results obtained from a rather short-timed simulation by López Cascales et al. (1996). Figure 8 shows typical coordinations of Na⁺ around PS molecules.

View larger version (0K): [in this window] [in a new window]   FIGURE 7   Density of various atoms in the simulation as function of the z coordinate.

View larger version (81K): [in this window] [in a new window]   FIGURE 8   Snapshot pictures showing the coordination of Na+ around the DPPS molecule. In both pictures, we see that the ion is coordinated with the serine and phosphate oxygens.

We noticed that the headgroup distribution peaks are sharper for DPPS than for DPPC. Assuming the peaks to be gaussian, we found the full width at half maxima of the peaks to be ~3-4 Å (for phosphorus it is 3.0 ± 0.5). These widths are ~6.6 Å (López Cascales et al., 1996; Egberts et al., 1994) in DPPC bilayer. This tells us that the headgroups in DPPS are less mobile, at least during the 2-ns time of the simulation.

We also studied the orientation of the headgroups with respect to the bilayer normal. Figure 9 shows the distribution of the cosine of the angle between the P-N vector and outward normal to the bilayer. Comparison of this distribution with the one for DPPC (Smondyrev and Berkowitz, 2000) reveals that DPPS headgroups are slightly more inclined with respect to the bilayer plane. The average angle was found to be ~68^°.

View larger version (2K): [in this window] [in a new window]   FIGURE 9   Distribution of the cosine of the angle between the vector joining phosphorus and nitrogen in DPPS and the outward normal to the bilayer.

Area per headgroup and hydrogen bonding

Experimental measurements performed on bilayers containing PE molecules in their liquid crystalline state show that the area per headgroup for PE is ~10-20% (Thurmond et al., 1991) smaller than the area per PC molecule. For DLPE, it is 50.6 Å² (McIntosh and Simon, 1986; Damodaran et al., 1992). The estimated area per molecule for DPPE is ~55.4 Ų (Thurmond et al., 1991). We note that the PE headgroup is just the PS headgroup without the carboxyl group. So we suspect that the behavior of PS may be similar to that of PE.

As we already mentioned, the Na⁺ counterions are closely located to carboxyl groups. To examine the location of counterions in more detail, we considered the radial distribution functions (rdf) between the Na⁺ and various atoms of PS. The rdf is defined as
where N(r) is the number of Na⁺ in the shell between r and r + r around the PS atoms.  is the number density of Na⁺, taken as the ratio of the number of atoms to the volume of the simulation cell. Figure 10 shows two rdf: the first one for Na⁺ and carbon in the carboxyl group and the second one for Na⁺ and phosphorus from the phosphate group. The sodium ions were found to be in close proximity to both negatively charged groups, thus shielding the charge of the headgroup. We also observe that the ions are closer to carbon, confirming our previous assessment that Na⁺ compensate the carboxyl group charge. We notice that the rdf between Na⁺ and carbon has two peaks, indicating that there are two most probable locations for the counter ions around the carboxyl group. Because Na⁺ ions compensate most of the charge on serine oxygen, it appears that, in PS, a dipole is formed between NH and PO, i.e., PS is essentially behaving like a PE and, therefore, one can expect that the area per headgroup in PS bilayer will be more like the one observed in PE bilayer.

View larger version (6K): [in this window] [in a new window]   FIGURE 10   Radial distribution functions of Na+ around carboxyl carbon in serine and around phosphorus.

The arrangement of dipoles in the headgroups and the pattern of intermolecular hydrogen bonding between NH and various oxygens in DPPS are closely related. Therefore, we first consider the rdf for H⁺ and various oxygens in the system. Figure 11 shows these rdf. From Fig. 11, a-c, it is quite clear that oxygen forms a coordination shell around NH hydrogens at distances from 1.5 to 2.3 Å, thus creating a potential for hydrogen bonding. We also calculated the coordination numbers for each of these rdf and presented them in Table 3. This table shows that the hydrogen in NH is more coordinated with the phosphate oxygen as compared to carboxyl and carbonyl oxygens (contrary to the observation by López Cascales et al., 1996). We used the criterion given by Raghavan et al. (1992) to get an estimate on the number of hydrogen bonds in the system. According to this criterion, the hydrogen bond between the donor and acceptor exists if the distance between the donor and acceptor is 3.5 Å and the angle between the vector joining the donor with the acceptor and the vector joining the donor with the hydrogen is 35^°. We find that the average number of intermolecular hydrogen bonds per molecule in the bilayer is 0.8 ± 0.3. Also from Table 3, we see that the number of hydrogen bonds between NH and phosphate group is almost twice the number of bonds between NH and carboxyl group. The observed intermolecular hydrogen bonding pattern is responsible for the reduction (compared to PC) in the area per head group of the PS bilayer. The gel phase of DPPE also shows a similar hydrogen-bonding pattern (Hauser et al., 1981). Figure 12 shows typical snapshot pictures of the intermolecular hydrogen bonds.

View larger version (4K): [in this window] [in a new window]   FIGURE 11   Radial distribution functions of various oxygens with hydrogen in NH. (a) The rdf of O in carboxyl group. (b) The rdf of O in phosphate group. (c) The rdf of O in carbonyl group. (d) The rdf of water oxygen.

View this table: [in this window] [in a new window]   TABLE 3   Coordination number and the average number of intermolecular hydrogen bonds between various oxygens in the system and hydrogens in NH

View larger version (57K): [in this window] [in a new window]   FIGURE 12   Snapshot pictures showing the intermolecular hydrogen bonds. Left: two hydrogen bonds are formed between oxygens of PO group and hydrogen of NH group. The third hydrogen bond is between carbonyl oxygen and another hydrogen of NH group. Right: two hydrogen bonds are shown. One is between oxygen of PO group and hydrogen of NH group, another is between hydrogen of NH group and serine oxygen.

Order parameters for the tails

The ordering of hydrocarbon tails can be studied with the help of NMR experiments by measuring the deuterium order parameters. The order-parameter tensor S, which is defined as
where _a is the angle made by ath molecular axis with the bilayer normal, and _ab is the Kronecker delta, can also be calculated in the simulation. In the simulations, the order parameter S_CD can be determined using the relation (Egberts and Berendsen, 1988)
Figure 13 shows the calculated S_CD order parameters for the Sn-2 hydrocarbon chain of DPPS. We also depict in Fig. 13 the experimental values for DPPC and DPPE molecules measured at roughly the same reduced temperature of ~0.07 (Thurmond et al., 1991; Douliez et al., 1995). Figure 13 shows that the order parameters for DPPS are very close to those of DPPE and the trend is consistent with the change in the area per headgroup. Notice that the order parameter profiles are similar for PS and PE.

View larger version (17K): [in this window] [in a new window]   FIGURE 13   Order parameter SCD of the Sn-2 hydrocarbon chain. A comparison of the order parameters for the DPPS, DPPE (experiment), and DPPC (experiment) at roughly the same reduced temperature of ~0.07.

Figure 14 shows the comparison between the calculated order parameter S_CD averaged over both tails of DPPS and the measured values, also averaged over both tails (Browning and Seelig, 1980). The agreement is reasonable, except for the third carbon. It is possible that the experimental value for this carbon is somewhat lower than it should be, because the depicted experimental value seems to be very close to the value for DPPC at the same reduced temperature. This is somewhat unexpected, because DPPC has a larger area per molecule, and we expect that the order parameters for DPPC should be lower than that of DPPS at the same reduced temperature. In general, one should be very careful comparing order parameters (and other quantities) for different bilayers, because they have different melting temperatures, and, moreover, in the case of DPPS, the state of the bilayer depends dramatically on the pH of the system (MacDonald et al., 1976).

View larger version (17K): [in this window] [in a new window]   FIGURE 14   Comparison between the calculated order parameters (SCD ) for DPPS averaged over both tails, experimentally measured ones and the experimentally measured order parameters for Sn-2 chain of DPPC.

Bilayer potential

To calculate the bilayer potential, we used the expression
where  is the local excess charge density in the system. We chose the zero of the potential at the center of the bilayer. The total potential has separate contributions coming from headgroups, ions, and water. To use the expression for the potential, we need to know the excess charge densities due to various components in the system. These are shown in Fig. 15 where the charge distributions due to the charges on the headgroups and Na⁺ counter ions are shown. Although our simulation is on the time scale of few nanoseconds, it is still relatively short and fails to produce a smooth curve for the charge distribution. Therefore, we fitted the calculated distributions to gaussian curves and calculated the total charge distributions from the gaussian fits. The total charge curve clearly indicates the dipolar character of the bilayer water interface. Figure 16 shows the electrostatic potential across the bilayer and the corresponding contributions. Similar to our previous observations made for hydrated DPPC bilayers, water overcompensates the contribution from headgroups and their counter ions. As a result, the total electrostatic potential in aqueous phase is negative with respect to the middle of the bilayer. This differs from the observation made in the previous simulations by López Cascales et al. (1996), and we trace the difference in potentials to the difference in counterion distributions in the two simulations. Although the total value of the potential is slightly larger in DPPS (~0.7 V in DPPS versus ~0.6 V in DPPC), the contributions from components is larger in DPPC. This means that dipole moments created by the DPPS headgroups are smaller, which is due to penetration of counterions toward the phosphate groups and the diffusional character of the counterion distribution inside membrane headgroup region.

View larger version (20K): [in this window] [in a new window]   FIGURE 15   Charge density due to DPPS and Na+.

View larger version (8K): [in this window] [in a new window]   FIGURE 16   Various contributions to the electrostatic potential across the bilayer.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION MODEL AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Our simulations on the bilayer-containing DPPS molecules in their liquid crystalline phase showed that Na⁺ counterions play an important role in determining the properties of the bilayer. The main motif appearing from our interpretation of the results is that counterions mostly screen the charge of the serine group. The remaining part of the molecule is analogous to DPPE, and, therefore, we claim that properties of DPPS bilayer in the presence of Na⁺ counterions are closer to the properties of DPPE bilayer. The reduction in the area per headgroup of DPPS compared to the area observed in DPPC bilayer is due to intermolecular hydrogen bonding between the NH group on one molecule and PO group on the other molecule. We observed that the electrostatic potential across the membrane/water interface is negative in the aqueous phase (with respect to the middle of the bilayer), as it was observed for DPPC. As in DPPC, this is due to water molecules overcompensating the contribution from the membrane. We also notice that the value of the potential is close to the one observed in DPPC membrane.

The DPPS molecule is charged to a value of 1 e.s.u. and a collection of these molecules in a membrane with the surface charge density of 1 e.s.u. per 54 Å² is expected to produce strong electric fields and potentials. Nevertheless, in the DPPS membrane, we do not observe very large electric fields or potentials, because the Na⁺ counterions are condensing on the surface of the membrane. An important question, therefore, is what fraction of counterions is condensed on the membrane? The answer depends on where we place the surface separating the membrane and water. If we place this surface at the peak of the distribution for serine oxygens (~22 Å) we observe that around of the Na ⁺ ions are inside the space defined as the "membrane." This intercalation of counterions into membrane plays a very important role in the electrostatics of the membrane and the nature of its interaction with peptides and proteins. It is therefore clear why membrane properties should depend on the kind of counterion, its charge, and the pH of the solution.