The Signature of Fluctuations of the Hydrogen Bond Network Formed by Water Molecules in the Interfacial Layer of Anionic Lipids

: As the water molecules found at the interface of lipid bilayers exhibit distinct structural and reorientation dynamics compared to water molecules found in bulk, the fluctuations in their hydrogen bond (HB) network are expected to be different from those generated by the bulk water molecules. The research presented here aims to gain an insight into temperature-dependent fluctuations of a HB network of water molecules found in an interfacial layer of multilamellar liposomes (MLVs) composed of anionic 1,2-dimyristoyl-sn -glycero-3-phospho-L-serine (DMPS) lipids. Besides suspending DMPS lipids in phosphate buffer saline (PBS) of different pH values (6.0, 7.4, and 8.0), the changes in HB network fluctuations were altered by the incorporation of a non-polar flavonoid molecule myricetin (MCE) within the hydrocarbon chain region. By performing a multivariate analysis on the water combination band observed in temperature-dependent FTIR spectra, the results of which were further mathematically analyzed, the temperature-dependent fluctuations of interfacial water molecules were captured; the latter were the greatest for DMPS in PBS with a pH value of 7.4 and in general were greater for DMPS multibilayers in the absence of MCE. The presence of MCE made DMPS lipids more separated, allowing deeper penetration of water molecules towards the non-polar region and their restricted motion that resulted in decreased fluctuations. The experimentally observed results were supported by MD simulations of DMPS (+MCE) lipid bilayers.


