A Coarse-Grained Molecular Dynamics Study of DLPC , DMPC , DPPC , and DSPC Mixtures in Aqueous Solution

1 Department of Chemistry, Faculty of Science, Universiti Putra Malaysia, 43400 Serdang, Selangor, Malaysia 2 Enzyme and Microbial Technology Research Centre, Faculty of Biotechnology and Biomolecular Sciences, Universiti Putra Malaysia, 43400 Serdang, Selangor, Malaysia 3 Laboratory of Molecular Biomedicine, Institute of Bioscience, Universiti Putra Malaysia, 43400 UPM Serdang, Selangor, Malaysia


Introduction
Phospholipid bilayers are vital components of life and play a crucial role in membrane-mediated cell signaling [1]. A typical class of lipids is the phosphatidylcholines (PCs), which contains a family of saturated and symmetric phospholipids with the same molecular structures and different alkyl chain lengths. The alkyl chain lengths of PCs range from 12 carbon atoms in DLPC (1,2-dilauroyl-sn-glycero-3-phosphocholine) to 14 carbon atoms in DMPC (1,2-dimyristoylsn-glycero-3-phosphocholine), 16 carbon atoms in DPPC (1,2-dipalmitoyl-sn-glycero-3-phosphocholine), and 18 carbon atoms in DSPC (1,2-distearoyl-sn-glycero-3-phosphocholine) [2]. Phosphatidylcholines are considered as zwitterionic surfactant molecules due to the existence of two charged groups with different sizes. The anionic part involves phosphates, whereas the positive charge is located on the ammonium group [3]. Lipid molecules can exhibit micellar shapes, rod-like structures, and bilayers in aqueous solution due to their surfactant-like properties. The shape of the structure is usually dependent on the solution's concentration, temperature, and physicochemical properties [4].
The self-assembly of PC molecules is important for nanoemulsion formulations [5,6]. Both theoretical and experimental studies have demonstrated that nanoemulsions have great potential for application as drug carriers for transdermal drug delivery [7]. Phosphatidylcholines can spontaneously organize into bilayers due to their cylindrical shapes [8]. Figure 1 shows the molecular structures of DLPC (a), DMPC (b), DPPC (c), and DSPC (d).
Traditional molecular dynamics (MD), which applies all-atom level or united-atom level force fields, normally requires a long time to equilibrate a real physical system. Therefore, a coarse-grained molecular dynamics (CG-MD) technique was developed during the last several years as a simplified model for performing molecular dynamics. The use of coarse-grained (CG) models in exploring a variety of structural and dynamic properties in large molecular systems 2 Journal of Chemistry has provided valuable results when probing time and length scales beyond what is feasible with traditional MD [9]. The CG-MD technique can provide a good approximation based on the selection of the interaction parameters, as well as filling the gap between theoretical and experimentally determined results. Although CG-MD has provided valuable information for understanding the phase behavior, there is still limited information available for a theoretical understanding of surfactants and lipids. DMPC lipid bilayers were first modeled using the Klein model [10]. The parameters of this model were chosen to reproduce the correct density of the real physical system in good agreement with hydrodynamics data without considering the electrostatic properties of the water model applied. Marrink et al. [11] developed a CG model for simulating lipids and surfactants with better efficiency, flexibility, and simplicity by coarse-graining the molecules. Based on their results, the structural properties (e.g., the area per headgroup and the phosphate-phosphate distance), the elastic properties (e.g., the bending modulus and the area compressibility), and the dynamic properties (e.g., the lipid lateral diffusion coefficient and the water permeation rate) were consistent with the experimentally obtained values. Faller and Marrink [12] found that a binary mixture of DLPC and DSPC using CG model could produce a bilayer state. The phase transition obtained was characterized by the area per lipid headgroup, the order parameter, and the dynamics. Their results indicated the phase separation into a gel and a liquid state in a small temperature window at concentrations, with a majority of the longer lipid having a strong dynamic heterogeneity. Marrink et al. [13] simulated the coarse-grained type of transformation between a gel and a fluid phase in DPPC bilayers by cooling bilayer patches, whereas the reverse process was observed by heating the gel phase. These researchers reported that their results are highly dependent on both the system size and the simulation time. Their observations were in agreement with experimental findings and showed that the nature of the ordered phase obtained using the CG model was a gel phase, rather than a crystalline phase. In this study, we present a CG-MD model for the mixture of DLPC, DMPC, DPPC, and DSPC in aqueous solution for a more detailed understanding of lipid mixture behavior in the bilayer state.

