Study of Lipid Heterogeneity on Bilayer Membranes Using Molecular Dynamics Simulations

Human cell membranes consist of various lipids that are essential for their structure and function. Several experimental techniques have been used to characterize the composition of human cell membranes; however, it is challenging task to depict in theoretical calculations. In this work, we have investigated the structure-function relationship of lipids in both homogeneous and heterogeneous bilayer models using Molecular-Dynamics to delineate the effect of heterogeneity on the biophysical properties of membranes. Results illustrated that the biophysical properties of heterogeneous membranes are dependent on the lipid composition and concentration. We observed that the presence of cholesterol in combination with other lipids, introduced compactness of the membrane, increasing the membrane thickness. The density of lipid head group, phosphate , and glycerol-ester in presence of cholesterol with or without sphingomyelin is an underlying reason for membrane ordering. The radial distribution function shows that the cholesterol, sphingomyelin and phosphatidylethanolamine self-interaction and the interaction between cholesterol and phosphatidylethanolamine determine the structure and function of the heterogeneous membrane. Our ndings provide a baseline for membrane heterogeneity that would help in understanding the physiological properties of membranes and may help to wisely select the heterogeneous bilayer model to mimic the realistic human cell membranes and the associated phenomenon.

study the uctuation among nanoscale assemblies of lipids [23,24]. The permeability barrier of the bilayer membrane impacts the treatment of many diseases in terms of challenges in drug discovery and the threat of drug resistance [25,26]. The molecular dynamics simulation has become an indispensable tool to comprehend the driving forces which leads to the structural and dynamics of the membrane [27].
Initially, simulation studies were carried out on the bilayer membrane containing a single type of lipid. The rst simulation study was done for the dipalmitoylphosphatidylcholine (DPPC) membrane [28,29]. The bilayer model containing a single type of lipid was very far from the real biological membrane. Later, substantial simulation studies were carried out on the heterogeneous membrane considering two or more lipid types [30][31][32]. Recently, various computational studies were performed on the heterogeneous bilayer membrane containing one or two lipids to delineate the phase separation [18,27,33], membrane permeability [19,34] and condensing effects [35,36]. The liquid-ordered and -disordered, and coexistence of the two phases by considering heterogeneous composition was depicted by Sodt et al [37]. The changes in the biophysical properties of the membranes were obtained depending on the composition as well as the concentration of lipid and provided valuable insights into the structural properties of mixed type bilayer model [27]. However, the effect of heterogeneity containing more than two lipids on the biophysical properties of bilayer membranes are poorly understood and the available information is just a tip of the iceberg.
In this study, various atomistic bilayer models i.e. single, binary, ternary, quaternary, and senary types have been constructed. The composition of bilayer models has been chosen in such a way that the effect of lipids can be distinctively explained. We have investigated the structure-function relationship of the  [43,44] for the lipids with TIP3P water model [45]. The protocol de ned by CHARMM-GUI has been used to prepare the simulation system [38][39][40], which includes steepest descent energy minimization and six-step of equilibration in which lipids restraints gradually decrease from 1000 kJ/mol×nm 2 . Twostep of NVT equilibration i.e. 2 ns constraint equilibration and 1 ns of relaxed NVT equilibration has been performed. The NVT equilibrated structure is subjected to 5 ns NPT equilibration in three constitutive phases i.e. 1ns constraint NPT equilibration then 4 ns NPT relax-equilibration. Finally, a total of 1.3 µs production run was carried out (100 ns for each bilayer models). The time step is 2 fs and saved the coordinates every 2 ps. The Nosé-Hoover thermostat (with a reference temperature of 310.15 K and a coupling constant of 0.2 ps) together with the semi-isotropic Parrinello-Rahman barostat (with a reference pressure of 1 atmosphere, the compressibility of 4.5 × 10 −5 bar −1 and coupling constant of 1 ps) is used [46,47]. H-bond of the lipid molecules is constrained using the LINCS algorithm [48]. The periodic boundary condition is used in all Cartesian directions. A switch function starting at 10 Å is used for the van der Waals (vdW) interactions with a cutoff of 12 Å. Particle-Mesh Ewald (PME) algorithm [49,50] is used to calculate the long-range electrostatic interactions where the cut-off of short-range electrostatic interactions is applied at 12 Å.
Membrane trajectories are analysed using Gromacs utilities and VMD software [51].

Results And Discussion
Surface area per lipid The surface area (SA) per lipid gives insights into the packing of the membrane, where a low and high value of SA usually indicates a membrane being more rigid and more uid respectively. In addition, it is commonly used to verify the system is equilibrated. The SA graph observed to be remained constant and there was no major drift throughout the simulation, as shown in Fig. 2. A snapshot of the representative structures of the bilayer model types from 80 to 100 ns frame is shown in Fig. 3

