Effects of lipid heterogeneity on model human brain lipid membranes.

Cell membranes naturally contain a heterogeneous lipid distribution. However, homogeneous bilayers are commonly preferred and utilised in computer simulations due to their relative simplicity, and the availability of lipid force field parameters. Recently, experimental lipidomics data for the human brain cell membranes under healthy and Alzheimer's disease (AD) conditions were investigated, since disruption to the lipid composition has been implicated in neurodegenerative disorders, including AD [R. B. Chan et al., J. Biol. Chem., 2012, 287, 2678-2688]. In order to observe the effects of lipid complexity on the various bilayer properties, molecular dynamics simulations were used to study four membranes with increasing heterogeneity: a pure POPC membrane, a POPC and cholesterol membrane in a 1 : 1 ratio (POPC-CHOL), and to our knowledge, the first realistic models of a healthy brain membrane and an Alzheimer's diseased brain membrane. Numerous structural, interfacial, and dynamical properties, including the area per lipid, interdigitation, dipole potential, and lateral diffusion of the two simple models, POPC and POPC-CHOL, were analysed and compared to those of the complex brain models consisting of 27 lipid components. As the membranes gain heterogeneity, a number of alterations were found in the structural and dynamical properties, and more significant differences were observed in the lateral diffusion. Additionally, we observed snorkeling behaviour of the lipid tails that may play a role in the permeation of small molecules across biological membranes. In this work, atomistic description of realistic brain membrane models is provided, which can add insight towards the permeability and transport pathways of small molecules across these membrane barriers.


Introduction
The cellular membrane is an essential component of every cell, forming barriers within the cell itself and with the external environment. These membranes are composed of a complex mixture of lipids and proteins that typically contain hundreds of different lipid species, which differ in headgroup and/or hydrocarbon tail length and degrees of saturation. Cell membranes have also been found to have approximately 30% of their interfacial area covered by a variety of proteins. [1][2][3][4] Normally, the lipid composition varies significantly between the various organelles and subcellular compartments 2,[5][6][7][8] and between different cells. Additionally, the composition of the plasma membrane has been found to be dependent on the organism, cell stage, environmental factors, and tissue types. [9][10][11][12][13][14][15][16][17][18] Changes in lipid composition have also been linked to many diseases. [19][20][21][22][23] For example, depression, anxiety, and drug addiction are strongly associated with changes in the lipid composition in the brain. [24][25][26] Furthermore, various neurological diseases including Alzheimer's disease (AD), Parkinson's disease, and epilepsy have been linked to lipid imbalances that result from the deregulation of biochemical pathways associated with phosphoinositide (PIPs) lipids. [27][28][29] The complex and evolving composition of cell membranes and the fact that they play significant roles in both healthy and diseased cells have led to an increased emphasis to understand how changes in the lipid composition affect the structural and dynamical properties of the membranes. Sophisticated experimental techniques, including single-particle tracking, fluorescence correlation spectroscopy, super-resolved imaging, scattering, solid-state NMR, and mass spectrometry have recently been used to provide insight into the lateral organisation of the various components within membranes, as well as the potential causes for this organisation. [30][31][32][33][34][35][36] Despite advances in these various techniques, the full molecular-scale details of the interactions that result in the various structural and dynamical properties of interest are currently not accessible by experiment alone.
Computational modelling, in principle, is capable of providing this added level of detail and has become an essential tool in the investigations of lipid membranes. 37,38 In recent investigations of brain lipid membranes, a range of different models have been used to represent the membranes and their interactions with small molecules or peptides. The models used range in lipid composition from pure 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) [39][40][41][42] to mixtures of up to five different lipid species (including cholesterol) 43,44 to complex realistic lipid mixtures. 45,46 In the current work, all-atom classical molecular dynamics (MD) simulations were used to investigate the effects of increasing lipid complexity on the bilayer structure, lipid dynamics, and dipole potential. We used two simple membranes consisting of pure POPC and POPC-CHOL in a 1 : 1 ratio. Meanwhile, we present the first heterogeneous brain cell models consisting of 27 lipid components in an effort to mimic the experimental compositions identified by Chan et al. 47 (Fig. 1 and Table 1).

