Biomembrane simulations of 12 lipid types using the general amber force field in a tensionless ensemble

The AMBER family of force fields is one of the most commonly used alternatives to describe proteins and drug-like molecules in molecular dynamics simulations. However, the absence of a specific set of parameters for lipids has been limiting the widespread application of this force field in biomembrane simulations, including membrane protein simulations and drug-membrane simulations. Here, we report the systematic parameterization of 12 common lipid types consistent with the General Amber Force Field (GAFF), with charge-parameters determined with RESP at the HF/6–31G(d) level of theory, to be consistent with AMBER. The accuracy of the scheme was evaluated by comparing predicted and experimental values for structural lipid properties in MD simulations in an NPT ensemble with explicit solvent in 100:100 bilayer systems. Globally, a consistent agreement with experimental reference data on membrane structures was achieved for some lipid types when using the typical MD conditions normally employed when handling membrane proteins and drug-membrane simulations (a tensionless NPT ensemble, 310 K), without the application of any of the constraints often used in other biomembrane simulations (such as the surface tension and the total simulation box area). The present set of parameters and the universal approach used in the parameterization of all the lipid types described here, as well as the consistency with the AMBER force field family, together with the tensionless NPT ensemble used, opens the door to systematic studies combining lipid components with small drug-like molecules or membrane proteins and show the potential of GAFF in dealing with biomembranes.

