Mechanical Properties Determination of DMPC, DPPC, DSPC, and HSPC Solid-Ordered Bilayers

Lipid bilayers are active participants in many crucial biological processes. They can be observed in different phases, liquid and solid, respectively. The liquid phase is predominant in biological systems. The solid phase, both crystalline and gel phases, is under investigation due to its resilience to mechanical stress and tight packing of lipids. The mechanical properties of lipids affect their dynamics, therefore influencing the transformation of cell plasma and the endomembrane. Mechanical properties of lipid bilayers are also an important parameter in the design and production of supramolecular lipid-based drug delivery systems. To this end, in this work, we focused on investigating the effect of solid phases of lipid bilayers on their structural parameters and mechanical properties using theoretical molecular dynamics studies on atomistic models of whole vesicles. Those include area per lipid, membrane thickness, density vesicle profiles, bending rigidity coefficient, and area compressibility. Additionally, the bending rigidity coefficient was measured using the flicker noise spectroscopy. The two approaches produced very similar and consistent results. We showed that, contrary to our expectations, bending rigidity coefficients of solid-ordered bilayers for vesicles decreased with an increase in lipid transition temperature. This tendency was reverse in planar systems. Additionally, we have observed an increase of membrane thickness and area compressibility and a decrease of area per lipid. We hope these results will provide valuable mechanical insight for the behavior in solid phases and differences between spherical and planar confirmations.


Calculation of order parameter for vesicle systems
Carbon-hydrogen order parameter of lipid tails is often used to assess force field accuracy. In this section we are presenting the change of order parameter throughout whole simulation for both inner and outer membranes. Three carbons were selected on both tails -second (C22, C32), eighth (C28,C38) and fifteenth (C215,C315). For DSPC vesicle system thirteen carbon (C213,C313) was selected instead of fifteenth for obvious reasons. Order parameter was calculated by established protocol for united-atom 4 with slight adaptation for vesicle systems. Namely, the vector between vesicle center and position of phosphorus atom for given lipid is treated as membrane normal. Carbon atoms for POPC system located closer to lipid head (namely C22,C32,C28 and C38) became stable after 20 ns of simulations. For C215 and C315 order parameter semi-stabilized after 20 ns, after which constant small drift could be observed. However such drift is to be expected in carbon atoms at the end of carbon tails. Parameter evolution is presented in Figure S4. Carbon atoms in DMPC vesicle system located closer to lipid head (namely C22,C32,C28 and C38) became stable after 10 ns of simulations for sn-2 tail and after 8 ns for sn-1 tail. For carbon atoms at the end of tails the order parameter semi-stabilized almost instantly after less than 2 ns, after which remain constant when a small drift could be observed around 28 ns. Additionally, significant difference was observed in parameter value between the inner and outer leaflets. This was not observed only for middle carbon atom. However it can be concluded that system is thermally stable. Order parameter evolution in time is presented in Figure S5. For DSPC system atoms closest to head, namely C22 and C23, obtained stability after 10 ns for both inner and outer leaflets. However the fluctuation of the parameter in time was much higher than observed for POPC or DMPC systems. Similar tendency was observed for C28 and C38 atoms. Additionally it takes longer for inner leaflet to obtain stability. A constant drift was observed for last carbon atoms in tails, the drift was higher in inner leaflet. Except for C22 atom, order parameter values were different between the inner and outer leaflet. Despite that the parameters and stable, therefore system can be treated as thermally stable as well. Order parameter evolution in time is presented in Figure S6. For DPPC system and carbon atoms closest to head, stability was obtained after 15 ns for inner leaflet and after 5ns for outer leaflet. For C28 and C38 stability was obtained relatively quick after 5 ns. However, either small drift or high fluctuation can be observed for carbon atom in inner leaflet. A constant drift was observed for last carbon atoms in tails, after obtaining semi-stability at around 3 ns. It can be concluded that system can be treated as thermally stable after 15 ns. Order parameter evolution in time is presented in Figure S7. For HSPC system two different lipid types can be found. It can be observed that order parameter value is the same for both DSPC and DPPC lipid molecules. Only exception from this tendency can be observed in C213 and C313 atoms in outer leaflet. For carbon atoms closes to head as well as for middle carbon atoms stability was observed after 8ns. For carbon atoms at the end of tail semi-stability was observed after 4 ns for outer and 10 ns for inner leaflet. However constant small drift was present in both cases. Nevertheless, it can be safely assumed that systems are thermally stable after 10 ns of simulations. Order parameter evolution in time is presented in Figure S8.

