Molecular dynamics simulations of charged and neutral lipid bilayers: treatment of electrostatic interactions

Molecular dynamics (MD) simulations complement experimental methods in studies of the structure and dynamics of lipid bilayers. The choice of algorithms employed in this computational method represents a trade-off between the accuracy and real calculation time. The largest portion of the simulation time is devoted to calculation of long-range electrostatic interactions. To speed-up evaluation of these interactions, various approximations have been used. The most common ones are the truncation of long-range interactions with the use of cut-offs, and the particle-mesh Ewald (PME) method. In this study, several multi-nanosecond cut-off and PME simulations were performed to establish the influence of the simulation protocol on the bilayer properties. Two bilayers were used. One consisted of neutral phosphatidylcholine molecules. The other was a mixed lipid bilayer consisting of neutral phosphatidylethanolamine and negatively charged phosphatidylglycerol molecules. The study shows that the cut-off simulation of a bilayer containing charge molecules generates artefacts; in par-ticular the mobility and order of the charged molecules are vastly different from those Presented at the Winter School of Faculty University, Poland, 28th 2003. determined experimentally. In the PME simulation, the bilayer properties are in general agreement with experimental data. The cut-off simulation of bilayers containing only uncharged molecules does not generate artefacts, nevertheless, the PME simulation gives generally better agreement with experimental data.

determined experimentally. In the PME simulation, the bilayer properties are in general agreement with experimental data. The cut-off simulation of bilayers containing only uncharged molecules does not generate artefacts, nevertheless, the PME simulation gives generally better agreement with experimental data.
Molecular dynamics (MD) simulation is a valuable method to study molecular systems of biological relevance. The method is particularly well suited to investigating the dynamical structure of biomembranes because it is characterised by atomic resolution and allows observation of the processes on the time scale of 10 -8 -10 -7 s with a time resolution of a few femtoseconds. To carry out MD simulations of a molecular system, relevant atom-atom interactions have to be defined. Of all inter-atomic interactions that are commonly included in the potential function of the system, the electrostatic ones have the largest contribution at longer distances. A correct treatment of long-range electrostatic interactions is essential for obtaining meaningful simulation results. Since full evaluation of Coulombic interactions requires enormous computational resources, various approximations have been used. The most common are (a) truncation of long-range interactions with the use of cut-offs and (b) particle-mesh Ewald (PME) summation (Essmann et al., 1995). The theoretical grounds of the two methods is discussed in depth in Sagui & Daren (1999). The use of cut-offs has been shown to lead to several artefacts in simulations of peptides (Schreiber & Steinhauser, 1992) and nucleic acids (York et al., 1995). In the case of lipid bilayers, cut-off simulations reproduce experimental data sufficiently well (e.g., Pasenkiewicz-Gierula et al., 1999;Tieleman & Berendsen, 1998). However, this observation is based on the studies of bilayers consisting of phospholipids with zero net charge.
In this paper, cut-off and PME MD simulations of bilayers formed of uncharged and charged phospholipids are compared. The study shows that MD simulation of a bilayer containing charged molecules gives reasonable results only when electrostatic interactions are calculated using the PME summa-tion method. Application of the PME method to uncharged bilayers does not significantly improve the results. The latter conclusion is at variance with the conclusion of Patra et al. (2003) who claim that bilayers made of neutral phospholipids generated in the cut-off and Ewald MD simulations were different. They used, however, an unjustified approach for the evaluation of electrostatic interactions.

Simulation systems
In this MD simulation study two phospholipid bilayers were used. One contained both neutral and charged phospholipids (POPE-POPG bilayer); the other contained only neutral phospholipids (POPC bilayer). The POPE-POPG bilayer consisted of 54 1-palmitoyl-2-oleoyl-phosphatidylethanolamine (POPE), 18 1-palmitoyl-2-oleoyl-phosphatidylglycerol (POPG) (25 mol% POPG), and 1955 water molecules, and of 18 sodium ions (Na + ). The cations were added to neutralise the net charge of -1.0 in units of electronic charge (e) on each POPG molecule. Details concerning the membrane construction are given in ). An initial 11-ns simulation of the POPE-POPG bilayer was carried out using AMBER 4.0 (Pearlman et al., 1991). In this part of the simulation, a residue based cut-off distance of 12 Å was used to evaluate the nonbonded interactions. The last 6-ns fragment of the generated trajectory was analysed. The simulation was continued for another 14 ns using AMBER 5.0 (Case et al., 1997). In this part of the simulation the electrostatic interactions were evaluated using the PME summation method. The last 10-ns fragment of this trajectory was analysed.
The POPC bilayer consisted of 72 1-palmitoyl-2-oleoyl-phosphatidylcholine (POPC) and 1955 water molecules. Details concerning the membrane construction are given in (Murzyn et al., 2001). An initial 11-ns simulation of the bilayer was carried out using AM-BER 4.0 and a 12 Å cut-off. The simulation was continued for 6 ns using AMBER 5.0 and a PME summation (Essmann et al., 1995) with the same parameters as in the case of the POPE-POPG bilayer (see below). Last 4-ns fragments of both trajectories (cut-off and PME) were used for analysis. Figure 1 shows the structure and numbering of atoms and torsion angles in POPC, POPE, and POPG molecules.