Methodology
The CG-MD model systems were prepared using MARTINI model version 2.0 [9]. The MARTINI force field employs a 4 : 1 mapping scheme for nonring structures, where on average 4 heavy atoms are represented by 1 interaction center, each with an effective size of 0.47 nm and a mass of 72 amu. The MARTINI model for the coarse-grained DPPC consisted of 12 beads where the choline (NC3) group was represented by type Q0 (blue), the phosphate group (PO4) was represented by type Qa (yellow), and the glycerol group (GL1, GL2) and the hydrocarbon tail were represented by type Na (magenta) and C1 (cyan), respectively. DLPC, DSPC, and DMPC had similar types of beads, except that the number of C tail beads was different due to the difference in the chain length ( Figure 2).

Simulation Parameters.
The model systems were composed of 32 DPPC, 32 DMPC, 32 DSPC, and 32 DLPC molecules in 1500 coarse grained water beads (6000 real water molecules). The MD production simulation was performed for 500 ns in an NPT ensemble using Gromacs 4.0 software package [14]. Periodic boundary conditions were applied with a nonbonded cut off of cut = 1.2 nm. Bonded interactions were defined by harmonic potentials, whilst nonbonded interactions were defined by the Lennard-Jones potential. The Lennard-Jones potential was shifted from shift = 0.9nm to cut , and the Coulombic potential was shifted from zero to cut . The DPPC and water particles were independently coupled to a 298 K and 323 K bath with a relaxation time of = 1.0 ps. The pressure was maintained semi-isotropically at 1 bar using a Berendsen thermostat [15] with a coupling time of = 5.0 ps and an isothermal compressibility of 3.0 × 10 −4 bar −1 . The equations of motion were integrated using the leap frog algorithm [16] with a timestep of 30 fs. The box sizes for all models were chosen as 8.0 × 8.0 × 8.0 nm 3 . The CG potentials are much smoother than the atomic potential, thus, by applying coarse-graining, the effective sampling time could increase four times faster than the actual simulation time. This idea has been confirmed by the comparison of diffusion constants for atomic and CG water molecules [9].

Results and Discussion
The lipid bilayer properties, such as area per lipid, might not be calculated accurately using experimental techniques, such as -ray and neutron scattering, whereas some properties, such as the order of the hydrocarbon chain and the structure of the lipid bilayer, can be estimated accurately [17,18]. Here we report the results of the CG-MD simulation calculations of some important properties of a mixture of DLPC, DMPC, DPPC, and DSPC, such as the area per lipid, bilayer thickness, and lateral diffusion coefficients at 298 K, for a better understanding of the lipid bilayer structure in atomic detail. The CG-MD simulations results are highly dependent on the force field parameters selected and the temperature used.

Bilayer Formation. Figures 3(a)-3(e) illustrate the snap-
shot pictures of the mixed bilayer model system under semiisotropic pressure at 298 K at (a) 1 ns, (b) 6.7 ns, (c) 10 ns (d) 12.4 ns, and (e) 500 ns. The self-assembly of PC molecules was visualized using VMD [19]. A similar bilayer formation pattern was obtained in all lipid molecules. The formation of the bilayer started from a random configuration. As shown in Figure 3(d), the bilayer was completely structured at 12.4 ns, with the head groups facing towards the water molecules and the tails inside the bilayer. The lamellar shape was retained until the end of simulation for 500 ns. Figure 3(f) shows the location of different PCs, including DPPC (blue), DMPC (purple), DSPC (cyan), and DLPC (yellow), in the bilayer leaflets. They were distributed evenly in the bilayer. The assembly of lipid molecules into bilayers originates from two key factors: the hydrophobic interactions between the hydrocarbon tails, which induce the molecules to congregate, and the hydrophilic nature of the headgroups, thus, the headgroups will remain in contact with water. There is a competition between the hydrophobic and hydrophilic effects. One type tends to decrease and the other to increase the headgroup area in contact with water. The overall interactions give the concept of an "optimal surface area" per headgroup with a minimum total free energy per lipid molecule. This minimum energy contributes to the self-assembly of the bilayer [20].