Biological membranes are crucial structures to living organisms. They mediate different physiological processes within cells, from transport phenomena to signalling events, and they serve as the interface between the intracellular and the extracellular environment. Biological membranes are also of high pharmacological interest, especially considering that more than 30% of the drug targets (such as enzymes, nuclear receptors or nucleic acids) used in current therapies (Drews, 2000) are located within the cell and across the plasma membrane. In the conventional pharmaceutical industry and in modern medicine, it is well known that many promising drugs that are discovered never make it to the market because of difficulties in delivery (Scherrmann, 2002), so the development of methodologies able to evaluate with atomistic detail their interaction with biomembranes is of great importance.
Although a large number of experimental studies involving biomembranes have been performed, yielding important information on a variety of properties and processes, for many biological problems a more detailed atomic-level understanding is often required. MD simulations offer an ideal complement to such experimental studies, providing detailed insight into the interaction and dynamics of lipids and of their interplay with drug-like molecules and biomembrane proteins. For this reason, MD simulations have become an important tool in the atomic-level study of biomembranes (Aller, Garnier, & Genest, 2006;Ash, Zlomislic, Oloo, & Tieleman, 2004;Beevers & Kukol, 2006;Cordomí, Caltabiano, & Pardo, 2012;Nitschke, Vequi-Suplicy, Coutinho, & Stassen, 2012;Raghunathan et al., 2012;Samanta, Hezaveh, Milano, & Roccatano, 2012;Stansfeld & Sansom, 2011;Tieleman, 2006;. Three main types of force fields have been available for the study of biomembranes: (1) coarse-graining force fields; (2) united-atoms force fields; and (3) all-atom force fields. Coarse-grained force fields describe large molecular units (such as amino acid side chains and whole or multiple water molecules) as a single particle, while united-atom force fields combine each aliphatic carbon and associated hydrogen atoms into single particles and are significantly more detailed. All-atom force fields (i.e. atomistic) treat every atom (including all hydrogen atoms) explicitly and are the most appealing type for force field when detailed atomic-level information is sought (Marrink, de Vries, & Tieleman, 2009).
Coarse-grained force fields (De Nicola, Zhao, Kawakatsu, Roccatano, & Milano, 2011;Marrink, Risselada, Yefimov, Tieleman, & de Vries, 2007;Wang & Deserno, 2010) were important in the onset of computational lipid simulations and continue to be a useful alternative to study effects that require a very long computational time or very large systems. Nevertheless, throughout the years a variety of united-atoms force fields have been developed. Examples include the Berger lipids (Berger, Edholm, & Jähnig 1997), the G53A6L FF (Poger, Van Gunsteren, & Mark 2010), the 43A1-S3 FF (Chiu, Pandit, Scott, & Jakobsson 2009) and the potential developed by Ulmschneider et al. (Ulmschneider & Ulmschneider, 2009) and Kukol (Kukol, 2009). More recently, however, considerable emphasis has been put into the development of all-atom force fields for lipid simulations. Atomistic molecular simulations retain virtually all atomic-level interactions and can use time-steps in the femtosecond range. While this makes them quite slow and expensive compared with the other alternatives, the wide range of carefully tested parameters available for these models, including proteins, nucleic acids and small organic molecules makes them reliable when it comes to quantitative prediction of properties such as motional timescales or interaction strengths (Lindahl & Sansom, 2008), showing that this type of simulations has great advantages, over those with lower detail level.
While fully atomistic parameters for different lipidic components have been made available over the last few years, particularly based on CHARMM (Högberg, Nikitin, & Lyubartsev, 2008;Hyvönen & Kovanen, 2005;Klauda, Brooks, MacKerell, Venable, & Pastor, 2005;Klauda et al., 2010;Sonne, Jensen, Hansen, Hemmingsen, & Peters, 2007) and in AMBER force fields (Jójárt & Martinek, 2007;Jämbeck & Lyubartsv, 2012a, 2012bRosso & Gould, 2008;Siu, Vácha, Jungwirth, & Böckmann, 2008;Skjevik, Madej, Walker, & Teigen, 2012), such attempts are often restricted to one or two individual lipid components, with different studies adopting quite different assumptions and different methodological variations. Also, an important number of these parameter sets was designed to describe lipid membranes when a positive surface tension is applied (NγPT ensemble) or when the total simulation box area is restricted (NPAT), solutions that restrain membrane undulations and stretching (Jämbeck & Lyubartsev, 2012a) and are clearly not ideal for simulating drugbilayer or protein-bilayer systems, as they can lead to severe artefacts. In fact, for such systems, the much more realistic NPT ensemble is the preferred choice. Although the application of a surface tension is not desirable, inaccuracies due to the inclusion of a surface tension are often favoured in detriment to simulating bilayers with parameters that lead to under or overestimation of structural parameters; therefore many studies have been conducted with this type of ensembles with success (Gullingsrud & Schulten, 2004;Leontiadou, Mark, & Marrink, 2007;Siu et al., 2008;Esteban-Martín, Risselada, Salgado, & Marrink, 2009;Vácha et al., 2009;Muddana, Gullapalli, Manias, & Butler, 2011;Skjevik et al., 2012).
The objective of this work hence lies in determination of a fully consistent set of parameters for a set of 12 representative lipid types, in a common framework with the General AMBER force field (GAFF) (which was developed for drug-like molecules) (Wang, Wolf, Caldwell, Kollman, & Case, 2004) and with the AMBER force fields that are specialized in proteins (Cheatham, Cieplak, & Kollman, 1999;Duan et al., 2003) allowing the simulation in a consistent fashion of biological systems comprised of structural elements from these different classes. Particular effort was dedicated to ensure that the present set of parameters is able to reproduce experimental reference data on lipid bilayers without the need to constraint the surface tension or area of the bilayer or the total simulation box area to improve the results, i.e. in the more realistic NPT tensionless ensemble.
The lipid models chosen are pure-species bilayers with variations in the most common type of phospholipids: phosphatidylethanolamine (PE) and phosphatidylcholine (PC). They vary in terms of the degree of saturation and length of the acyl chains. We have performed the biomembrane simulations at 310 K, since this is the temperature at which biological systems are most frequently simulated, including proteins, nucleic acids, drug-like compounds and the combination of all these components with biomembranes. Properties analysed included the volume per lipid, the area per lipid and the bilayer thickness. A comparison of the electron density profiles, the order parameters for the acyl chains and an analysis of the hydration at the water/lipid interface are also performed. In addition, the dynamic properties of the systems, evaluated through the coefficients of lateral diffusion of the lipid molecules, are explored too.

Force field parameters
In this work, we have used the GAFF as a basis for the parameterization of the different lipid components. This force field was not developed specifically for lipid systems; however, it is advantageous since it can be used later on in a consistent fashion in the treatment of mixed systems such as drug-membrane or even protein-membrane simulations. The attribution of GAFF (Wang et al., 2004) atom types was performed using the Antechamber (Wang, Wang, Kollman, & Case, 2006) module in the AMBER 10.0 package (Case et al., 2008) for all the 12 types of phospholipids used ( Figure 1). The point charges for each of the studied phospholipids were determined using GAUSSIAN09 (Frisch et al., 2009), at the HF/6-31G(d) level, considering a Restrained Electrostatic Potential (RESP) (Bayly, Cieplak, Cornell, & Kollman, 1993) protocol, starting from B3LYP/6-31G(d) initially optimized structures. The use of HF/6-31G(d) and of RESP charges can be justified by the need for consistency with the majority of the AMBER force fields, namely GAFF for small drug-like molecules, and with the protein force field for amino acids (Cheatham et al., 1999;Cornell et al., 1996;Duan et al., 2003).

Initial structures
The bilayer systems used in this work are fully atomistic systems with 200 phospholipids (100:100) plus explicit water molecules, in a total of 12 different bilayer setups. These bilayer systems and their composition are described in Figure 1 and Table 1. The initial structure was generated using the packmol program (Martínez, Andrade, Birgin, & Martínez, 2009). This program packs molecules in defined regions of space and guarantees that short repulsive interactions do not disturb the initial stage of simulation. Basically, through packmol, a single phospholipid was replicated along the x and y axes and flipped along the z axis, until a rectangular grid of 200 phospholipid molecules was obtained. The user has to specify the size of the grid and some spatial restrictions for a proper orientation of the molecules within that grid.
The program then distributes and pre-equilibrates the molecules under the conditions imposed. LEAP was then used to add a 15 Å layer of TIP3P type water molecules above and below the bilayer (considering the z axis). In the acyl chain's double bond containing phospholipids (PLPC, POPC, DOPC, PLPE and POPE), the double bond is in a cis conformation which is the more relevant isomer in a biological environment (van Meer, Voelker, & Feigenson, 2008). Different approaches to generate the initial lipid conformation were also tested, including biomembranes generated with VMD's (Humphrey, Dalke, & Schulten, 1996) membrane builder and preequilibrated membranes from studies of other authors (Heller, Schaefer, & Schulten, 1993;Rosso & Gould, 2008).

Simulation conditions
After the initial structure has been set up, a geometry optimization was performed in order to eliminate bad contacts. For this purpose, the SANDER module from the AMBER package 10.0 was used. The minimization was performed in three stages, combining both steepest descent and conjugate gradient minimization algorithms. In a first stage (10,000 steps) only the molecules of water were minimized. In the second stage (10,000 steps) the minimization of the hydrogen atoms was performed. The remaining of the system was maintained rigid with 50 kcal mol À1 Å À2 harmonic forces. In the third and final stage (80,000 steps) all the system was minimized. Molecular dynamics simulations were performed with the PMEMD (Particle Mesh Ewald Molecular Dynamics) module from the AMBER package 10.0. Protonation states at pH = 7.0 were considered for all the phospholipids. All the simulations were performed using an NPT ensemble, and for these simulations a constant temperature of 310.0 K was set throughout the simulation, using the weak-coupling algorithm (Berendsen, Postma, Vangunsteren, Dinola, & Haak, 1984). Although this temperature is below some of the phase transition temperatures (T m ) of some of the glycerophospholipids used (DPPC, DPPE, DMPE, and DSPE), it is frequently the temperature of choice when performing MD simulations involving biological systems, particularly those involving the behaviour of proteins or drug-like molecules interacting with biomembranes. Also, despite the fact that this temperature is different from the available experimental data (by approximately 7°C), the typical variation associated to this difference is not very high (ca. 0.08 Å/°C for the bilayer thickness and ca. 0.20 Å 2 /°C for the area per lipid) (Nagle & Tristram-Nagle, 2000a;Kučerka, Nieh, & Katsaras, 2011), when dealing with temperatures far from phase transition regions. For the pressure regulation we have used anisotropic conditions and the reference pressure was set to 1 atm, using a pressure relaxation time of 0.2 ps. The same variety of weak-coupling algorithms was used as in the temperature regulation, using 1.0 ps of relaxation time. The Langevin temperature scheme (Pastor, Brooks, & Szabo, 1988), with a collision frequency of 1.0 ps À1 , was also tested. The SHAKE algorithm (Ryckaert, Ciccotti, & Berendsen, 1977) was employed in all simulations, in order to constraint only bonds involving hydrogen, allowing for the use of an integration time-step of 2 fs. Periodic boundary conditions were used and the cut-off radius for non-bonded interactions under PME (Essmann et al., 1995) protocol was set to 10 Å. The length of the MD simulations comprised 40 ns. In selected cases these simulations were extended to 80 ns (Table 1). The MD trajectory was saved every 2 ps. All of the MD results were analysed with the PTRAJ module of AMBER 10.0 and with the VMD program version 1.8.6 graphical interface (Humphrey et al., 1996).

Results and discussion
The performance of the lipid parameters developed here was evaluated through a structural analysis at the atomic level for all the bilayer systems described in Table 1, taking into particular attention three important quantities: volume per lipid, area per lipid and bilayer thickness. These properties are very sensitive to the inter and intramolecular interactions, which makes them a good reference to evaluate the quality of the force fields. Moreover, there are experimental values in the literature that can be used to benchmark the computational results. For the six systems for which experimental reference data was available (DLPC, DLPE, DOPC, POPC, POPE, DPPC, and DMC) Kučerka et al., 2011;McIntosh & Simon 1986;Nagle & Tristram-Nagle, 2000b;Rappolt, Hickel, Bringezu, & Lohner 2003), direct comparisons were performed in order to validate our models. In addition, an analysis of the hydration at the water/lipid interface is pursued, as well as an assessment on the lateral diffusion of the phospholipid molecules. The inherent disorder in the acyl chains was also analysed through the determination of computed deuterium order parameters (S CD ). For DOPC and DMPC, for which the AMBER parameters have been recently made available in the literature in studies adopting a similar protocol (Rosso & Gould, 2008;Siu et al., 2008;Tessier, DeMarco, Yongye, & Woods, 2008;Skjevik et al., 2012), RESP charges and the structural properties were also compared. Additionally, Lennard-Jones (LJ) parameters were related with the recent work conducted by Jämbeck & Lyubartsev (Jämbeck & Lyubartsev, 2012a) and the GLYCAM06 force field (Tessier et al., 2008). Finally, the influence of using different membrane starting structures and different thermostats was also evaluated with results provided with the supporting information of this article.

