A mechanistic perspective on the effect of cholesterol in phospholipid bilayers

Cholesterol (CHOL) is one of the most important components of plasma membranes of higher cells and one of the main factors for the formation of (nano)domains. In this work, molecular dynamics simulations of mixtures of CHOL with DPPC (saturated lipid) and DLiPC (unsaturated lipid) as standard phospholipids (PLs) are presented in a wide range of CHOL concentrations. The key idea is to systematically extract all structural and enthalpic properties relevant to the formulation of a lattice model of these systems and express them in dependence of the acyl chain order parameters. Detailed interpretation is simplified by the observation that, to a good approximation, the interaction effects do not depend on the total CHOL concentration, but only on the local CHOL arrangement. The resulting information can be used to motivate the agglomeration of CHOL molecules, the relevance of entropic rather than enthalpic effects for understanding the stronger influence of CHOL on DPPC compared to DLiPC, or the thermodynamic background of raft formation. It is verified that the interaction functions hardly change during the transition from binary to ternary mixtures, suggesting the applicability of the concepts to more complex mixtures.


Introduction
The prominent role of CHOL in the formation of functional nanodomains has been widely accepted. [1][2][3][4][5] Since its initial formulation as lipid rafts, 6 many efforts have been made to gain a mechanistic understanding of their properties in live cells and their associated biological function. The nanodomains are associated with several vital membrane functions, like trafficking and signaling, activation, and quenching of proteins (through enrichment), and are proposed to be an attack site for pathogens. [7][8][9] Surprisingly, nanodomains could, to date, not be directly observed in cell membranes due to resolution limits in time and length. 4,[10][11][12][13] Computational models are, therefore, of great value for the understanding of the underlying mechanisms of nanodomain formation. The problem of molecular dynamics simulations is the limitation in the compositional complexity of live cells, even though some breakthroughs have been made. 14,15 Additionally, the nonequilibrium state that characterizes life poses a problem for the formulation of model systems. 16 Simplified (lattice) models reduce a complex relation into simpler sub-issues, each reproducing a specific property of the complex relationship. Thus, each model adds a facet to the big picture. Concepts to be elucidated included the determination of coexisting domains as a function of composition and temperature, domain sizes and their fluctuations, the understanding of the properties and the resulting impact of CHOL in membranes, or the lateral isotropy of lipid interactions (e.g., for hybrid lipids, or CHOL's smooth and rough face).
Lattice models can be characterized by their design choices for the lattice geometry, the interaction parameters, and the representation of each site. Neutron scattering experiments 17,18 and atomistic MD simulations [19][20][21][22] have shown that acyl chains form hexagonal patterns which might guide the choice of the lattice geometry. The hexagonal/triangular geometry, indeed, was utilized in many models that were formulated to reproduce properties of CHOL-like phase behavior, [23][24][25][26] nano domain size, [27][28][29] and segregation behavior. [30][31][32][33] Most of these lattice models use pairwise interactions. The parameters for the interaction functions are derived from experiments or tuned to give correct system behavior. To the best of our knowledge, all these models require phenomenological input to properly fit CHOL characteristics and typically use a constant PL-CHOL interaction.
In previous work, our group has developed a lattice model for mixtures of the phospholipids (PL) DPPC (saturated PL) and DLiPC (unsaturated PL). 34,35 One key aspect was the consideration of the PL-chain order parameter S as one degree of freedom. The model development occurred in two steps. At level I, the neighbor functions, expressing the number of neighbors as a function of S and the interaction functions (equivalently: enthalpy function), reflecting the pairwise lipid-lipid interactions, again in dependence of the chain order parameters, were extracted from molecular dynamics (MD) simulations of the DPPC systems, the DLiPC systems, and the DPPC-DLiPC mixtures. Note that these energies incorporate all structural effects that describe lipid ordering (e.g. spontaneous curvature, which was reported to be a significant factor for lateral structuring asymmetric bilayers 36 ), and represent the mean properties of such systems. This information alone provided a new perspective on the properties of these systems and a qualitative understanding of the different aspects of the system. At level II, the S-dependent chain entropy was adjusted in an iterative procedure in such a way that the Monte Carlo simulations of the lattice model finally yielded the same distribution of order parameters as in the MD simulations. With this complete model different aspects could be positively compared with the results of MD simulations, such as the reproduction of the phase transition temperature or the time dependence of the phase separation dynamics in case of DPPC-DLiPC mixtures. Furthermore, a detailed comparison between atomistic and coarse-grained force fields was possible, revealing the fine-tuned interplay of enthalpy and entropy. 35 In recent work it has been shown that the properties of PL-CHOL mixtures are mainly driven by the the nearest neighbor PL-PL and PL-CHOL interactions as well as entropic effects. 37 In contrast, interleaflet and water interactions hardly varied with the CHOL concentration, indicating that interleaflet do not play a role in the ordering effect of CHOL.
Based on this insight, in this work we formulate a lattice model for PL-CHOL mixtures which, on the one hand, conceptually is similar to our previous work on pure PL-systems but, on the other hand, due to the major structural differences of CHOL and PLs, involve major modifications as compared to the pure PL-based system. CHOL concentrations up to 30% are analyzed. Here we restrict ourselves to the level I which already allows a closer understanding of the specific impact of CHOL on DPPC or DLiPC. It provides answers to, e.g., the following questions: To which degree can the well-known condensing effect of CHOL be taken into account? Is the PL-PL interaction significantly modified by the presence of CHOL? Is it possible to understand the ordering effect of PL due to the presence of CHOL?
Is it possible to characterize the system by neighbor and interaction functions which can be used for all CHOL concentrations? The last question is of particular relevance because a positive answer would strongly increase the applicability of a lattice model, e.g., to study domain formation on the mesoscale for which locally averaged CHOL concentrations may vary quite significantly within the overall system.

Computational details MD simulations
All bilayers containing random mixtures of DPPC, DLiPC, and CHOL were prepared using the CHARMM-GUI web-based graphical interface. 38 The simulations of the DPPC bilayers, the DLiPC bilayers, and the PL/CHOL bilayers were run for 1 µs. The initial 50 ns of each trajectory were discarded to assure the exclusive analysis of equilibrium data. The simulations were run at temperatures from 290 K to 350 K in 10 K steps. Bilayer compositions and trajectory lengths are listed in table S1 in the SI.
All simulations were conducted using the CHARMM36 force-field 39,40 and Gromacs 2020 software package. 41,42 The systems were equilibrated using the 6-step established CHARMM-GUI parameter set, gradually decreasing restraints on lipid head group positions and chain dihedral angles (see table S2). The TIP3P model was used as water model. 43 The Nosé-Hoover algorithm was used to maintain the temperature with a coupling constant of 1 ps, coupling bilayer and solvent separately. 44 To maintain the pressure at 1 bar the Parrinello-Rahman barostat was used with a semiisotropric coupling scheme, a coupling constant of 5 ps and compressibility of 4.5x10 −5 bar −1 . 45 Particle mesh Ewald electrostatics were used with a real-space cutoff of 1.2 nm. 46 The Lennard-Jones potential was shifted to zero between 1.0 and 1.2 nm, with a cutoff of 1.2 nm and the nonbonded interaction neighbor list was updated every 20 steps with a cutoff of 1.2 nm.

Parameter definitions and software
Handling and analysis of MD data heavily relied on the MDAnalysis package for Python. 47,48 All molecular visualizations in this work were created using VMD. 49 The radial distribution functions (RDFs) were calculated using the Gromacs analysis tools. 42 All RDFs were calculated individually for each leaflet using in-plane distances and if not stated otherwise were averaged. For the PL molecules, the phosphor atoms and, for CHOL, the hydroxyl oxygen atoms were used as reference positions.
The lipid chain order parameter S is defined as S = 1.5 cos 2 θ − 0.5 , where θ denotes the angle between the vector spanned by every second carbon atom in a lipid chain and a the bilayer normal (z-axis). This definition of calculating the order parameter differs from the commonly used S CH order parameter, as we do not consider hydrogen atoms in this calculation. Thus S = −0.5 is the lowest order state, while S = 1.0 is the highest order state (linear chains). This has the advantage of easier comparability between the parameters of saturated and unsaturated PLs.

Results and discussion
Lateral bilayer geometry of binary PL/CHOL mixtures  The structuring is stronger with increasing CHOL concentration, especially for DPPC. In contrast, the CHOL concentration exhibits no influence on the head group structure (figure 1 top) and accordingly we find that the chain level is more sensitive to changes of the bilayer composition.
The DPPC-CHOL RDFs (figure 1 bottom) show very similar behavior compared to the head and tail levels of the PL-PL. Generally, the structuring between DPPC and CHOL is weaker than between DPPC and DPPC. The chain level RDFs also show increased structuring by CHOL, while CHOL does not affect the head level structure. Interestingly, the maxima of the DPPC-CHOL chain level RDFs roughly coincide with the DPPC chain level  The PL neighborhood around CHOL resembles the neighborhood of DPPC and DLiPC, except that CHOL has, on average, one additional PL neighbor. This increase may be related to the fact that there exists a finite probability that the projection of the CHOL hydroxy oxygen atom on the membrane plane may be close to the projection of the PL head-group.
However, all three conclusions, drawn for the PL-PL neighborhood, hold here as well.
Using the presented results, we are in the position to find a lattice geometry, based on the lipid head group region, that can incorporate the average lipid neighborhood. The PL head group region is not affected by the presence of CHOL. Therefore, PLs and CHOL should reside on separate square lattices. In that framework, the CHOL lattice is a sub-lattice with void sites, enabling the simulation of different CHOL concentrations. A minimum version of the resulting lattice model could be formulated in the way depicted in figure   3. The slight deviation from an optimal number of 4 neighbors on a square lattice can be handled by including an appropriate scaling factor N/4 (N : total number of neighbors for the specific CHOL environment) for the individual Pl-PL interaction functions, thereby correctly accounting for the total interaction. This minor correction was already successfully used in previous work 51 and brings in a very minor mean-field character into the theoretical description.
From the RDF plots of the chain level it is obvious, that the number of PL (chain) neighbors around PL increases with increasing CHOL concentration. Thus, it would no longer be possible to define neighbor functions of chains without taking into account the overall CHOL concentration.

PL-PL related interaction functions
The interaction function between lipids, which form the key aspect of the lattice model, is influenced by many factors (e.g., variation of lipid head groups, chain length, and saturation).
These factors all influence the inherent order of a lipid. In previous work, 37 we evaluated the interaction of a PL molecule with its environment and expressed the total interaction as a function of the PL's order parameter S. Indeed, a significant dependence on S was observed.
Here we go one step further and characterize the interaction between the individual adjacent lipid pairs with indices i and j as a function of the pair's chain order parameters S i and S j . In particular, we are interested in the effect of CHOL on this interaction which should reflect CHOL's unique ability to alter the phase behavior in lipid mixtures. The interaction functions are determined as averages over the temperature ranges from 330 to 350 K (DPPC) and 290 K to 350 K (DLiPC), respectively.
The interaction energy matrix E(S i , S j ) of a pair ij, shown in figure 4, is necessarily symmetric. The interaction increases with increasing order parameters, and the CHOL concentration does not change the general features of the energy matrix. Following our previous work, we aim to express the interaction matrix as a 1D function where the argument depends on both order parameters. This approach is convenient from a statistical perspective (less noisy) and allows a more direct view on the dependence of the interaction energy on the order parameter. Of course, when actually approaching level II (employing neighbor and enthalpy functions to get the entropy) of the lattice modeling, one might still have the option to work with the full 2D interaction matrix. We now have multiple options to simplify the interaction description (i.e., reducing the two-dimensional to a one-dimensional function). In previous work 51 we chose the pair's average order parameter S avg for PL mixtures. For this purpose we calculated the weighted average for given average order parameter. Here we propose another scheme, using the minimum order parameter S min averaging over horizontal/vertical entries in the interaction matrix and show that this choice is superior to our previous one. An optimal choice minimizes the error of using the 1D variant as compared to the complete 2D interaction matrix. Each CHOL concentration is analyzed separately. Specifically, we calculated for each value of the respective order parameter (S avg or S min ) the weighted variance, see figure 5, where its square root, i.e. the standard deviation, is displayed. It shows that in the relevant regime of positive order parameters, S min is a better choice than S avg . Furthermore, when weighting each value with the occurrence probabilities of the individual order parameters, one obtains an overall quality (figure S1) which also shows the superiority of the S min -scheme, in particular for DPPC.
In figure 6 top we show the interaction enthalpy between the PL-PL pairs as a function of S min for different CHOL concentrations. In agreement with previous work on DPPC bilayers without CHOL, 51 the interaction of DPPC molecules is stronger for chains with high order (figure 6 top). The interaction energy reaches a plateau at around 60 kJ/mol for S min smaller than a threshold value that depends on the CHOL concentration. The DLiPC-DLiPC pair interactions behave quite similarly. Importantly, the PL-PL interaction becomes weaker for higher CHOL concentration. The interpretation is straightforward: since one PL also wants to optimize its interaction with adjacent CHOL molecules (particularly relevant for high CHOL overall concentration), the efficiency of its interaction with the adjacent PL is reduced.
To get more specific information about the effect of CHOL on the PL-PL interaction, we additionally record the number of direct CHOL neighbors N c and display the PL-PL This energy difference for DPPC and DLiPC pair interactions is depicted in figure S2 in the SI. The differences between the energies at 330 K and 350 K lie below 2 kJ/mol. Only the interaction energy for the unlikely case of 30% CHOL and 0 CHOL neighbors shows a higher temperature dependence. Consequently, we infer that the temperature dependence of the interaction as a function of S min and N c is basically temperature independent in the investigated temperature range. right.
For small order parameters the absolute interaction values differ by about 2 kJ/mol between profiles from ternary and binary bilayers whereas they are basically identical in the limit of large order parameters. The interaction in the binary mixtures is stronger for DPPC (decreased packing in ternary mixtures) and weaker for DLiPC (increased packing in ternary mixtures). Thus, the energy gain for ordering DPPC acyl chains is slightly weakened in the ternary mixtures. However, due to the smallness of these effects we may state the PL-PL interaction is only slightly affected when going from the binary to the ternary bilayer. Thus the present results suggest that the interaction functions, derived from pure PL mixtures, can be, to a good approximation, used also in more complex mixtures.

Cholesterol-related interaction functions
Next, we investigate the characteristics of the PL-CHOL and CHOL-CHOL interactions.
The interaction profile of a PL with CHOL as a function of the PL order S P L is shown in figure 8 top. Equivalent to the discussion of the PL-PL interaction functions, we also calculated the interaction for a given number N c of CHOL neighbors of the respective PL- For DPPC and for all CHOL concentrations and N c -values the strength of the DPPC-CHOL interaction energy displays a maximum for S P L ≈ 0.65. Thus, the optimal PL chain order for interaction with CHOL is not an all-trans but a slightly disordered chain. The DLiPC-CHOL interaction displays a monotonous dependence on S P L . However, this may be related to the fact that for the unsaturated DLiPC-lipid no order parameters beyond 0.7 are accessible. Interestingly, DLiPC shows a stronger dependence on the order parameter than DPPC which results in a stronger interaction by approx. 4 kJ/mol at S P L =0.6. The differences between profiles of different CHOL concentrations lie between 0 and -5 kJ/mol and, importantly, hardly depend on S P L (figure 8 bottom). In analogy to the PL-PL interaction, the variation with order parameter is less than 2kJ/mol. The CHOL-CHOL interaction does not directly affect the PL order parameter, but it determines how favorable an agglomeration of CHOL is from an energetic point of view. We show the CHOL pair interaction as a function of the pair's CHOL neighbors N c in DPPC and DLiPC bilayers and the overall CHOL concentration in figure 9. Qualitatively, we see the same tendency as for the other two interactions: stronger interaction for less CHOL neighbor and, as a much smaller effect, higher CHOL content. In DPPC mixtures, the concentration dependence is slightly higher with approximately 2 kJ/mol per 10% CHOL.
In DLiPC mixtures, the concentration dependence is roughly half as strong as in DPPC. The interaction decreases with 1 to 1.5 kJ/mol per N c , both in the DPPC and DLiPC mixtures, and the overall interaction strength between CHOL is slightly higher in DLiPC mixtures.

Comparison with MD results
Single PL-CHOL mixture: from the lattice model to MD results Here we discuss for interesting properties of CHOL-containing membranes, as obtained from MD simulations, how the results can be interpreted from the properties of the neighbor and interaction functions.
We start with the evaluation of the order parameters in the different CHOL-containing bilayers, shown in figure 10. The average order of the DLiPC-containing bilayers increases linearly with CHOL concentration and decreasing temperature. The CHOL ordering is roughly twice as strong for DPPC as it is for DLiPC. For DPPC, we see a phase transition from gel to the liquid phase at around 320 K. The properties of the PL-CHOL interaction function suggest a possible explanation. In the relevant range of order parameters below S P L ≈ 0.7 the PL-CHOL interaction strength increases with increasing order parameter S P L . Thus, addition of CHOL and the concomitant increasing relevance of the PL-CHOL interaction for the overall energy balance prefers higher PL order parameters. Since the interaction functions are only determined for temperatures in the liquid phase, we restrict our discussion to this regime.
However, a closer analysis of the MD results reveals several important features, for which a fresh perspective is provided based on the results of this work. (1) Why is the order parameter for the CHOL-free case higher for DPPC (≈ 0.35) than for DLiPC (≈ 0.2)? The PL-PL interaction functions suggest that the enthalpic gain upon ordering is not stronger than for DLiPC (compare, e.g., the increase of energy between order parameters of 0 and 0.4 in figure 6 which is approx. a factor of 1.5 larger for DLiPC). Although this effect is slightly counterbalanced by the fact that the number of nearest PL neighbors is smaller by a factor of ≈ 0.75 for DLiPC bilayers compared to the DPPC bilayers, these considerations do not explain the significant difference between the order parameters of DPPC vs. DLiPC.
Thus, this observation suggests that the difference in chain entropy must be used to explain this difference. Indeed, due to the disordered nature of DLiPC one may envisage that for DLiPC there is a significant entropy penalty to reach more ordered states. Among others, this also reflected by the observation that for DLiPC order parameters S min above 0.5 are hardly visited whereas for DPPC order parameters up to 0.8 can be observed. (2) Why is the impact of CHOL on the chain order much larger for DPPC than for DLiPC? Again, this cannot be explained by the specific form of the PL-CHOL interaction function since, if at all, it would predict a stronger impact of CHOL on the DLiPC order. However, the chain entropy argument, already used in (1), would again explain the additional suppression of states with higher order parameters for DLiPC. In contrast, the enthalpic gain at higher order parameters upon addition of CHOL is less suppressed by the chain entropy, which naturally favors smaller order parameters. (3) Why is the impact of CHOL on the DPPC order particularly strong at T = 320 K? We envision that this effect is related to the critical effects close to the transition to the gel state. As shown in figure 11 the fluctuations of the order parameter for pure DPPC is particularly large at T = 320 K. In analogy to the general system behavior at critical points (e.g. for the Ising model close to T c ), the system is particularly susceptible to external perturbations such as the addition of CHOL.
Finally, one may wonder whether this insight may help to rationalize the emergence of rafts for the ternary mixtures, i.e. the presence of a DLiPC-domain with low order parameters and a DPPC-CHOL-domain with higher order parameters. We have seen that the impact of CHOL on DLiPC domains is very small (see figure 10). As a consequence there is hardly any gain in free energy due to the mutual interaction. In contrast, for DPPC-CHOL domains the DPPC order parameter is strongly increased which implies a significant free energy gain as compared to the pure DPPC domain. Thus, these observations suggests that the free energy is lower if CHOL forms the liquid ordered phase together with DPPC, which is exactly what happens upon raft formation.
As another observable we characterize the agglomeration of CHOL. It is reflected by the average number of CHOL neighbors around CHOL which we extract from the CHOL-CHOL RDF for distances smaller than 1 nm. When starting from the RDF, this value is automatically normalized by the number of nearest CHOL neighbors in case of a statistical distribution of CHOL molecules. The results are shown in figure 12. For all temperatures Naturally, the significant interaction between adjacent CHOL molecules is likely to drive this agglomeration effect which in the considered range of temperatures is basically Tindependent. Since this interaction is stronger for the DLiPC-CHOL systems, it naturally emerges that the agglomeration effect is somewhat larger for that case. Closer inspection of the PL-CHOL interaction shows that for fixed overall concentration it (approximately) shows a linear dependence on the number N c of CHOL neighbors. Thus, the overall contribution of the PL-CHOL interaction to the total interaction energy does not depend on the arrangement of CHOL molecules.

Conclusion
In this work, we quantified the role of cholesterol in bilayer mixtures with DPPC and DLiPC by studying the neighbor and interaction functions as extracted from MD simulations. It was essential to express these functions in terms of the acyl chain order parameter. They form the core ingredients of a lattice model of PL-CHOL mixtures, reflecting the level I contributions and, thus, in particular the enthalpic part. Furthermore, from the properties of the neighbor functions we came up with a straightforward suggestion for the topology of the lattice model. In this way a minimalistic visualization of the structure of PL-CHOL mixtures is possible.
It could be shown that the pair interaction of two phospholipids, which in general depends on both order parameters independently, can be expressed as a function of an effective order parameter which was chosen as the minimum S min of both order parameters. Among others the properties of the interaction functions clearly rationalize clustering of CHOL for small CHOL concentrations, as seen in the MD simulations.
Surprisingly, the PL-PL and PL-CHOL interaction functions would suggest that, first, the tendency of the lipid chains to attain higher order parameters is stronger for DLiPC and, second, the impact of CHOL should be stronger for DLiPC as well. This contradicts the actual behavior, seen from the MD simulations. The present analysis now clearly enables the relatation of the observations to the entropic contribution of the resulting free energy, reflecting the tendency of the acyl chain to acquire disordered configurations.
A priori it was not evident that the information about the neighbor and interaction functions, extracted from the MD simulations, can be formulated, within a very good approximation (typically less than 2 kJ/mol), in such general terms. This analysis may naturally serve as the basis for the actual analysis of the lattice model. However, this requires the additional determination of the acyl chain entropy via an iterative procedure in analogy to previous work on cholesterol-free mixtures which requires techniques such as used in previous work on CHOL-free lipid membranes and which, thus, brings in additional degrees of complexity. In this way, a more detailed thermodynamic perspective on the nature of raft formation in these model membranes would be conceivable, since the individual enthalpic and entropic contributions would be explicitly known.

Supporting Information Available
The following files are available free of charge.   1000  1000  2  25  NVT  1000  400  3  25  NVT  400  200  4  100  NpT  200  200  5  100  NpT  40  100  6 100 NpT 0 0 Figure S1: Sum of variances per order parameter bin (see figure 5) weighted by the respective order parameter distributions for the S min and S avg scheme.