Simulation parameters
For POPC, POPE, and POPG, optimised potentials for liquid simulations (OPLS) parameters (Jorgensen & Tirado-Rives, 1988), and for water, transferable intermolecular potential three-site model (TIP3P) parameters (Jorgensen et al., 1983) were used. The unitedatom approximation was applied to the CH, CH 2 , and CH 3 groups of the phospholipids. The numerical values for the atomic charges of POPE and POPG are given in  and of POPC follow those in Charifson et al. (1990). Other details concerning the POPE-POPG bilayer can be found in  and the POPC bilayer in Murzyn et al. (2001).

Simulation conditions
Three-dimensional periodic boundary conditions with the usual minimum image convention were used. The SHAKE algorithm (Ryckaert et al., 1977) was used to preserve the bond lengths of the OH and NH 3 groups of water, POPE, and POPG molecules, and the time step was set at 2 fs (Egberts et al., 1994). The list of nonbonded pairs was updated every 25 steps. Restraints of a flat-bottom harmonic potential as defined and implemented in the AMBER package (Case et al., 1997) were imposed on the double bonds to prevent cis-trans isomerisation. To evaluate the nonbonded interactions at first a residue base 12 Å cut-off was used (cut-off simulation) and then a PME summation (Essmann et al., 1995) with a real cut-off of 12 Å, b spline interpolation order of 5, and direct sum tolerance of 10 -6 was used (PME simulation).
The MD simulations were carried out at constant pressure (1 atm) and at the temperature of 310 K (37°C), which is above the main phase transition temperature for POPC (-5°C) (Seelig & Waespe-Šarèeviae, 1978), POPG (-4°C) (LIPIDAT), and POPE (26.1°C) (Huang & Li, 1999) bilayers. The temperatures of the solute and solvent were controlled Vol. 50 MD simulations of lipid bilayers 791 independently. Both the temperature and pressure of the system were controlled by the Berendsen method (Berendsen et al., 1984). The relaxation times for temperatures and pressure were set at 0.4 and 0.6 ps, respectively. The applied pressure was controlled anisotropically, each direction being treated independently with the trace of the pressure tensor kept constant at 1 atm.

RESULTS
The average values given below are ensemble and time averages obtained from the block averaging procedure. Errors in the derived average values are standard error estimates.

Surface area per phospholipid and bilayer thickness
A two-dimensional Voronoi tessellation method (Shinoda & Okazaki, 1998) allows one to compute the surface area of a single molecule in the bilayer made of lipids of a single type (cf. Murzyn et al., 2001). In the case of a mixed lipid bilayer, the method is less precise, nevertheless it was applied to estimate the surface area of each POPE and POPG molecule in the POPE-POPG bilayer. From the estimated single molecular areas, the average values for POPE and POPG in both cut-off and Ewald simulations were obtained (Table 1). The average surface area per POPC in the POPC bilayer was obtained in a straightforward way by dividing the total surface area of the bilayer by 36, which is the number of molecules in each leaflet of the bilayer. The values for both cut-off and Ewald simulations are given in Table 1.
The bilayer thickness was estimated from the distance between average positions of the phosphorus atoms in the two leaflets of a bilayer (P-P distance). The values for both cut-off and Ewald simulations are given in Table 1.
Both the surface area per phospholipid and the bilayer thickness are more affected by the simulation protocol in the case of the POPE-POPG bilayer than the POPC bilayer, but the effect is not very large.

Bilayer/water interface
Intermolecular interactions that take place at the bilayer/water interface include hydrogen bonding between water molecules and polar groups of POPC, POPE, and POPG; direct H-bonding and charge pairing among POPC, POPE, and POPG headgroups; and clathrating of water molecules around the choline group of POPC. A water molecule that is simultaneously H-bonded to two phospholipid molecules forms a water bridge between them, whereas a charge pair is formed between a negatively charged oxygen atom of one POPC and a positively charged choline moiety of a neighbouring POPC. In the following analyses, the same geometrical definitions of H-bonding, water bridging, and charge pairing as in our previous papers Murzyn et al., 2001) were used. In addition to these short distance Table 1