Volume per lipid
The volume per lipid is a very important quantity to take into consideration when comparing simulations to experimental data, as it shows fewer fluctuations and converges faster than the area per lipid (Anézo, de Vries, Höltje, Tieleman, & Marrink, 2003;Rosso & Gould, 2008). The volume per lipid (V L ) was calculated using Equation (1): where V box is the volume of the simulation box, V W is the volume per water molecule, n W is the number of water molecules used in the simulation and N Lipids is the number of lipid molecules. Regarding the V W constant, we have used the value of 30.53 Å 2 , as reported by Rosso and Gould (Rosso & Gould, 2008) from a simulation of 361 TIP3P water molecules at T = 310 K and p = 1 atm. In Figure 2a and 2b, we can observe the evolution of this quantity for the different bilayer systems. This quantity shows great stability throughout the simulation for each of the tested systems. Considering the convergence time, this is reached within the first 10 ns of the MD simulations. When comparing with available experimental data, we have found great similarity between values. The smallest variation from experimental techniques observed, when comparing the average of the last 20 ns of the MD simulations with the experimental data, is for the DOPC and POPC systems (below 2.5% difference), while for the others (DLPC, DLPE, DMPC, DPPC and POPE) a difference between 3 and 9% is obtained. For all the systems this quantity is below the experimental value, as can be observed in Table 2.
Comparing the PC (DLPC, DPPC, DMPC, PLPC and POPC) and PE systems (DLPE, DPPE, DMPE, DSPE, PLPE and POPE), we perceive that with the increment on the acyl chains' carbon length, the volume per lipid increases, as expected. We also observe identical values for the PLPC and POPC systems, as well as for the PLPE, POPE and DSPE systems, justified by their greater similarity on the acyl chains' composition. The trend is also for a smaller volume per lipid for the PE systems, when comparing with analogues (identical acyl chains) in the PC systems, due to the smaller size Table 1. Overview of the developed simulation systems concerning the two glycerophospholipid acyl chains' carbon skeletons (sn-1 and sn-2); the total number of atoms in the system; and the length of the simulations. of the PE headgroup. Hence, these parameters are able to reproduce the relative volume per lipid of the different lipid types. Interestingly, a greater similarity of the MD results to the experimental values at 20°C for the DPPC and DLPE systems was obtained (ca 2% and 1% difference, respectively). This last phenomenon will be discussed ahead.