Model membranes
Two brain cell membranes representing the healthy (BH800) and Alzheimer's diseased (BAD800) state were modelled based on experimental lipidomics data of the prefrontal cortex cells 47 ( Fig. 1 and Table 1). BH800 and BAD800 each consists of 800 lipids per leaflet. We also conducted simulations of membranes with the same lipid composition as BH800 and BAD800 but with only 200 lipids per leaflet, BH and BAD, respectively, to investigate system size effects. Additionally, we simulated a pure 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) membrane, and a membrane consisting of a 50 : 50 mixture of POPC and cholesterol to investigate the effects of lipid heterogeneity on the membrane properties. Each of these simpler membranes contain 200 lipids per leaflet. Table 1 reports their symmetric lipid distributions in the upper and lower leaflets. The chemical structures of the various lipids used are shown in Fig. 2. All systems were generated using the CHARMM-GUI Membrane Builder, using the CHARMM36 parameters for lipids and the CHARMM TIP3P model for water. [48][49][50] As a result, each lipid molecule is randomly placed within the two leaflets of the bilayer. In doing so the structure of each lipid molecule is randomly chosen from a lipid structural library which contains 2000 conformations for each lipid type that have been produced from a combination of homogenous and heterogeneous lipid membrane simulations. Then each membrane was solvated with a water thickness of at least 22.5 Å (default setting) and neutralised to a salt concentration of 150 mM NaCl.

MD simulation protocol
GROMACS version 2018 was used to perform all MD simulations in this manuscript. 51 Initial structures were energy-minimised using steepest descent to remove steric clashes between atoms. Then, a series of simulations were carried out using the NVT (constant number of particles, volume, and temperature) and then the NPT (constant number of particles, pressure, and temperature) ensembles for at least 0.45 ns in order to equilibrate the temperature and density of the systems. In doing so, the prescribed simulation protocol by CHARMM-GUI was followed. 49,50 All equilibrated systems were simulated under NPT conditions, to a total length of 300 ns for smaller systems (POPC, POPC-CHOL, BH, and BAD), and to a total length of 200 ns for larger systems (BH800 and BAD800). In the production simulation, the temperature was maintained at 310 K by the Nosé-Hoover thermostat, 52 with a time constant of 5 ps. Semiisotropic pressure coupling was maintained at 1 bar by the Parrinello-Rahman barostat, 53 with a time constant of 5 ps and standard compressibility of 4.5 Â 10 À5 bar À1 . The Verlet cut-off scheme was employed. Electrostatic interactions were calculated using the Particle-Mesh Ewald algorithm. Both electrostatic and van der Waals interactions were cut off beyond 1.2 nm. All bonds involving hydrogen atoms were constrained using the LINCS algorithm. 54