. Bulk properties of POPC and POPE-POPG bilayers in Ewald and cut-off simulations.
Surface area per phospholipid and bilayer thickness (P-P distance) are listed with standard error estimates. intermolecular interactions, in the POPE-POPG bilayer we define charge pairs between Na + ions and oxygen atoms of phosphate, carbonyl, or glycerol groups (cf. Fig. 1), when the ion and the oxygen atom are no more than 3.0 Å apart.

POPC bilayer
Numbers of water molecules that interact with POPC as well as of POPC-POPC charge pairs for the two simulation protocols are given in Table 2. The numbers are practically unaffected by the simulation protocol, as differences between the corresponding numbers for the cut-off and Ewald simulations are within estimated errors.

POPE-POPG bilayer
Numbers of water molecules that interact with POPE and POPG, of direct H-bonds between phospholipids, and of charge pairs between phospholipids and Na + for the cut-off and Ewald simulations are given in Table 3. In contrast to the POPC bilayer, the numbers of intermolecular interactions at the bilayer/water interface of the POPE-POPG bilayer strongly depend on the simulation protocol. The number of water molecules H-bonded to a phospholipid is significantly higher, whereas the number of the lipid-lipid H-bonds is smaller in the cut-off simulation. A considerable difference in the number of Na + -lipid charge pairs for the two protocols is noticeable. In the cut-off simulation, Na + practically do not interact with POPE but do very frequently interact with POPG. On average, 50% of all Na + ions are bound to POPG at any instant. In the Ewald simulation, only 30% of all Na + ions are bound to lipids, 20% to POPG and 10% to POPE (Table 3).

Bilayer core
The molecular order parameter for the nth segment of an alkyl chain, S mol , is defined through (Hubbell & McConnell, 1971) S mol = 1 2 3 1 2 cos -q n q n is the instantaneous angle between the nth segmental vector, i.e., the (C n-1 , C n+1 ) vector linking the n -1 and n + 1 carbon atoms in the alkyl chain and the bilayer normal; áñ denotes the ensemble and time average. S mol profiles along the b-chain of POPC, POPE, and POPG molecules are shown in Fig. 2. Details concerning calculation of conformation-related quantities are given in Róg and Pasenkiewicz-Gierula (2001) and . Mean values of the order parameter for the b-chain for the cut-off and Ewald simulations are given in Table 4. The order of the b-chain is, in general, higher for the PME than the cut-off simulation. A particularly large difference is observed for the b-chain of POPG (Fig.  2 c, and Table 4).
Other parameters characterising the hydrocarbon chain region of the bilayer, like the b-chain tilt angle (for definition see Róg et al., 2003*) and the number of gauche rotamers in Vol. 50 MD simulations of lipid bilayers 793 the b-chain are also given in Table 4. These parameters are affected only mildly by the simulation protocol.

Lateral diffusion of phospholipids
Diffusion coefficients for the lateral self-diffusion of POPC, POPE, and POPG molecules in the bilayer plane (D) were calculated from linear parts of the respective mean square displacement (MSD) curves for the centre of mass of all lipids of a given type in the membrane (Fig. 3). The coefficients for the cut-off and Ewald simulations are compared in Table  4. For POPC and POPE, the differences between the respective values of D for the two protocols are small. For POPG, the value of D in the cut-off simulation is 35 times smaller than in the Ewald simulation. A good illustration of the huge difference between the diffusion of POPG in the cut-off and Ewald simulations are the trajectories of the centre of mass of a selected POPG molecule displayed in Fig. 4. In the cut-off simulation, the centre of mass of the POPG practically remains in the same position during a period of 4 ns, whereas 794 T. Róg and others 2003 Table 3.

Intermolecular interactions at the bilayer/water interface for Ewald and cut-off simulations of the POPE-POPG bilayer.
Listed are average numbers of water molecules H-bonded to a lipid molecule, of direct lipid-lipid H-bonds, of water bridges, and of Na + -lipid charge pairs per POPE and POPG molecule. Errors are standard error estimates.

Table 4. Average values of parameters for Ewald and cut-off simulations of POPC and POPE-POPG bilayers.
Average values of molecular order parameter (S mol ), tilt angle, and number of gauche conformations per chain, for the b-chain as well as the lateral diffusion coefficient (D) for the centre of mass of POPC, POPE, and POPG are given with standard error estimates.
in the Ewald simulation, in this period it covers a large area.