Area per lipid
The area per lipid is another useful quantity to evaluate the equilibration of the bilayer systems or to monitor phase transitions (Siu et al., 2008). In this case, it was determined through the rectangle covering the lipid's nitrogen atoms projected onto the xy plane. The area per lipid is given by dividing the area of the xy rectangle by the total number of lipids in one layer (100). In Figure 2c and 2d, we have a representation of the variation of this quantity throughout the entire simulation time. When analysing Figure 2c and 2d, we observe a very stable behaviour of this quantity for all the studied systems. However, when comparing with the volume per lipid entity, the convergence time is much slower (near 20 ns). This supports the need of increasing the time length of the MD simulations for this type of systems (Anézo et al., 2003).
Regarding the comparison of the calculated area per lipid with the experimentally determined data, the value for this structural parameter is consistently being underestimated, which is coherent with literature observations that have used GAFF (Jójárt & Martinek, 2007;Rosso & Gould, 2008;Siu et al., 2008;Skjevik et al., 2012) and simulations in an NPT ensemble without restraining the surface tension or the area per lipid. What is also very interesting is the fact that the results obtained for DOPC show greater consistency with experimental data (ca. 2 to 8% variation, depending on the literature data considered) and improved accuracy when compared to other MD studies with this type of system and under similar protocols (Rosso & Gould, 2008;Siu et al., 2008). The consistency with experimental data in the DOPC situation could probably be justified by the increasing disorder of the hydrocarbon skeleton, more consistent with L α phases, as can be seen in Figure 3. However, we could not observe the same tendency for the mono-acyl unsaturated POPC, POPE, PLPC and PLPE systems.
When comparing the PE systems, we can see lower results for the area per lipid quantity justified by the smaller headgroup region, when compared with the PC systems. Also, we observe that with the unsaturation of the acyl chains the area per lipid increases, reaching a maximum result in the DOPC situation (66 Å 2 ), where both acyl chains are unsaturated. These observations are in line with the experimental data available, as seen in Table 2.
This underestimation is also consistent with a probable transition of the bilayer to a gel-like state that is also supported by the greater consistency of values, when considering experimental data, determined at lower temperatures (20°C), for the DPPC and DLPE systems, with a difference to the experimental data of ca. 6% and ca. 2%, respectively.