Deuterium order parameters
The order parameter (S CD ) was calculated to quantify the orientation and exibility of the lipid tails.
Generally, carbon atoms adjacent to the lipid head group have a higher S CD value and decreases down with the length of the carbon chain. This is because the lipid head group restricts the motion whereas the tail regions move relatively free in the bilayer. The S CD values of both sn-1 and sn-2 chains of POPC, POPE, POPI, POPS, and PSM were higher in the case of heterogeneous bilayer models than the pure POPC bilayer models, as shown in Fig. 4, which correlated (inversely proportional) with the lower SA ( Table 2). The S CD of POPC increased and more ordered whereas S CD of POPE dropped down when added into the POPC bilayer model (M3). A double bond between C9 and C10 of sn-2 chain of POPC, POPE, POPI, and POPS restricts the rotation of the fatty acid chain that resultant to decrease in S CD of carbon number 6 to 9. The obtained S CD for pure POPC and pure POPE bilayer models demonstrated an excellent agreement with the NMR experiment [58,59]. Therefore, we expected a similar agreement for the heterogeneous bilayer models. The S CD of both sn-1 and sn-2 chain of POPC, POPE, POPI and POPS were higher in CHL and PSM/CHL containing bilayers (M7, M9-M13). The introduction of PSM to the POPC/POPE bilayer seemed not to be noteworthy. The average S CD value of PSM in the POPC/POPE/PSM bilayer was 0.28 which could be correlated with the experimentally obtained value [60].
Addition of PSM along with CHL into the POPC/POPE bilayer improved the ordering of lipid chain signi cantly, even more than the CHL containing POPC/POPE bilayer. This suggests that the ordering of lipid chain improves the SA of PSM/CHL containing bilayer (M9), as shown in Fig. 4 and Table 2. On the contrary, the addition of POPI and POPS into the POPC/POPE/PSM/CHL bilayers caused an adequate drop in the ordering of lipid chain, but as the concentration of PSM and CHL increased (M10-M13), lipid chain got more ordered (Fig. 4) by letting the lipids alined better and packed closer, as indicated by the lowering of SA of M13 (Table 2).