DISCUSSION
In this study, the effect of the MD simulation protocol on the properties of a bilayer made of uncharged and charged phospholipids has been evaluated. Two simulation protocols were applied. In one, long-range electrostatic interactions were truncated with the use of a residue-based 12 Å cut-off, in the other, the interactions were evaluated with the use of the PME summation method. Two bilayers were simulated and analysed. One consisted of neutral POPC molecules, the other consisted of neutral POPE and charged (-1.0 e) POPG molecules.
The study clearly shows a strong effect of the simulation protocol on the behaviour of POPG but not POPE, although both POPGs and POPEs in the proportion 1:3 comprise one bilayer. The most striking result is the almost complete immobilisation of POPG molecules in the cut-off simulation as compared to POPE molecules in the same simulation and to their mobility in the PME simulation. In the cut-off simulation, the diffusion coefficient for POPG self-diffusion on the nanosecond time scale is over 30 times smaller than that of POPE and POPC, and than that of POPG in the PME simulation. The diffusion coefficients of PCs measured experimentally by neutron scattering (Tabony & Perly, 1991) and electric dipole relaxation (Haibel et al., 1998) methods, which report events on a similar time scale as MD simulations, are similar to those of POPC and POPE obtained in both cut-off and PME simulations as well as to the value for POPG obtained in the PME simulation.
The simulation protocol also affects the molecular order parameter profile of POPG rather strongly but not POPE. S mol profiles for hydrocarbon chains are available in the literature for phosphatidylcholines (PC) and phosphatidylethanolamines (PE) but not for Vol. 50 MD simulations of lipid bilayers 795  phosphatidylglycerols (PG). The S mol profile of the POPG b-chain obtained in the cut-off simulation does not resemble the S mol profiles derived experimentally for saturated (Trouard et al., 1999) or unsaturated (Seelig & Waespe-Šarèeviae, 1978) PC's or for PE's (Marsh et al., 1983). However, that obtained in the PME simulation does resemble these profiles. The S mol profiles for POPE obtained in both cut-off and PME simulations show a similar resemblance to experimental data. An effect of the simulation protocol on the intermolecular interactions in the interfacial region of the POPE-POPG bilayer is also visible. In this region, both POPE and POPG are similarly affected. The cut-off simulation "over-produces" lipid-water H-bonds and "under-produces" lipid-lipid H-bonds. Consistently with the immobilisation of POPG in the cut-off simulation, the Na + ions are predominantly in the bound state with POPG. In the PME simulation, the Na + ions are more mo-bile with fewer bound to the lipids; moreover, they interact with both POPE and POPG.
The effect of the simulation protocol on the POPC bilayer is much smaller than in the case of the POPE-POPG bilayer. Intermolecular interactions at the bilayer/water interface are practically unaffected. The order of the POPC alkyl chains better reproduces experimental data (Seelig & Waespe-Šarèeviae, 1978) in the PME simulation but the difference is not large. In a recent paper, results of 12 Å cut-off MD simulations of bilayers made of saturated, trans-, and cis-unsaturated lipids were compared with available experimental data (Róg et al., 2003). The comparison indicated a good agreement between both sets of data. The first and main conclusion of our study is that cut-off simulation of a bilayer containing, among others, molecules of non-zero net charge, generates artefacts. These artefacts concern the mobility and order of the charged molecules but not of the uncharged ones within the same bilayer. In the PME simulation, the bilayer properties are in general agreement with experimental data. Thus, whenever a bilayer system contains charged molecules, reasonable results of MD simulation can only be obtained when the long-range electrostatic interactions are evaluated using the PME summation method.
The second conclusion of our study it that cut-off simulation of a bilayer containing only uncharged molecules does not generate artefacts. The parameters that describe bilayer properties, derived from both cut-off and PME simulations are similar, although in some cases the results of the PME simulation better reproduce the relevant experimental data. The MD simulation studies of uncharged bilayers that we have undertaken so far provided qualitatively correct results even though neither cut-off nor PME simulations reproduce faithfully all relevant experimental results. To minimise the effects of methodological inaccuracies and other uncertainties, it is advisable always to conduct comparative MD simulation studies in which two or more bilayer systems are simulated using the same protocol with one of the systems earlier having been compared successfully with appropriate experimental data. As the possible errors in parameters or methodology are most likely to be the same for all systems, such studies can be expected to give a faithful account of both similarities and differences among the systems.
In a recent paper, Patra et al. (2003) concluded that the properties of a bilayer made of neutral dipalmitoylphosphatidylcholine molecules generated in the cut-off and Ewald simulations were different. In those simulations, several parameters and definitions differed from those used by us. Particularly affecting the time development of the simulated system was their definition of a group in the groupbased (our residue-based) cut-off simulation, where the sum of charges of all atoms comprising the group was far from zero, as well as the use of cut-off distances much larger than those for which OPLS parameters were derived. These unjustified choices caused significant discrepancies in the behaviour of their and our systems and resulted in opposite conclusions.