Data analysis details
Ensembled averages were analysed over the full length of the trajectories. Standard errors were calculated as standard deviations of the mean. The results of the two different sized healthy (BH800 Fig. 1 (a) General compositions of the lipid types used in the BH800 and BAD800 brain cell models (see Table 1 for full detail of different membranes). (b) Side view of a complex brain membrane system. Lipids are coloured according to the pie charts. and BH) and diseased (BAD800 and BAD) membranes are generally similar. Data for the smaller membranes (BH and BAD) are mainly included in the ESI, † unless stated otherwise.
Area per lipid (APL) was calculated by Voronoi analysis in the membrane analysis tool, MEMBPLUGIN. 55 The oxygen atom (name O3) was used as the key atom to represent CHOL, while a triad of glycerol carbons were used to represent non-CHOL lipids (names C2, C21, C31 for PC, PE, PS, and PI; names C1F, C2S, C3S for SM and CER) ( Fig. S1a and b, ESI †). These key atoms were selected because they lie at a similar depth in the membrane. By projecting the coordinates of the key atoms onto the xy-plane, a series of Delaunay triangulations and Voronoi diagrams can be constructed 56 to give the APL (Fig. S1c, ESI †).
Bilayer thickness was defined as the average phosphate-tophosphate positions, and was calculated using MDAnalysis. 57 Surface roughness (R) was calculated using MDAnalysis. 57 The C2 and C2S atoms of all non-CHOL lipids were used as reference atoms. The method by Plesnar et al. 58 was applied: where n is total number of reference atoms, z i is the z-coordinate of reference atom i, and z m is the mean z-position of all reference atoms. Order parameter (S CH ) quantifies the time-averaged C-H bond angle, y, with respect to the bilayer normal. Using MEMBPLUGIN, 55 S CH was calculated according to: Tilt angles for the various lipid species were defined according to the vectors in Table 2. For the CER lipid group, only the ring tilts in the headgroup of the cerebrosides were considered. All angles were measured with respect to the bilayer normal, and computed using MEMBPLUGIN. 55 Lipid interdigitation was calculated using MEMBPLUGIN, which measures the degree of acyl chain interdigitation using a correlation-based fraction that gives the mass overlap of two leaflets (I r ) and its corresponding width of the overlapping regions (w r ). I r is measured so that 0 r I r r 1 indicates between zero to complete mass overlap. In addition, a coordination-based fraction (I C ) evaluates the number of heavy Table 1 Compositions of the membrane models (per leaflet) corresponding to the CHARMM-GUI lipid information. Complex healthy and diseased brain models, BH800/BH and BAD800/BAD, respectively, contain up to 27 lipid components, to include cholesterol (CHOL), phosphatidylcholines (PC), phosphatidylethanolamines (PE), phosphatidylserines (PS), phosphatidylinositols (PI), sphingomyelins (SM), and ceramides (CER) lipids. Lipid distributions are symmetric in the upper and lower leaflets atoms that are in contact with the opposing leaflet within the recommended default cutoff at 4 Å. 55 Lipid tail snorkeling was determined by the terminal carbon z-positions of each hydrocarbon sn chain. The absolute values of z from both leaflets were taken, and binned to a histogram in order to obtain a normalised probability density, using MDAnalysis. 57 Dipole potential (C(z)) was computed using the gmx potential module from the GROMACS package. 51 The simulation box was divided along the bilayer normal to obtain the charge density per slab. Double integration of the charge density, r(z), then yields the dipole potential, C(z): 59 where z 0 is the z-position at the bulk water centre at which the potential is set to zero. Lateral diffusion was calculated using the gmx msd module from the GROMACS package. 51 Diffusion coefficients (D) are calculated from the slopes by fitting a straight line (Dt + c) through the plots of the mean square displacements (MSDs), 51 between 30-200 ns for POPC, POPC-CHOL, BH, and BAD, and between 20-150 ns for BH800 and BAD800. However, the lipids in each membrane experience subdiffusion since productions were simulated on a sub-microsecond timescale. 60,61 Thus, the coefficients obtained in this work are merely estimated values of D, which we use because we are mainly interested in the qualitative comparisons of the lipid mobility between the various bilayer environments, as was similarly done in previous work. 62 3 Results and discussion