Introduction
Fundamental biological phenomena such as enzyme activity [1,2], cell adhesion [3], membrane fusion [4], and lipid diffusion [5] are intimately related to the structural features and reorientation dynamics of the water molecules found at the biomembrane interface.In exclusively lipid membranes, often investigated with the aim of simplifying and understanding events on the surface of biomembranes, the arrangement and orientation of water molecules within the interfacial aqueous layer are primarily a function of the magnitude and of the charge distribution on lipid molecule polar head-groups [6,7].In bulk water, a single water molecule is hydrogen-bonded to tetrahedrally surrounded water molecules from the first solvation shell that can be exchanged with those from a second solvation shell owing to the concerted reorientation and cleavage of the existing hydrogen bond (HB) [8].However, a disruption of the tetrahedral ordering of water molecules due to the proximity of the lipid membrane interface [9] imposes significant limitations on water reorientation dynamics.
The investigation of structural and reorientation dynamics of interfacial water usually involves ultrafast vibrational spectroscopic techniques such as vibrational sum frequency Biophysica 2024, 4 93 generation (VSFG) and two-dimensional IR spectroscopy (2D-IR) [10][11][12][13].Despite the inherent limitation of its application to exclusively lipid monolayers, VSFG showed that H-atoms of water molecules near choline (positive) [14] and phosphate (negative) [15,16] moieties are oppositely oriented, and these findings are further corroborated by independent computational studies [17][18][19].Besides identifying differently oriented H-atoms, VSFG showed that oppositely oriented water molecules in the vicinity of choline (positive) and phosphate (negative) groups obey different HB dynamics [20].Specifically, as water molecules near positively charged functional groups display the maintenance of tetrahedral arrangement, their reorientation dynamics resemble those of bulk water [14], whereas near negatively charged functional groups the opposite scenario is found; the loss of the tetrahedral ordering of water molecules [9] is accompanied by the retardation of HB exchange dynamics due to the reduced number of partners in a HB arrangement [9,13,[21][22][23][24][25][26].Unlike VSFG, 2D IR is not an intrinsically surface-sensitive technique, but nevertheless it allows distinguishing the events that occur at different time scales [27]; in particular, by using the appropriate vibrational echo pulse sequence it becomes possible to decipher the position of the oscillator, commonly of an environmentally sensitive C=O moiety, within the heterogeneous sample such as a hydrated bilayer [28].
Along with the mentioned and briefly introduced state-of-the-art spectroscopic techniques, whose main disadvantage is the price, limiting their use to a small number of laboratories in the world, the use of widely available conventional linear FTIR spectroscopy [29] occasionally produces a spark of progress in the characterization of water-based HB systems.Although it contains the averaged signal of all water molecules present in the lipid suspension (bulk and interfacial water), the analysis of the band(s) loaded with water HB network fingerprints [30,31], as well as defining/introducing an arbitrary metric, can enable the expression of the acquired spectroscopic data in terms of water HB network fluctuations [32].
Accordingly, Tong et al. managed to distinguish bulk water molecules from those at the air-water interface by inspecting water librations (δ L ), i.e., intermolecular vibrations between water molecules, the existence of which is the result of HBs [31].Further, Verma et al. conducted a comprehensive concentration and temperature study of aqueous solutions of various salts in order to elucidate the impact of a particular anion/cation on the water combination band (usually centered at about 2130 cm −1 ) that originates from the fundamental transition of water bending (δHOH) and ρ L [30].Because the direction and the magnitude of the water combination band maximum displacement to higher/lower frequencies in the presence of ions coincided with their position in a Hoffmeister series, this approach provided a platform for monitoring both local and global water dynamics through the impact of kosmotropes/chaotropes on a HB network.
The quantification of temperature-dependent fluctuations in a water-molecule-based HB network in the presence of anions of Na + and K + salts (Cl − , Br − ) was presented for the first time in the study conducted by Brkljača et al. [32].Briefly, the water combination band was subjected to a multivariate analysis in which the response of water molecules in a certain temperature range (25-75 • C) was divided into two contributions: the one produced by low-temperature (high-density) water and the other produced by high-temperature (low-density) water [33,34].The observed hysteresis between the responses of these two kinds of water, expressed in terms of mathematical quantities that reflect HB-network-associated phenomena, showed that temperature-dependent fluctuations in a HB network are actually ion-invariant.The validity of the obtained results was confirmed by independent study conducted on equivalent systems investigated using a femtosecond elastic second harmonic scattering (fs-ESHS) technique that produced analogous findings [35].
As the synergy of linear FTIR spectroscopy and multivariate analysis managed to provide an insight into the fluctuations in the water HB network in the presence of anions, this study aimed to apply the same approach in characterizing the fluctuations in the HB network of water molecules found at the interface of negatively charged supramolecular aggregates.Accordingly, multilamellar liposomes (MLVs) of 1,2-dimyristoyl-sn-glycero-3phospho-L-serine (DMPS) were prepared in a phosphate buffer saline (PBS) and examined by FTIR spectroscopy in a temperature range that spanned about ±10 • C around their melting point (T m ).The latter is the temperature at which van der Waals interactions between hydrocarbon chains are reduced, causing lipids to transition from a more ordered gel (L β ) phase to a disordered fluid (L α ) phase, and which for DMPS lipids is reported to be between 37 and 39 • C [36,37].Overall, this L β →L α transition is accompanied by thinning of the lipid bilayer, loss of lipid lateral arrangement, and consequently greater hydration of polar groups and the glycerol backbone that encompasses the HB network between the water molecules in the interfacial layer [38].
By knowing that negatively charged phosphatidylserine (PS) lipids are generally very unstable in terms of curvature variations and local pH change [39], especially in the form of large unilamellar liposomes (LUVs) [40], the presumed fluctuations in the interfacial HB network of DMPS in the L β and in L α phase were enhanced in two ways: (i) MLVs of DMPS were prepared in PBS with three pH values: 6.0, 7.4, and 8.0, at which all DMPS lipids are negatively charged; (ii) the incorporation of flavonoid myricetin (MCE), a non-polar molecule that prefers accommodation in a hydrocarbon chain region [41], which should induce the separation of lipid molecules and provide deeper penetration of interfacial water that consequently modulates fluctuations in their HB network (structural formulas are displayed in Figure 1).The results obtained in this study are aimed to be applied in further studies addressing the fluctuations in the interfacial water layer of other supramolecular aggregates carrying the charges of different magnitudes and distributions.
As the synergy of linear FTIR spectroscopy and multivariate analysis managed to provide an insight into the fluctuations in the water HB network in the presence of anions, this study aimed to apply the same approach in characterizing the fluctuations in the HB network of water molecules found at the interface of negatively charged supramolecular aggregates.Accordingly, multilamellar liposomes (MLVs) of 1,2-dimyristoyl-sn-glycero-3-phospho-L-serine (DMPS) were prepared in a phosphate buffer saline (PBS) and examined by FTIR spectroscopy in a temperature range that spanned about ± 10 °C around their melting point (Tm).The latter is the temperature at which van der Waals interactions between hydrocarbon chains are reduced, causing lipids to transition from a more ordered gel (Lβ) phase to a disordered fluid (Lα) phase, and which for DMPS lipids is reported to be between 37 and 39 °C [36,37].Overall, this Lβ→Lα transition is accompanied by thinning of the lipid bilayer, loss of lipid lateral arrangement, and consequently greater hydration of polar groups and the glycerol backbone that encompasses the HB network between the water molecules in the interfacial layer [38].
By knowing that negatively charged phosphatidylserine (PS) lipids are generally very unstable in terms of curvature variations and local pH change [39], especially in the form of large unilamellar liposomes (LUVs) [40], the presumed fluctuations in the interfacial HB network of DMPS in the Lβ and in Lα phase were enhanced in two ways: (i) MLVs of DMPS were prepared in PBS with three pH values: 6.0, 7.4, and 8.0, at which all DMPS lipids are negatively charged; (ii) the incorporation of flavonoid myricetin (MCE), a nonpolar molecule that prefers accommodation in a hydrocarbon chain region [41], which should induce the separation of lipid molecules and provide deeper penetration of interfacial water that consequently modulates fluctuations in their HB network (structural formulas are displayed in Figure 1).The results obtained in this study are aimed to be applied in further studies addressing the fluctuations in the interfacial water layer of other supramolecular aggregates carrying the charges of different magnitudes and distributions.

FT-IR Spectra Acquisition
FT-IR spectra were measured on an ABB Bomem MB102 spectrometer (ABB Bomem factory, Québec, QC, Canada), equipped with CsI optics and a DTGS detector in a transmission.FT-IR spectra of DMPS (+MCE) lipid suspensions in PBS of three different pH values (6.0, 7.4, and 8.0) were acquired in a sealed cell equipped with CaF 2 windows with the path length of the cell d = 27.37 µm [42].Temperature was regulated by a Specac 3000 Series high-stability temperature controller with heating jacket (Specac, Kent, UK).The spectra of all suspensions were acquired as duplicates in the temperature range 30-50 • C using a heating rate of 1 • C min −1 , at nominal resolution of 2 cm −1 and 10 scans.

FT-IR Spectra Analysis
After measuring FT-IR spectra of DMPS (+MCE) suspensions in PBS with three pH values, the obtained spectra were divided into two separate spectral regions (D 1 , D 2 ): (i) 2980-2820 cm −1 , displaying the bands originated from the (anti)symmetric stretching of methylene (CH 2 ) groups of hydrocarbon chains (ν (a)s CH 2 ), and (ii) 2650-1800 cm −1 where only water absorbs (a combination band that originates from water bending (δHOH) and water libration (ρ L H 2 O); H 2 O comb).The analysis of spectral region D 1 is thoroughly presented in Supplementary Materials, Section S1 (Figure S1), and here the analysis of D 2 spectral region is presented.Two-point baseline correction at the band minima in spectral region D 2 was followed after smoothing (Savitzky-Golay; polynomial of 3rd degree through 50 points [43]).Additionally, the spectral range 2650-2300 cm −1 was eliminated from the analysis because the band originated from ν as CO 2 , the amount of which varies during the measurement, obscures the results [32].The information obtained by analyzing the first spectral region (D 1 ) concerned the determination of the melting temperature (T m ) [44,45]i.e., the temperature at which van der Waals interactions between hydrocarbon chains are weakened-and its change in the presence of MCE (Table S1 in Supplementary Materials), whereas the second one addressed temperature-dependent fluctuations in a HB network formed by water molecules, both bulk ones and those that are found in the interfacial water layers [32].
The prepared spectra underwent further analysis using the Multivariate Curve Resolution-Alternating Least Squares (MCR-ALS) method.This mathematical method examines a data set for main sources of variation and, when applied to spectroscopic data, identifies the spectra of system components with no prior knowledge of their identity, position, or signal shape [46,47].The requirements of MCR-ALS are for data to be structured as a two-way matrix, and for that matrix to be bilinear (i.e., a product of two matrices).In the case of FTIR spectra in this study, raw data form a matrix describing signal intensity at particular temperatures and wavenumbers.The matrix is the result of contributions of smaller matrices describing individual species and their concentrations.This relationship is described mathematically as where D is the raw data matrix, S T is a spectral matrix of individual components, C is the concentration matrix of those components, and E is the matrix of unexplained residuals.Before starting the MCR-ALS analysis, initial guesses of contributions are required.First, the Singular Value Decomposition (SVD) method was applied to estimate the number of components.As Libnau et al. demonstrated [33], the spectral evolution of the water signals as a function of temperature can be presented in the form of a matrix whose rank is two; in other words, two components contain all the information about the system under investigation.Thus, two components were expected to meaningfully contribute to our spectral data [48].Evolving Factor Analysis (EFA) was then utilized to obtain initial guesses of the concentration profiles of those components [46,47].Subsequently, the MCR-ALS equation with the provided initial estimates can be solved iteratively to produce the final S T and C. The MCR-ALS procedure with EFA is publicly available and implemented in MATLAB [49], and further explanations of the theory may be found elsewhere [46,50,51].It is important to note that constraints need to be applied when analyzing FTIR data: all spectral profiles and concentration values must be kept positive (non-negativity constraint), and the individual pure spectra are regarded as unimodal (unimodality constraint).
The spectral and concentration profiles obtained will depend on the examined spectral region and the nature of the chemical species producing signals.Since the examined system contained species that undergo a phase transition in the explored temperature range, a sigmoid profile was expected, with the inflection point around the transition temperature [48].Thus, the contributions of the species existing below T m will be referred to as the low-temperature component (LTC), and those above T m will be referred to as the high-temperature component (HTC).
The concentration profiles generated by MCR-ALS were of sigmoid character for the D 1 spectral range (examined in Supplementary Materials, Figure S2), and of polynomial for the spectral range D 2 .The obtained curves in spectral range D 2 were fitted on a 2nddegree polynomial (parabolic fit; R 2 (DMPS/DMPS + MCE) ≥ 0.998/0.999)expressed in the form y = B 0 + B 1 x + B 2 x 2 .In a mathematical context, the obtained coefficients B 0 (dimensionless quantity), B 1 (in • C −1 ), and B 2 (in • C −2 ) indicate relative concentrations of the particular component (LTC or HTC) at 0 • C, the rate of a decrease in LTC/HTC in favor of an increase in HTC/LTC, and a profile curvature, respectively.As thoroughly described in a previous paper [32], B 2 values obtained for LTC and HTC reflect the hysteresis between the temperature-dependent fluctuations in a HB network meshed by water molecules.

Molecular Dynamics Simulations
Simulations of DMPS bilayers with and without the presence of MCE were conducted using classical molecular dynamics (MD).Homogeneous DMPS bilayers were prepared with the CHARMM-GUI Membrane Builder module [52] and consisted of 200 lipid molecules.DMPS + MCE bilayers were generated by randomly adding 60 MCE molecules to the bilayer (MCE-to-lipid ratio: 0.3).Bilayers were solvated with 50 water molecules per bilayer component.We added 21 Na + and Cl − ions to achieve the experimental ion concentration of 100 mM, along with an additional 200 Na + ions to neutralize the bilayer charge.All simulations were conducted using the GROMACS 2020.0 software [53], CHARMM36m force field [54], and TIP3P water model [55].MCE was parametrized using CGenFF [56].NBFIX corrections for Na-carbonyl interactions have been included in all simulations, as provided in the July 2021 version of the CHARMM36m force field [57].The systems underwent energy minimization followed by 200 ps of equilibration in the NVT ensemble with the V-rescale algorithm.The production run was conducted for 300 ns in the NpT ensemble using the Nosé-Hoover thermostat [58] with a time constant of 1 ps and the Parrinello-Rahman barostat [59] with semi-isotropic pressure coupling and a time constant of 2 ps (target pressure: 1 bar).When conducting the trajectory analysis, the first half of the production run was discarded as equilibration time, and only the last 150 ns were examined.In all simulations, the cutoff for short-range Coulomb interactions and van der Waals interactions was set to 1.2 nm with a switching function turned on after 1 nm.The Particle Mesh Ewald (PME) procedure [60] was used to handle long-range Coulomb interactions.The LINCS algorithm was employed to constrain bonds involving hydrogen.Three-dimensional periodic boundary conditions were used throughout for all systems.And the time step was 2 fs.Visualization and analysis were conducted with Gromacs modules and Visual Molecular Dynamics (VMD) [61].

Results and Discussion
The multivariate analysis of the D 2 spectral region (Figure 2) produced spectral (Figure 2, dotted curves) and, more importantly, concentration profiles of the LTC and HTC (Figure 3).The expected low-frequency shift of the combination band maximum upon heating was found in all systems (~2131 cm −1 at 30 • C and ~2117 cm −1 at 50 • C), regardless of the absence/presence of MCE or the pH (see Table 1).However, it appears that for the DMPS suspension in the absence of MCE there is a slightly greater shift in the band maximum than in its presence, especially for the pH value of 7.4 (2130 cm −1 at 30 • C and 2114 cm −1 at 50 • C), which is the closest to the physiologically relevant pH value.
nm.The Particle Mesh Ewald (PME) procedure [60] was used to handle long-range Coulomb interactions.The LINCS algorithm was employed to constrain bonds involving hydrogen.Three-dimensional periodic boundary conditions were used throughout for all systems.And the time step was 2 fs.Visualization and analysis were conducted with Gromacs modules and Visual Molecular Dynamics (VMD) [61].

Results and Discussion
The multivariate analysis of the D2 spectral region (Figure 2) produced spectral (Figure 2, dotted curves) and, more importantly, concentration profiles of the LTC and HTC (Figure 3).The expected low-frequency shift of the combination band maximum upon heating was found in all systems (~2131 cm −1 at 30 °C and ~2117 cm −1 at 50 °C), regardless of the absence/presence of MCE or the pH (see Table 1).However, it appears that for the DMPS suspension in the absence of MCE there is a slightly greater shift in the band maximum than in its presence, especially for the pH value of 7.4 (2130 cm −1 at 30 °C and 2114 cm −1 at 50 °C), which is the closest to the physiologically relevant pH value.Following the analogy developed in a previous study, these results can be examined in the context of a decrease in temperature-dependent fluctuations in the water HB network when it surrounds MLVs of DMPC+MCE suspensions, which is not the case for pure DMPS lipids.As this approach catches the response of all water molecules, i.e., both bulk (the contribution of which is dominant) and interfacial ones, it seems reasonable to assume that the bulk water contribution is nearly the same for all suspensions, regardless of the pH value or the presence of MCE.The major variable that changes in this consideration is the amount of interfacial HB water and how its HB pattern is changed when MCE is inserted into the DMPS lipids.A potential explanation for this phenomenon is that the motions (primarily librations) of water molecules in the interfacial regions are more restricted when MCE is incorporated in MLVs.If the DMPS lipids are more separated even at lower temperatures (due to the inclusion of MCE), then more water will be able to penetrate into the interfacial region, that is, it will penetrate deeper into the liposome itself and will certainly show different dynamics compared to interfacial water which does not penetrate so deeply into the lipid double layer (as with DMPS itself).
At first glance, polynomial fits of the LTC and HTC (Figure 3) reveal deviations from the experimental values for DMPS suspensions in the absence of MCE.Although the deviations are observed at all pH values, for the DMPS suspension in PBS of pH 6.0, 7. Finally, it should be noted that the incorporation of MCE, in addition to affecting the change in the dynamics of the interfacial layer of water molecules, also affects the separation of or interaction between lipid molecules both in the gel phase and in the fluid phase.The most obvious indicators of the consequences of incorporating MCE into the DMPS lipid bilayer are changes in the lipid T m and the cooperativity of the phase transition [38], both of which are available and can be determined from the temperature-dependent FTIR spectra in the spectral region D 1 (Figures S1 and S2 and Table S1) [48,62].Although the reduction in T m due to the incorporation of MCE was observed only in the DMPS lipid bilayers hydrated with PBS of pH 8.0 (Table S1, Figure S2), which might be associated with greater rigidity and compactness of anionic lipid bilayers in comparison with zwitterionic ones [40,41,44,63,64], in all three lipid suspensions (i.e., at all three pH values) it is evident that the melting is not characterized by such a steep sigmoidal transition in the presence of MCE as it is in its absence, which is a clear indication of a decrease in the cooperativity of the transition [62,65].
The structure and hydration of DMPS bilayers at different temperatures with or without MCE were also examined using MD simulations.Examination of bilayer structural parameters revealed larger spatial separation of DMPS lipids following MCE inclusion, as demonstrated in Table 2.Additional information on the structural analysis can be found in Supplementary Materials.APL and membrane thickness are common parameters describing lipid organization: smaller APLs indicate tighter packing, while larger thicknesses are due to ordered lipid acyl chains with less interdigitation.Thus, the values obtained for pure DMPS bilayers reflect the tendency of lipids to pack more tightly at lower temperatures (gel phase) compared to higher temperatures (fluid phase) [66,67], and are consistent with reports in the literature for similar systems.Wide-angle X-ray scattering data on DMPS at 30 • C pointed to an APL of 0.415 nm 2 [68], which is close to the calculated figure.Also, pressure isotherm studies have shown that DMPS monolayers collapse at an APL of ~0.40 nm 2 [69].Computational studies provided estimations of 47.6 nm 2 in the gel (30 • C) and 55.6 nm 2 in the fluid phase (60 • C) for pure bilayers or ~0.60 nm 2 for bilayers containing small molecules [70].The respective studies also provided bilayer thickness estimates of 4.0-4.4nm, or ~3.5 nm.Experimental data on DMPS bilayer thicknesses are lacking.X-ray diffraction of DMPS MLVs showed the lamellar repeat distance in phosphate buffer to decrease from 6.2 nm to 5.2 nm following an increase in temperature from 20 • C to 53 • C; however, those values include the layer of interlamellar water [71].MCE addition disrupted the ability of DMPS to pack tightly and interfered with the close interactions between lipid chains that are crucial for bilayer organization, resulting in significant APL increase and a decrease in thickness.Greater fluidity of acyl chains is shown by analyzing acyl chain order parameters (Figure S3 in Supplementary Materials), which are significantly decreased in DMPS + MCE systems.Notably, the effect of MCE is larger at 30 • C, where structural parameters are changed enough to be consistent with the fluid phase despite the low temperature, thus confirming membrane disruption.Even visually, the equilibrated DMPS + MCE system at 30 • C is more like the fluid phase of pure DMPS (Figure S4 in Supplementary Materials), and the distribution of chain tilt angles is broadened (Figure S5 in Supplementary Materials).The change in structural parameters followed the same trend in the simulations of MCE + 1,2-dimyristoyl-sn-glycerol-phosphocholine by Sadžak et al. [41].Other bioflavonoids such as quercetin, genistein, or daidzein were also shown to raise APL and lower bilayer thickness [72,73].However, the existence of a fluid phase at 30 • C is in contrast with the experimental results, which have not suggested a meaningful reduction in T m.Still, some changes in the transition process were noted, and it is known from the literature that hydrophobic flavonoids may influence the melting [41,[74][75][76].Considering the relatively small gap of ~8 • C between the simulation temperature and T m , and the known tendency of simulated membranes to melt at temperatures up to 5 • C below experimental values [66], this discrepancy is not unusual.Nevertheless, since the gel phase was not successfully achieved in the DMPS + MCE systems, the comparison of DMPS and DMPS + MCE at 30 • C is considered unreliable.
Since bilayer disruption may lead to changes in water penetration and interfacial interactions, membrane hydration was examined.Density profiles of water (Figure 4) were used to visualize the presence of water molecules in the bilayer.
sual.Nevertheless, since the gel phase was not successfully achieved in the DMPS + MCE systems, the comparison of DMPS and DMPS + MCE at 30 °C is considered unreliable.
Since bilayer disruption may lead to changes in water penetration and interfacial interactions, membrane hydration was examined.Density profiles of water (Figure 4) were used to visualize the presence of water molecules in the bilayer.Water was able to penetrate the surface of the bilayers in all cases, showing some presence among the polar headgroups oriented outwards.In the center of the bilayer, water was excluded due to the hydrophobicity of the acyl chains.Then, the HB network was analyzed in order to confirm lipid-water interactions.As expected, MD simulations showed better hydration and more HBs established in the fluid-phase systems (Table 3), since the separation of lipids due to the expanding surface area allowed for more water to penetrate and bind to the lipids.The impact was notably greater for the DMPS functional groups that are situated deeper in the bilayer, in accordance with findings showing solvent accessibility to be the greatest barrier to lipid group hydration [17].Thus, the COO − group, which sits at the bilayer surface exposed to water, showed no major change in Hbonding, while PO2 − and carbonyl C=O were significantly better hydrated in the fluid phase.Considering the aforementioned unreliability of the DMPS + MCE results at lower temperatures, we cannot draw conclusions on the MCE impact.However, the comparison of higher-temperature simulations with and without MCE suggests an effect of increased water interactions, particularly with the deeply situated C=O group.Water was able to penetrate the surface of the bilayers in all cases, showing some presence among the polar headgroups oriented outwards.In the center of the bilayer, water was excluded due to the hydrophobicity of the acyl chains.Then, the HB network was analyzed in order to confirm lipid-water interactions.As expected, MD simulations showed better hydration and more HBs established in the fluid-phase systems (Table 3), since the separation of lipids due to the expanding surface area allowed for more water to penetrate and bind to the lipids.The impact was notably greater for the DMPS functional groups that are situated deeper in the bilayer, in accordance with findings showing solvent accessibility to be the greatest barrier to lipid group hydration [17].Thus, the COO − group, which sits at the bilayer surface exposed to water, showed no major change in H-bonding, while PO 2 − and carbonyl C=O were significantly better hydrated in the fluid phase.Considering the aforementioned unreliability of the DMPS + MCE results at lower temperatures, we cannot draw conclusions on the MCE impact.However, the comparison of higher-temperature simulations with and without MCE suggests an effect of increased water interactions, particularly with the deeply situated C=O group.
A deeper analysis of the HB networks examined the average number of HBs an interfacial water molecule may form with lipids and surrounding water molecules.The DMPS-water binding patterns (Table S2) show that most water molecules establish 1 HB with 1 DMPS molecule (1-1 pattern), but 15-20% of the water molecules form "bridges", linking 2 lipids together (1-2 pattern).There is also a small amount of water molecules that establish two HBs with a single lipid, while other patterns are negligibly rare.When comparing the systems, it is notable that gel-to-fluid transition increases the percentage of 1-1 patterns from 75.8 ± 1.9% to 80.1 ± 1.9%, which is accompanied by a reduction in 1-2 patterns from 20.1 ± 1.2% to 16.0 ± 1.0%.On the other hand, the water-water binding patterns revealed the amount of additional water-water HBs that each lipid-bound water may establish (Table S3).In all systems, the majority of bonded water molecules at the interface established one or two more HBs with other water molecules, with less than 10% fulfilling their maximum of three additional HBs.Roughly a quarter of water molecules were not bonded to any other water molecules, having been isolated within the bilayer, with all donor or acceptor sites occupied by DMPS or MCE.This is consistent with observations of the loss of the tetrahedral HB network in interfacial water [77].However, the fluctuations in the water network appeared to be insignificant between the systems.There is a small trend of a decreasing number of isolated water molecules (zero extra HBs), accompanied by a small increase in water molecules with three extra HBs when going from gel to fluid.Again, MCE effects at lower temperatures must be dismissed as unreliable.Thus, DMPS-water and to a lesser extent water-water networks show some changes as a consequence of increased water penetration, though the MCE effect is inconclusive.Despite the small scale of the effects on the simulated bilayer segments, in combination with the experimental work they serve to confirm the temperature-dependent change in water behavior at the interface exacerbated by the presence of MCE.Ultimately, as demonstrated by Schönfeldová et al. [36], the melting itself induces redistribution of water molecules in the interfacial water layer.As this feature is changed upon the presence of MCE, as well as the characteristics of the melting itself (see Figure S2), it is somewhat expected that the water HB network features are changed in DMPS ± MCE situations.In the approach applied in this study, this reflects the difference in the fluctuations of a water HB network.Due to the great potential of the former, further development of its methodology is underway.

Conclusions
The research presented in this paper tackles the characterization of the fluctuations of a HB network of water molecules found in the interfacial region of liposomes constituted by negative DMPS lipids.Using the presented multivariate-based approach, two important fluctuation-related findings were discovered: (i) when the pH value of PBS is changed, it appears that for pH 7.4 the fluctuations are slightly larger than for pH values of 6.0 or 8.0; (ii) when MCE, a non-polar flavonoid molecule, is inserted in the hydrocarbon region of DMPS multibilayers, the fluctuations are significantly reduced.The suggested increase in the number of water molecules in the interfacial region upon incorporation of MCE is also suggested by the MD simulations, as well as a perturbation of the membrane ordering.
FTIR spectra are designated with solid curves, whereas spectral profiles of low-and high-temperature components obtained from MCR-ALS with EFA using dotted curves; Figure S2.The concentrational profiles of low-and high-temperature component of DMPS in the absence (left column) and the presence of MCE (right column) suspended in PBS at pH 6.0 (upper row), 7.4 (middle row), and 8.0 (bottom row).Experimentally obtained low-and high-temperature components were labeled with solid curves colored red/wine (DMPS, pH 6.0), orange/dark yellow (DMPS + MCE, pH 6.0), purple/violet (DMPS, pH 7.4), purple/pink (DMPS + MCE, pH 7.4), blue/navy (DMPS, pH 8.0) and cyan/dark cyan (DMPS + MCE, pH 8.0).Their fit on a single Boltzmann profile is designated with dotted curves of the corresponding color;  S1.T m values (T m,1 and T m,2 ) obtained from fitting LTC (T m,1 ) and HTC (T m,2 ) on a single Boltzmann profile and from the intersection point (T m,i ); Table S2.The percentage of H-bonded water molecules that form a particular bonding pattern with DMPS lipids, for studied systems.S3.The percentage of water molecules bonded to DMPS lipids, that also establish additional HBs with other waters.References [36][37][38]48,62,63,65,70,78] are cited in the supplementary materials.

Figure 2 .
Figure 2. Baseline-corrected and Savitzky-Golay-smoothed, temperature-dependent FT-IR spectra of DMPS in the absence (a,c,e) and the presence of MCE (b,d,f) suspended in PBS at pH 6.0 (a,b), 7.4
4, and 8.0 they are the most significant in the temperature range 43-45 • C, 42-48 • C, and 37-44 • C, respectively.Interestingly, this is not observed for DMPS + MCE suspensions, regardless of the pH value of PBS.Another interesting phenomenon observed from the FTIR data of the DMPS suspensions concerns the difference in absolute values of B 1 (∆|B 1 |) and especially of B 2 (∆|B 2 |) coefficients of the LTC and HTC (Table 1); ∆|B 1 | values for the DMPS in PBS with pH values of 6.0, 7.4, and 8.0 are about 0.003 C −1 , 0.009 C −1 , and 0.002 C −1 , whereas for the DMPS + MCE in PBS with pH values of 6.0, 7.4, and 8.0 they are about 0.003 C −1 , 0.001 C −1 , and 0.002 C −1 .Regarding ∆|B 2 | values, for the DMPS in PBS of pH 6.0, 7.4, and 8.0 they are about 0.3-0.5 C −2 , 1 C −2 , and 0.1 C −2 , whereas for the DMPS + MCE in PBS of pH 6.0, 7.4, and 8.0 they are about 0.3 C −2 , 0.1 C −2 , and 0.3 C −2 .When the uncertainty is taken into consideration, it can be seen that ∆|B 1 | and ∆|B 2 | for all suspensions are the same value.However, it cannot be ignored that dissipation in ∆|B 1 | and especially in ∆|B 2 | is much larger for DMPS only, which contrasts the situation when MCE is incorporated in the bilayers.

Figure 4 .
Figure 4. Number density profiles of DMPS, MCE, and water molecules along the z axis of the simulation box, obtained from molecular dynamics simulations: (a) DMPS at 30 °C, (b) DMPS at 50 °C, (c) DMPS + MCE at 30 °C, and (d) DMPS + MCE at 50 °C.

Figure 4 .
Figure 4. Number density profiles of DMPS, MCE, and water molecules along the z axis of the simulation box, obtained from molecular dynamics simulations: (a) DMPS at 30 • C, (b) DMPS at 50 • C, (c) DMPS + MCE at 30 • C, and (d) DMPS + MCE at 50 • C.
Figure S3.Order parameters of acyl chains for DMPS and DMPS + MCE systems at 30 • C and 50 • C, obtained from molecular dynamics simulations; Figure S4.Representative snapshots of lipid bilayers in the production runs of molecular dynamics simulations: (a) DMPS at 30 • C, (b) DMPS at 50 • C, (c) DMPS + MCE at 30 • C and (d) DMPS + MCE at 50 • C. Lipids are shown in Line mode with the acyl chains painted in gray.Phosphorus atoms are shown as gold spheres.MCE is shown in Licorice mode in violet, and water molecules as cyan dots.Black outlines represent the edges of the simulation boxes; Figure S5.The histograms showing the distribution of chain tilt angles of DMPS and DMPS + MCE systems for sn-1 chain (a) and sn-2 chain (b).; Table Pattern 1-1 indicates one water molecule bound to one lipid.Pattern 1-2 indicates one water molecule bound to two different lipids.Pattern 1-3 indicates one water molecule bound to three lipids.Pattern 1=1 indicates one water molecule forming two bonds with the same lipid molecule.Pattern 1-12 indicates one water molecule bound to one lipid through one bond, and to one another lipid through two bonds; Table

Table 1 .
Parameters obtained from polynomial regression (B 0 , B 1 , B 2 ) of the results produced by multivariate analysis conducted on spectral range 2300-1800 cm −1 .

Table 2 .
Structural parameters of DMPS and DMPS + MCE bilayers obtained from molecular dynamics simulations.
a In nm 2 .b In nm.

Table 3 .
Average number of hydrogen bonds (HBs) between select molecules or functional groups for DMPS and DMPS + MCE systems at 30 • C and 50 • C.