Electron density profiles and bilayer thickness
The bilayer thickness can be calculated through the peak-peak distance of the density profile considering the hydrated system, defined as D HH . However, since these peaks are due to the highly dense phosphate groups (Kučerka et al., 2008;Rosso & Gould, 2008), the bilayer thickness was determined here through the distance of the opposing phosphorus atoms in each layer. This value was used in the comparison with the experimental data. The bilayer thickness shows great stability throughout the simulation protocol with a convergence time close to that observed for the area per lipid quantity (ca. 20 ns), as can be seen in Figure 2e and 2f.
Again, comparing the PC (DLPC, DPPC, PLPC and POPC) and PE systems (DLPE, DPPE, DMPE, DSPE, PLPE and POPE) we notice that with the increment on the acyl chains' carbon length the bilayer thickness also increases, in similarity with the volume per lipid. The exception lies in the DMPC case. We also observe identical values for the PLPC and POPC systems, as well as for the PLPE and POPE systems, justified by their greater similarity on the acyl chains' composition.
Contrarily to the area per lipid determination, when applying this formalism to estimate the bilayer thickness there is a consistently overestimation of this quantity, regarding the experimental data, as can be observed in Table 2. Our best result is again the DOPC situation (ca. 6% difference from experiment).
Electron density profiles (EDP) were calculated using a VMD analysis plugin (Giorgino). Regarding the electronic density profiles, we found a maximum density in the 0.43-0.52 e.Å À3 interval and a minimum density in the 0.18-0.24 e.Å À3 interval. When comparing the PE saturated systems (Figure 4a), we see an expected behaviour for the EDP profile, as the PE systems with longer carbon skeletons of the acyl chains show a higher peak-peak separation. Moreover, the electron density profiles show great similarity within the hydrophobic region of the profile, except for the DLPE system, where the hydrophobic region appears with a lower electron density, and the peaks with higher electron densities. Regarding the saturated PC systems (DLPC, DMPC and DPPC), the same behaviour is detected as in the PE situations, where we also observe the same hydrophobic region's lower density and higher peak densities for the DLPC system (Figure 4b), analogous to what is observed for the DLPE system.
When comparing the unsaturated systems, Figure 4c, we do not see the same differences as for the previous situations, probably due to the higher similarity of the acyl chains (variation only in the degree of saturation). Focusing on the DOPC system, we can observe the peak-peak distance contraction, in agreement with the lower bilayer thickness determined for this case, as well as a wider distribution of the electron density profile throughout, when compared with the mono-acyl unsaturated phospholipids (PLPC, PLPE, POPC and POPE), translating into lower maximum densities and higher lower densities, respectively, 0.43 and 0.24 e.Å À3 .
Considering identical phospholipids only with different head groups (PE vs. PC), we have found concise variations between them, for each of the structural results. For the volume per lipid, the variation between them is 5-8%; for the area per lipid, the difference is among the interval 11-17%; and for the bilayer thickness this difference is in the order of 3-10%. The area per lipid variation between PC and PE is also consistent with the literature data, suggesting a variation of ca. 10-20% Figure 3. Representations of the last frame of the 80 ns dynamics for: (a) DLPE system and (b) DOPC system. Hydrogen atoms are omitted for clarity. The carbon atoms from the acyl chains are represented in yellow, while the carbon atoms from the polar headgroups in green. Water molecules are shown in blue. These representations were produced with the PYMOL software (Schrödinger, 2010). (Pandit & Berkowitz, 2002). This supports the consistency of the obtained parameters.