Effects of lipid heterogeneity on the structural properties of lipid membranes
In the homogeneous POPC membrane, the bilayer thickness and APL equilibrate to 38.8 Å (Fig. S4, ESI †) and 64.8 Å2 (Fig. S2, ESI †),  Table 1. Table 2 Summary of key membrane properties: (mean AE standard error). Area per lipid (APL). Bilayer thickness measured from the average phosphorusphosphorus distance between leaflets. Interdigitation parameters where I r is the interdigitated mass overlap, and w r is the corresponding width region of the mass overlap. I r is measured so that 0 r I r r 1 indicates zero to complete mass overlap. I C is the fraction of heavy atoms that are in contact with the opposing leaflet. Surface roughness parameters of the upper and lower leaflets. Tilt angles for the different lipid headgroups present in the various membranes, and the defining vectors used to determine this angle  58,63,64 The addition of 50% cholesterol in POPC-CHOL causes the bilayer to thicken to 45.8 Å (Fig. S4, ESI †), and the APL of PC lipids to reduce to 54.8 Å 2 (Fig. S2, ESI †), characteristic of the well-established condensing effect of cholesterol. 58,65,66 As the membrane gains lipid heterogeneity in BH800/BAD800, the thickness increases to 47 Å (Fig. S4, ESI †), but interestingly, the APL of PC lipids also increases slightly to 56 Å 2 (Fig. S2, ESI †) in comparison to the POPC-CHOL membrane (Table 2). However, PE and PS lipids do not seem to feel the condensing effect of cholesterol to the same extent as PC, since both lipid types maintain APLs (Table 2 and Fig. S2 Inversely correlated to the bilayer thickness is the lipid interdigitation, whereby the lipid chains in the upper and lower leaflets overlap at the bilayer centre. Such is the case where an increase in thickness from POPC to POPC-CHOL is accompanied by a decrease in the interdigitation parameters (Table 2 and Fig. S5, ESI †). Our findings are in line with the fact that cholesterol is proposed to reduce lipid interdigitation. 70 However, complex membranes BH800/BAD800 at similar cholesterol levels to POPC-CHOL depart from this trend by showing a greater extent of interdigitation even as the bilayer thickens (Table 2 and Fig. S5, ESI †). This likely arises from the membrane containing a range of lipids that are diverse in chain lengths as well as degrees of saturation. Changes in the interdigitation have been suggested to have numerous implications, including on the interfacial tension of lipid membranes, 71 the formation of membrane micro-domains, 72 and the modulation of lipid droplets. 73 The condensing effect of cholesterol is also closely related to its ordering effect on the lipid chains in a membrane, 66 which is observable using the order parameter (S CH ). Typically, lipids with a smaller APL will see an increase in the order of their hydrocarbon tails. At first, only the sn tails of POPC lipids in the various membranes were evaluated (Fig. S6, ESI †), and effectively displays a higher order in the POPC-CHOL, BH800, and BAD800 systems, relative to the pure POPC system. To approximate the order of the overall membrane, the order parameters were averaged across the various lipid types containing either a sn16:0, sn18:0, or sn18:1 tail (Fig. S7, ESI †), all of which show a similar trend to the order in Fig. S6 (ESI †). However, in BH800/ BAD800, the B5 Å 2 difference in APL between the SM and CER lipids and the PC, PE, PS, and PI lipids ( Table 2) appears to not have been resulted in a greater amount of order in the SM and CER lipid tails (Fig. S7c compared to Fig. S7a and b, ESI †). Instead, in the membranes that contain 50% cholesterol, the value of S CH reaches a maximum (S CH B 0.4) in the middle of the tails. Generally, the order profiles are consistent with existing literature that has reported studies of bilayers containing B50% cholesterol. [74][75][76][77] The values of the order parameter observed in the lipids within BH800/BAD800 indicate that brain cell membranes realistically adopt a liquid-ordered state, since usually S CH o 0.3 in liquid-disordered membranes.
The tilt angle of cholesterol is another quantity that further explains how cholesterol imparts order on the packing environment. The orientation of CHOL shows the narrowest distribution in BH800/BAD800 (Fig. S10, ESI †), with the smallest average tilt angle of 111 (Table 2). This value is smaller compared to those previously measured in experiment [78][79][80] and simulation 69,81,82 (B14-251), all of which studied membranes with lower percentages of cholesterol within them, but the tilt is consistent with the angle observed in a simulation of a lipid bilayer with a 1 : 1 ratio of POPC and cholesterol. 62 Otherwise, tilt angles that were measured in the headgroups of PC and PI lipids 69,83 (Table 2 and Fig. S10, ESI †) are similar to those in simple model bilayers (where comparable). The orientation of the CER sugar headgroups seem to correspond to the tilt angles reported for membranes containing glycolipids (B401), 84-86 even though they have a glycerol backbone (glyceroglycolipids) instead of the sphingosine backbone of the CER lipids used in this study (glycosphingolipids). The broad distributions in the tilt angles demonstrate a degree of freedom on the membrane surface that seems to be independent of the lipid packing. This is further evidenced by the similar hydration features of waters surrounding CHOL (Fig. S11, ESI †) and PC (Fig. S13, ESI †) lipids that are shared across the various membranes, which is unaffected by changes in the physical properties (Table 2). Previous simulation studies have found similar distances that define the first neighbor hydration shell around the oxygen atom in cholesterol, 62 the nitrogen atom in the headgroup of the PC lipids, 87,88 and the nitrogen atom in the headgroup of PE lipids. 89 (Further details on the hydration of lipids are discussed in the ESI. †) Incidentally, the steadiness in the APL, thickness, and roughness parameters (Fig. S3, S4 and S8, ESI †) show an effectiveness of the CHARMM force field to quickly stabilise not only simple membranes, including POPC and POPC-CHOL, but more importantly complex membranes like BH800 and BAD800 that contain a cocktail of up to 27 lipid components. The average surface roughness of POPC (1.1 Å) is higher than POPC-CHOL (0.7 Å) ( Table 2), whereby the relationship between the two is consistent with the findings by Plesnar et al. 58 As lipid tails were expected to interdigitate vertically, we found instances in the simulation where some tail ends would snorkel towards the headgroup interface as in Fig. 3. In existing literature, Feix et al. has demonstrated flexibility in the tails of nitroxides to infer vertical fluctuations of terminal methyl groups in dimyristoylphosphatidylcholine (DMPC) bilayers. 90 In addition, the snorkeling phenomena has been ascribed to the side chains of amino acids in transmembrane proteins to assist in protein-membrane interactions. [91][92][93] Meanwhile, a membrane-only simulation has observed the sn24 : 1 tail to adopt a hooked conformation, 94 but not quite bending to the extent of the U-shapes seen in Fig. 3. It is curious that lipid tails would snorkel deeply to an area of the membrane with high order (Fig. S7, ESI †). Thus, we evaluated the z-positions of the terminal carbons, and used a threshold of z 4 10 Å to indicate any snorkeling events. Similar trends are expressed in the hydrocarbon chains across all systems (Table 3 and Table S2, ESI †). Mainly, the tails tend to equilibrate around the bilayer centre (z o 10 Å) (Table 3  and Table S2, ESI †). Otherwise, the probability of snorkeling increases with increasing degrees of unsaturation in the hydrocarbon chain, for instance, in the order of sn16:0/sn18:0/sn20:0 o sn18:1 o sn18:2 o sn20:4 o sn22:6 in the case of BH800 (Table 3). Additionally, for tails with the same degree of unsaturation, we have found that the likelihood of snorkeling increases with the length of the hydrocarbon tail, for instance, sn24:1 4 sn18:1 (Table 3). Compared to POPC, it is also immediately obvious that the snorkeling events become much less probable in POPC-CHOL and BH800/ BAD800 (Table 3 and Table S2, ESI †), which also seems to correlate to a decreasing number of water permeation events (Table S3, ESI †).

Changes in transmembrane dipole potential resulting from lipid heterogeneity
The distribution of electric charge in a bilayer is important because it affects the mechanisms of interaction and permeability of drugs with the lipid bilayer and/or the surrounding membrane proteins. 95,96 Thus, the dipole potential profile of membranes is particularly crucial in the process of novel drug design. For example, the electrostatics of bacterial membranes is especially relevant for the design of novel antimicrobial drugs since the electrostatic interactions during drug binding is what provides the driving force for the insertion process of the drug into the membrane. 97,98 Poisson's equation, as we have used here, is typically used to calculate the dipole potential from MD trajectories. 63 Fig. 4 shows the dipole potential across each of the membranes we have simulated here. Generally, the dipole profiles across all systems have matching qualitative features. Initially, the potential starts at the reference value in the aqueous phase at 0 V. Then, the potential increases sharply over the lipid headgroups to experience potential barriers at the membrane interface at z B AE15 Å, before reaching a global maximum at the bilayer centre (Fig. 4). Between POPC and POPC-CHOL, the presence of cholesterol is demonstrated to raise the potential barriers most substantially at the centre by 0.3 V (Fig. 4). In the realistic membranes BH800/BAD800, lipid diversity in the headgroup as well as tail regions lower potential barriers at the interfaces to 0.42-0.45 V, and at the centre to 0.89-0.90 V (Fig. 4). Thus, the dipole potential is sensitive to the membrane composition and environment. And this reiterates that the selection of an appropriate model is important as it could affect the relevant mechanisms of interaction in such permeation studies.

Changes in lipid diffusion with increased lipid heterogeneity
Lateral diffusion evaluates the mobility of lipids in the xy-plane of the membrane, and is therefore used as a gauge for the Fig. 3 Example snapshot of lipid tails snorkeling in the BH membrane all occurring in one time frame, in the chains sn18:2 (top), sn20:4 (lower left), and sn24:1 (lower right), with the carbon tails in grey. Hydrogen, water, and ion atoms were removed for clarity.   (table below) dynamics of lipid mixing. However, the systems here were only simulated on a sub-microsecond timescale, which means that lipids in the bilayer experience subdiffusion. [60][61][62] Thus, the diffusion coefficients (D) obtained from the slopes of the mean square displacements (MSDs) (Fig. 5) are estimated values of D.
As we are mostly interested in qualitatively comparing the lateral movement of the different lipid types within the various bilayer environments that we have studied, we will use the values of D determined from the slopes of the MSD plots in order to do so. This has been done previously for similar purposes. 62 As cholesterol is well-known to retard diffusion rates, 99 the diffusion constant of POPC in the pure membrane consequently reduces by half, from 9.1 Â 10 À8 cm 2 s À1 (consistent with reported value 100 ) to 4.1 Â 10 À8 cm 2 s À1 in the POPC-CHOL membrane ( Fig. 5a and d). The diffusion of PC lipids (and all other lipids) slow down further in BH800 and BAD800 as the lipid composition becomes more heterogeneous (Fig. 5c and f). We believe that the greater extent of interdigitation in the lipid packing of BH800 and BAD800 compared to POPC-CHOL is able to affect the fluidity of the respective membranes consisting of B50% cholesterol ( Table 2). The only significant difference we observed between the large and small brain membranes was when comparing the lipid diffusion in BH800 and BH ( Fig. 5c and e). In these two healthy systems, the lipids diffuse approximately twice as fast in the smaller BH membrane than in the larger BH800 membrane. Interestingly, no such difference was observed in the two diseased BAD800 and BAD membranes ( Fig. 5f and e). One difference in the lipid compositions between the healthy BH and diseased BAD membranes is that BH contains lipids with generally more polyunsaturated fatty acid (PUFA) chains than BAD. And it has been previously shown that increasing degrees of chain unsaturation increases the rate of diffusion, 101,102 which seems to describe the case between BH and BAD ( Fig. 5b and e). Thus, we believe that the difference between BH800 and BH results from system size effects as well as short simulation time lengths. As the larger membranes of course have more of each lipid species, we will therefore achieve better sampling of the lipid diffusion of a given species in the larger membranes. This would be particularly true for the lipid species which are minor components of the membrane.

Conclusions
In this manuscript, we have conducted atomistic simulations to compare between bilayers containing simplistic (POPC and POPC-CHOL) and complex (BH800/BH and BAD800/BAD) lipid compositions. We present the first atomistic simulations of realistic human brain membrane models (under healthy and Alzheimer's diseased conditions) incorporating 27 lipid types that range in lipid headgroups, lipid tail lengths, degrees of saturation, as well as symmetry and asymmetry in the lipid tails.
Our studies demonstrated effective differences in the properties between the simple and heterogeneous membranes, especially in the dipole potential and lateral diffusion. As homogeneous PC bilayers are commonly used to investigate the permeability of small molecules across the membrane, [103][104][105] it is worth noting that the decreased potential barriers in simple models could alter the permeation mechanisms of small molecules compared to membranes of similar complexity as found in actual biological cells. But of course, other factors at play include the permeate size, membrane organisation phase, etc. 106 Whilst POPC-CHOL seems a good approximation to the complex membranes, there is however a significant qualitative difference in the diffusion dynamics, which could also affect studies of permeability. 107,108 The decreased diffusion of lipid molecules in the model BAD membrane is consistent with the findings that decreased lipid fluidity leads to accelerated amylodogenic processing of the Alzheimer's Precursor Protein. 109 In addition, we have demonstrated that the hydrocarbon tails of the lipid molecules have a tendency to snorkel towards the membrane's interface with the aqueous surroundings. We correlated the likelihood of a lipid tail to snorkel with increasing degrees of unsaturation in the hydrocarbon chains as well as longer tail length. Also, we observe that lipid snorkeling is more common in the more disordered membrane (POPC) where there is less steric hindrance to the motion of the lipid tails. We have also shown that there is a direct correlation between the amount of snorkelling observed in a lipid bilayer and the amount of water that is able to permeate into the membrane. This suggests that lipid snorkeling may play a role in the ability of small molecules to permeate across biological membranes. Continuing to understand the mechanisms that induce changes to the membrane mechanics and biophysics, and therefore cell functions will prove beneficial in developing new therapeutic strategies for the prevention and treatment of Alzheimer's disease and other neurological diseases.

Conflicts of interest
There are no conflicts to declare.