Mass density and membrane thickness
The mass density pro le (MDP) is an analysis of how the lipid changes from bulk water solution to the centre of the bilayer. The calculated density pro le of lipids, water, lipid tails, head groups, P-atom of head group, and glycerol ester along the z-axis is shown in Fig. 5 and Supplementary Fig. S1 online. The section with the lowest density in the pro le was the centre of bilayer structure, which suggested the ends of lipid tails. The distance between two density peaks of symmetric lipids indicated the membrane thickness. The general position and the shape of the density pro le of pure POPC and pure POPE bilayer models were nearly identical, except for the density amplitude, where the density amplitude of POPE was higher than POPC, as shown in Fig. 5 (a). The mass density pro le of the lipid components, i.e. water, lipid tails, head groups, P-atom of head group, and glycerol ester exhibited that the higher density of POPE arose from PO 4 and glycerol ester group, as the density of PO 4 and glycerol-ester was observed slightly higher in POPE than POPC (see Supplementary Fig. S1 online). The mass density of binary type bilayer, i.e. M3 was higher than pure POPC (M2) and lower than pure POPE (M2) bilayers. This might be because of the increase of head density [see Supplementary Fig. S1 (a) online] and ordering of POPC (Fig. 4). The total density of trimer type of bilayer models, i.e. M4, M5, M6, and M7, remained nearly identical whereas, the density of POPE and head group and glycerol-ester of POPC increased, and the density of tails overlapped. It signi es that lipid tails occupy the same space in absence of CHL [see Supplementary Fig.  S1 (b) online]. The lipid density increased towards the exterior and decreased in the bilayer centre in the presence of CHL due to the shrinking of lipid tails, as suggested by S CD parameters (Fig. 4). In the case of the quaternary type of bilayers, the density curve of POPC/POPE/POPI/POPS and POPC/POPE/PSM/CHL remained distinct from each other, where the lipids and their components density increased in presence of PSM and CHL. Moreover, the density of the lipid components showed that lipid tails overlapped with each other except for PSM tail, which corresponded to the positioning of lipid tails in the same space. The density of PSM tails seemed to shift away from the centre and was similar to the density of glycerol-ester, as shown in Supplementary Fig. S1(c). The total mass density of the senary type of bilayers, i.e. M10-M13 seemed to be nearly identical. The amplitude of density pro le was minutely higher in the M13 than other senary bilayers, which might be because of higher the concentration of PSM, and CHL [ Fig. 5 (a)]. Further, the density of lipid components exhibited that the density of lipid head groups increased with the relative concentration of PSM and CHL, and therefore, the total density of M13 was minutely increased [see Supplementary Fig. S1 (d) online] However, the density of PSM and CHL containing bilayers suggested that CHL always remained near and below the glycerol-ester group of lipids, indicating that the increased in the lipid order parameters might be because of inserting the hydroxyl in the free space. A cursory view of the calculated bilayer thickness (  Table 3. Thus, the results suggested that ordering of lipid tails and the increment of lipid headgroup, PO 4, and glycerol ester density in the presence of PSM and CHL with signi cant concentration were the underlying reason for the lowest SA and ordering of bilayer membrane.

Radial distribution functions
The radial distribution function (RDFs) is calculated to get further insights into the local packing and lateral ordering of lipid bilayers and to explore the membrane heterogeneity. The 2D RDFs for the phosphates and carbonyl oxygen of phospholipid, viz. POPC, POPE, POPI, POPS, and PSM, and the oxygen atom of CHL were calculated for all the possible combinations and shown in Fig. 6 and Supplementary Fig. S2. In the 2D RDFs plot, the rst peak with higher g(r) suggested the probable distance between phosphates and the oxygen of phosphate groups, which was pointing toward each other. The second peak with lower g(r) corresponded to the conformation when the oxygen atoms were pointing away from each other. The RDFs of phosphate-phosphate in all the combinations of lipids showed two sharp peaks except for POPC/POPE, which seemed to be exhibiting two different conformations of phosphate oxygen. The RDFs of POPC self-interaction in both pure POPC and mixed POPC/POPE bilayers showed two major peaks with a consistent shape. However, POPE self-interaction in POPC/POPE bilayer exhibited greater peak amplitudes than pure POPE bilayer model, suggesting the dominance of POPE self-interaction in POPC/POPE bilayer [ Fig. 6 (a)], and the physical interaction is shown in Fig. 7. In the case of trimer bilayers, two major peaks were obtained except for POPC/POPE where a larger peak amplitude was obtained for POPE/POPE and POPC/POPE in the CHL containing bilayer (M7). It revealed a greater association of POPE with POPE and POPC in the presence of CHL, as shown in Fig. 6 (b), where g(r) values of POPC/CHL and POPE/CHL indicated a strong interaction between CHL oxygen and the carbonyl oxygen of POPC [ Fig. 6 (e) and Supplementary Fig. S4 online]. In case of ternary systems, i.e. M8 and M9, the calculated RDFs exhibited a strong interaction of POPC with POPC and POPE in the presence of PSM and CHL as shown in Fig. 6 (c), where the CHL seemed to be bonded strongly with POPE than POPC [ Fig. 6 (f)]. Comparing the RDFs of M7 and M9, it could be seen that the g(r) value of POPC/POPC and CHL/POPE (Fig. 6) increased signi cantly in the presence of PSM (M9), where CHL seemed to be favouring POPE a little over POPC in the presence of PSM [ Fig. 6 (f) and Supplementary Fig. S5 online], suggesting a reason for the decreased SA of M9 compared to M7. The RDFs of quaternary type bilayers indicated that the association between POPE/POPE and POPC/POPS increased with the increasing concentration of PSM and CHL into the bilayers. Notably, CHL in quaternary type bilayers tend to be associated poorly with POPE in low concentration of CHL and PSM, which was re ected in the RDFs that were usually below the values observed in the CHL/POPC bilayer. However, CHL associated strongly with the POPE in the high concentration of CHL and PSM, as suggested by the g(r) value, [Fig. 6 (f) and Supplementary Fig. S3

Conclusions
The present MD simulation study explains the effect of lipid heterogeneity on the biophysical properties of bilayer membranes containing 1 to 6 lipids (i.e. POPC, POPE, POPI, POPS, PSM and CHL). It showed that the compactness and ordering of the membranes are highly dependent on the degree of heterogeneity. The calculated SA and order parameters for the pure POPC and POPE bilayer models showed an excellent agreement with the experimental results. The presence of CHL with or without PSM in a signi cant concentration and in combination with other lipids resulted in a more compact and ordered bilayer membrane. The introduction of POPI, POPS and PSM into the POPC/POPE bilayers did not signi cantly improve the ordering of lipid tails and compactness of bilayer membrane. The order parameters, mass density, membrane thickness and lipid interactions in the studied bilayer models were dictated by the lipid composition and concentration. These parameters were higher for the bilayers with a higher concentration of CHL and/or combination of CHL and PSM. The study signi es the role of CHL self-interaction, CHL/PSM, CHL/POPC, CHL/POPE and POPE/POPE interaction, and dominance of CHL and POPE self-interaction, and CHL/POPE interaction in lowering the SA, and improving the compactness of the membrane. Since it is a common practice to use a simple one or two lipid bilayer in membraneassociated studies, and therefore, the effect of heterogeneity on the protein dynamics, protein-drug interactions across membranes, drug permeability, and membrane-nanoparticle interactions could not have been explained. Therefore, the membrane composition should be wisely considered during the simulation, when investigating the function of the protein, drug and nanoparticle that reside or interact with the membrane. The CHL/PSM or CHL rich heterogeneous bilayer may act as a representative computational model of the human cell membrane and can be used for further studies.

Declarations
Data availability All data generated or analysed during this study are included in this published article (and its Supplementary Information le).
Author contributions G.N.S conceived and supervised the study. N.K. performed modelling, simulations and data analysis, and interpreted the results. Writing of the manuscript is done by both authors.

Competing interests
The authors declare no competing interests.   The surface area per lipid of all the considered bilayer models at 310 K.