Radial distribution functions
Water molecules play a key role in the assembly of biological membranes (Berkowitz, Bostick, & Pandit, 2006). In fact, the water molecules at the water/lipid interface are even considered an integral part of a biomembrane (Jójárt & Martinek, 2007). We have performed an analysis of the water/lipid interface through the determination of radial distribution functions (RDF). Integrating the RDF curves, we can determine the number of water molecules in the hydration shells for specified atoms in the system.
In Figure 5a, we can see the graphical representation of the RDF for the distance between the oxygen atoms of the TIP3P water molecules and the closest atom in the DOPC headgroup. The first peak of the RDF was found at 2.55 Å and a second small peak can be observed in the 5 Å region. Integrating the RDF curve, we have found a number of ca.12.7 water molecules surrounding each DOPC headgroup. This result is in agreement with experimental values (Hristova & White, 1998) and previous MD simulations (Mashl, Scott, Subramaniam, & Jakobsson, 2001;Rosso & Gould, 2008). For the other cases, a hydration result of 6-8 water molecules for the first hydration shell was found when considering the PE systems, while for the PC systems an hydration shell of 8-11 water molecules was determined (when using an integration radius of ca. 3.3 Å). In the literature, the suggested number of water molecules per lipid is in the range of 10-30 water molecules, depending on the state of the alkyl chains (solid or liquid) (Wennerstrom & Sparr, 2003). The lower hydration for the PE systems is expected, when compared with the PC systems, due to the lower numbers obtained for the area per lipid quantity.  We have also determined the RDFs for the water molecules and the closest atoms of the different chemical groups of the phospholipids (not shown), where we have found consistent results between the systems. In the PE systems, the number of water molecules in the first hydration shell for the ethanolamine, phosphate, glycerol and carbonyl groups is, respectively, 6.6, 3.5, 2.7 and 1.1 (averaged between the PE systems), while in the PC systems (except in the DOPC and DMPC situation), for the choline, phosphate, glycerol and carbonyl groups, the averaged water molecules in the first hydration shell are, respectively, 9.3, 3.8, 2.7 and 1.0. Considering the DOPC system the results are correspondingly 11.5, 4.2, 3.6 and 1.3, while for the DMPC system they are 11.0, 4.7, 3.2 and 1.1; consistent with the data obtained by Rosso and Gould (Rosso & Gould, 2008). The integration radius used for both the PE and the PC systems were, respectively, ca. 3.5, 3.2, 3.6 and 3.2 Å.
Additionally, RDFs were constructed for the closest atom of the amine group of ethanolamine containing phospholipids (PE) and the oxygen atoms of the phosphate and carbonyl groups (data not shown) of surrounding phospholipids, where we have found peaks for the profiles in the 1.75-1.90 Å interval, with a minimum at ca. 1.6 Å.
Interestingly, we have found a consistent behaviour for the RDF's of the water molecules and the phosphorus and nitrogen atoms regarding the PE and PC systems, with PLPC and PLPE systems depicted as an example in Figure 5b. For the PC systems, the peak of the RDF for the phosphorus atoms appears sooner than for the nitrogen atoms and the width is smaller, suggesting a higher ordering of the water molecules around the phosphate group, possibly due to hydrogen bonding between the water molecules and the oxygen atoms of the phosphate group. The wider RDF profile for the nitrogen atom, despite the positive charge, could be due to the hydrophobic character provided by the surrounding bulky methyl groups (Hyvönen, Öörni, Kovanen, & Ala-Korpela, 2001). For the PE systems we observe a narrower profile for the nitrogen atoms, suggesting the existence of hydrogen bonds with water molecules, while for the phosphorus atom the width is identical to the PC systems.