Area per Lipid.
The area per lipid is one of the important physical properties of a lipid bilayer structure that determines the lateral organization of molecules in a bilayer system consisting of a single lipid molecule, as well as a homogeneous mixture of lipids [21]. The average area per lipid can be calculated from the equation = / , where is the total area of the dimension of the simulation box and is the number of lipid molecules in each leaflet [22]. was 64 in our model system. Therefore, the area per lipid was 6.48 nm in our coarse-grained bilayer mixture. The area per lipid obtained from our simulation was higher than the value of single DMPC, DLPC, or DSPC at 298 K due to the existence of different PC molecules with different chain lengths. Petrache et al. [23] performed experimental 2 H-NMR measurements for DLPC, DMPC, DPPC, and DSPC on randomly oriented, fully hydrated multilamellar samples. They found that at any given temperature, the area decreased with increasing acyl length, indicating increased van der Waals attraction for longer lipid chains. Kučerka et al. [24] reported several structural parameters, such as area/lipid, bilayer thickness, and hydrocarbon chain thickness, relating   to fluid phase bilayers of three different classes of PC bilayers at different temperatures, which were determined through the simultaneous analysis of small angle neutron and X-ray scattering data. Table 1 shows a summary of their experimental results. However, this could be a result of the calculation of the area per lipid projected onto the plane. This type of calculation may underestimate the actual area due to undulations that might occur in the bilayer model, especially in heterogeneous systems [25].

Bilayer Thickness.
Thickness of a lipid bilayer can strongly affect the properties of transmembrane proteins because the activity of membrane proteins, as well as the insertion and orientation of foreign particles in the membrane, are critically dependent on the bilayer thickness [26]. Some amount of water molecules can intercalate into the lipid headgroup region. Thus, the thickness of the bilayer can be demonstrated by plotting the density profile based on the calculation of the distance between two peaks [27]. Figure 4 shows the density profile of the headgroups of our mixed bilayer model. We determined the bilayer thickness from the peak-to-peak distance in the total profile. The value was 3.93 nm, which was in the experimental bilayer thickness range of the single DLPC, DMPC, DPPC, and DSPC values calculated experimentally. Petrache et al. [23] reported that hydrocarbons with longer chains might induce a greater thickness. Additionally, in the study of a binary mixture of DSPC-DLPC, Faller and Marrink  Figure 5: The MSD results of the mixed bilayer system over the time of simulation at 298 K. [12] showed that at the phase transition, the bilayer became thicker. They found that the bilayer thickness increased from 4.37 nm at 285 K to 5.04 nm at 265 K. From our results, the thickness of the mixture was greater than single DLPC and DMPC bilayers at 298 K but lower than single DPPC and DSPC at 323 K and 333 K [28]. Our results suggest that the existence of lipids with long chains in the mixed model can increase the bilayer thickness at room temperature.

Mean Square Displacement (MSD)
. MSD is used for calculating the value of the diffusion coefficient, which can then be used for exploring the lateral diffusion on interfaces. Figure 5 shows the MSD fluctuation results of DLPC, DMPC, DPPC, and DSPC, respectively, in the bilayer. From our results, all lipids in the mixed model diffused in a similar way. Hofsäß et al. [29] studied microscopically the interactions between cholesterol and lipids in biological membranes. They found that the diffusion of the DPPC molecules was slower than that of the cholesterol. Thus, longer phospholipid fatty acid chains with entanglement possibilities could slow down the diffusion. This conclusion is consistent with our MSD analysis. The values of the diffusion coefficients were 5.62 (±0.7) × 10 −7 cm 2 s −1 , 4.85 (±0.1) × 10 −7 cm 2 s −1 , 6.15 (±0.0) × 10 −7 cm 2 s −1 , and 4.53 (±0.4) × 10 −7 cm 2 s −1 for DLPC, DMPC, DPPC, and DSPC, respectively. It might be expected that the diffusion coefficients would decrease as the alkyl chain increases because the molecular weight and chain length could cause entanglements. Our CG-MD simulations results followed the same reported pattern, except for the DPPC molecule. This was a result of the DPPC dipole effect, which is smaller than the dipole values of the other PC molecules.

Conclusions
The self-assembly of a mixed bilayer system, consisting of 32 DLPC, 32 DMPC, 32 DSPC, and 32 DSPC, at room temperature was successfully explored using a CG-MD simulation 6 Journal of Chemistry technique. From our results, the mixing of phospholipids with similar heads but different acyl chain lengths can significantly affect several bilayer properties, such as the area per lipid, the bilayer thickness, and the lateral diffusion. This type of theoretical study is important for achieving a detailed understanding of the membrane signaling environment. At the same time, the formulation of nanoemulsions as drug carriers may be improved to maximize drug permeability through the skin barrier.