Flicker-noise analysis approach for Molecular Dynamics simulation
For bending rigidity determination in case of molecular dynamics study the fluctuation analysis is done on whole liposome. However, in case of flicker-noise spectroscopy, only cross-section is analysed and used to mathematically re-establish fluctuations on whole vesicle. To compare accuracy of such approach, fluctuation contour of cross-section of vesicle in molecular dynamics studies was extracted and analysed similarly to flicker-noise images. Fluctuation spectra was collected for range of 3 nm and under five different angles as visualized in Figure S9. Bending rigidity calculated using Braun and Sachs 5 approach was equal to κ=17,86·k B T. When analysed the slice under 0 degree angle it was equal to 17.9±0.6·k B T and 17±3·k B T from average-based and statistical approaches respectively. In other angles result were within margin of error: for 30 degree slice it was equal to 17.2±0.7·k B T and 16±4·k B T, for 45 degree slice -16.6±0.6·k B T and 15.5±4.3·k B T and for 90 degree slice -17.7±0.7·k B T and 17±4·k B T. In each case first result from average-based approach is presented followed by statistical approach. Only in case of 60 degree angle slice result was slightly different -13.9±0.5·k B T and 14.3±3.4·k B T. Despite this single discrepancy, it can be concluded that determination of fluctuation spectra from cross section is accurate.

Planar bilayer MD simulations details and mechanical parameters determination
Mechanical parameters for investigated lipids were also determined in planar lipid bilayer configuration. All full-atomic simulation simulations were performed with the NAMD 1 software and united-atom CHARMM36 force field 2 under NPT conditions. Specifically, planar lipid bilayers were generated using CHARMM-GUI membrane builder, which was followed with hydrogen removal to reflect united-atom force field. Each investigated bilayer consisted of 648 lipids (324 for each leaflet). For HSPC system the bilayer consisted of 574 DSPC molecules and 74 DPPC molecules. Other options were the same as in vesicle system simulations. Planar bilayer simulations were run for at least 100 ns. Last 50 ns were used for mechanical parameter determination.
To determine mechanical parameters (focusing on bending rigidity coefficient, but also tilt modulus) a real-space fluctuation (RSF) method was used 6 . Specifically, a probability distribution for both tilt and splay is determined for all lipids over all analyzed time steps. Tilt θ is defined as an angle between the lipid director (vector between lipid head -midpoint between C2 and P atoms -and lipid tail -midpoint between last carbon atoms) and bilayer normal. Lipid splay S r is defined as divergence of an angle formed by the directors of neighboring lipids providing that they are weakly correlated. The method, along with equations and calculations, is thoroughly described in given references. APL in planar bilayer simulation is simply determined by dividing box area over number of lipids, which is followed by averaging over analyzed time steps. In table 1 obtained parameters are presented. In figures S10-S14 obtained probabilities P(θ) and P(S r ) along with model fit to the potential of mean force (PMF) for each bilayers are presented. Area compressibility is determined using same algorithm as for vesicle systems 7 .

Flicker-noise detailed results
Presented in main paper bending rigidity coefficients determined in flicker noise spectroscopy were averaged values. They are averaged over at least 10 individual vesicles. However average value do not fully show the diversity of individual measurements. To this end the values of bending rigidity coefficient for individual vesicles as well as their average values are presented in this paragraph ( Figures S15-18). They are presented for both average-based approach and statistical one. Furthermore image of vesicles are shown to further emphasize the difference in their shape, which is stated in main paper.

Time-stability of solid-ordered vesicles shapes
Vesicles created from lipids with T m higher than room temperature exhibited oddly-rectangular shape rather than typical quasi-spherical one. This was more visible the higher was the T m of lipid -namely observed 'bilayer wrinkles' were common view in DSPC, DPPC and HSPC, while they were less visible in DMPC. In this section time stability of this peculiar bilayer shape are shown. As one can see in Figures S9 and S10, visible 'wrinkles' can be seen stable for over a minute. Furthermore, the changes in the vesicle shape is mostly due to rotation of vesicle rather than change of individual wrinkles.