Diffusion coefficients
We have also determined the diffusion coefficients for the lipidic systems, averaged for the last 4 ns of the MD simulations, considering the plane of the bilayer, defined as the lateral diffusion coefficient. The movements along the z axis (normal to the bilayer) are not considered as diffusion but as oscillations (Zlenko & Krasil'nikov, 2008). The diffusion coefficients can be determined from the slope of the linear part of the plots of mean square displacement (MSD) of the centre of mass of the lipids vs. time, as depicted in Figure 6. The lateral diffusion coeffi-cient values for all the systems, except for the DOPC and DPPE bilayers were in the 10 À8 cm 2 s À1 order of magnitude. The DOPC bilayer showed a higher lateral diffusion coefficient in the 10 À7 cm 2 s À1 order of magnitude, while the DPPE system a lower lateral diffusion coefficient in the 10 À9 cm 2 s À1 order of magnitude, as can be observed in Table 3. Considering the DOPC situation, we have found a very good agreement with experimental data at 35°C (1.37 Â 10 À7 ) (Filippov, Orädd, & Lindblom, 2003). Lower lateral diffusion coefficients were observed for the other situations, consistent with other simulations that have used GAFF (Jójárt & Martinek, 2007;Rosso & Gould, 2008;Siu et al., 2008), but approximately one order of magnitude lower than experimental determinations (Tabony & Perly, 1991;Konig, Pfeiffer, Bayerl, Richter, & Sackmann, 1992).
Our calculations show distinct diffusion coefficient patterns for saturated and unsaturated glycerophospholipids, with the unsaturated glycerophospholipids exhibiting Figure 6. MSD of the lipid parallel to the lipid surface (x and y dimensions) vs. the simulation time, for the last 4 ns of (a) the DOPC and (b) POPC MD simulation. A linear fit for the part of the plot that corresponds to the diffusion regime in both systems is also displayed. higher values for this quantity. This observation is in agreement with the expected behaviour for these two classes and to the fluidity of the corresponding systems.

Computed deuterium order parameters (S CD )
Order parameters are an important property for the evaluation of bilayer ensembles (Vermeer, de Groot, Réat, Milon, & Czaplicki, 2007;Leftin & Brown, 2011). The computed order parameters (S CD ) were determined by following the Equation (2), and using an extension to the ptraj program written by Hannes Loeffler. In Equation (2), h represents the angle between the C-H vector of an acyl chain carbon (sn-1 or sn-2) and the bilayer normal vector (in this case, the z axis); and the brackets represent an ensemble average in the present cases over the last 20 ns of the MD simulations and over all the lipids.
The higher the S CD values, the higher the ordering and rigidity of the acyl chains. Lower values characterize the opposite behaviour. Results for DOPC, POPC and DMPC systems are presented in Figure 7 and are compared with experimental available data. As we can see, for the DOPC system, there is a great proximity of the computed S CD parameters to the experimentally determined S CD values, especially for the unsaturated bond between carbon atoms numbered 9 and 10. This is consistent with the good agreement obtained for structural and dynamic properties for this specific system. Considering the POPC system, there is a greater consistency when near the unsaturated region (Carbon atoms 9 and 10), contrasting with overestimation of order parameters for the regions immediately above and below the unsaturated portion. We also observe an overestimation of the S CD values considering the sn-1 chain. Focusing on the saturated phospholipids, we consistently observe an overestimation of S CD values when comparing with structural data; as depicted in example in Figure 7, for the DMPC situation. These general tendencies are in agreement with previous results by other authors when using GAFF in NPT lipid bilayer simulation (Jójárt & Martinek, 2007;Rosso & Gould, 2008;Siu et al., 2008;Skjevik et al., 2012).
If we compare the order parameters from the analogous PC and PE systems (data not showed), we can also gather that the PE systems have consistently higher S CD values for the two acyl chains (sn-1 and sn-2). This correlates with the structural information for these systems, since in PE lipids there is an higher packing efficiency due to the presence of hydrogen bonds between headgroups, provided by the single amine group rather than the trimethylammonium functional group of PC systems (Leftin & Brown, 2011). This translates into a more ordered bilayer and consequently into higher S CD values.

RESP charges and comparison of structural properties with MD data
In order to compare the RESP charges obtained with the protocol used in the study with similar ones described in the literature for DOPC and DMPC (Rosso & Gould, 2008;Siu et al., 2008;Skjevik et al., 2012;Tessier et al., 2008), a detailed comparison of the charge values reported for the different chemical groups ( Figure S3) was performed. A comparison with Jämbeck & Lyubartsev (Jämbeck & Lyubartsev, 2012a) charges for the PC phospholipids was also performed. These studies are based in a RESP charge methodology protocol, consistent with the GAFF. Results are presented in Table 4.
The results presented in Table 4 show, in general, a great similarity between the values obtained in the different schemes, with a near + 1 charge for the ethanolamine/ choline group and a À1 charge for the phosphate group, consistent with physiological pH protonation states. The charge of the acyl chains is near a value of 0 in most of the cases. With this comparison, we can see that our RESP charge's determination protocol leads to values very similar to other protocols and consistent between themselves.
When comparing the structural properties of volume per lipid, area per lipid and bilayer thickness with the work conducted by Rosso and Gould for the DOPC and DMPC phospholipids (Rosso & Gould, 2008), in the case of the DOPC system we have improved the approximation to the experimental data by 1.7% (from 2.4 to 0.7%) in the case of the volume per lipid, and we have obtained slightly better results when considering the area per lipid and the bilayer thickness; as for the DMPC system the accuracy to the experimental values obtained by Rosso and Gould is slightly better. When considering the work shown by Siu et al. for the DOPC phospholipid (Siu et al., 2008), we have obtained better results for the area per lipid and the bilayer thickness. Again it is important to stress out the fact that the simulations were     Recently, an AMBER modular framework for the simulation of lipids (Skjevik et al., 2012), presently available in the new AMBER 12.0 program, has been developed. Still, it is not capable of modulating biological membranes in the NPT ensemble, requiring the application of a surface tension (at times very high) for the modulation of correct experimental available data.

Comparison of LJ parameters
A comparison with recently made available LJ parameters from the work of Jämbeck & Lyubartsev (Jämbeck & Lyubartsev, 2012a) and also with Tessier et al. (Tessier et al., 2008) was also conducted, with the results presented in Table 5.
As depicted in Table 5, GLYCAM06 force field (Tessier et al., 2008) employs the same LJ parameters than GAFF, considering the equivalent atom types (CG and HC). However, LJ parameters are quite different when compared with the work performed by Jämbeck & Lyubartsev (Jämbeck & Lyubartsev, 2012a), mainly for the carbon atoms of the chemical groups -CH 2 and CH 3 .

Conclusions
With this work, we determined parameters for 12 lipid types consistent with the GAFF that enable to carry out molecular dynamics simulations on 12 different phospholipidic bilayer systems. Previous parameterization and performance evaluations were made, considering only one or two phospholipids. We have addressed the accuracy of the parameters based on the ability to reproduce differently experimental structural parameters: (i) volume per lipid, (ii) area per lipid, (iii) bilayer thickness and (iv) order parameters for the acyl chains. A consistent agreement with experimental data was obtained, for some of the studied systems, and when considering the volume per lipid property. For the area per lipid and bilayer thickness there was, respectively, a moderate underestimation at times and overestimation of the results at times too, consistent with other MD studies in the literature. In consistency with the previous observations, we also observe moderate overestimation of order parameters for some of the acyl chains. Some of the latter have been improved, as exemplified in some studies, by the application of an artificial surface tension to the bilayer system (e.g. with the use of an NPγT or NPAT ensemble). However, these alternatives impose significant constraints in the model system, compromising the structural knowledge that we would like to obtain out of the calculations and hampering the insertion of molecules/proteins to the membranes. Nevertheless our results show the potential for the simulations performed in a tensionless NPT ensemble, a feature that could favour future studies involving proteins inserted in biomembranes and the interaction of drug-like molecules with biomembranes.

Supplementary material
The supplementary material for this paper is available online at http://dx.doi.10.1080/07391102.2012.750250.