Influence of cholesterol on hydrogen-bond dynamics of water molecules in lipid-bilayer systems at varying temperatures

Abstract

Cholesterol (Chol) plays a crucial role in shaping the intricate physicochemical attributes of biomembranes, exerting considerable influence on water molecules proximal to the membrane interface. In this study, we conducted molecular dynamics simulations on the bilayers of two lipid species, dipalmitoyl phosphatidylcholine (DPPC) and palmitoyl sphingomyelin (PSM); they are distinct with respect to the structures of the hydrogen-bond (H-bond) acceptors. Our investigation focuses on the dynamic properties and H-bonds of water molecules in the lipid-membrane systems, with particular emphasis on the influence of Chol at varying temperatures. Notably, in the gel phase at 303 K, the presence of Chol extends the lifetimes of H-bonds of the oxygen atoms acting as H-bond acceptors within DPPC with water molecules by a factor of 1.5 to 2.5. In the liquid-crystalline phase at 323 K, on the other hand, H-bonding dynamics with lipid membranes remain largely unaffected by Chol. This observed shift in H-bonding states serves as a crucial key to unraveling the subtle control mechanisms governing water dynamics in lipid-membrane systems.

I Introduction

Lipid bilayers, serving as the fundamental architectural frameworks of biological membranes, form stable aggregates through the amphiphilic effect inherent to lipid molecules. The attributes of lipid bilayers vary diversely with the specific type and composition of lipid molecules, giving rise to distinctive structures, such as gel and liquid-crystalline phases. [1] The interactions with water are also a key to self-organization,[2] and water properties are affected by the states of lipid molecules in turn.

Cholesterol (Chol), extensively investigated as a constituent of bio-related membranes, exhibits significant influence on structures of lipid bilayers. [3, 4] In particular, Chol enhances the packing density and rigidity of the lipid, thereby modulating membrane fluidity. Furthermore, water molecules play a crucial role in the structure and function of biological membranes, influencing electrostatic properties, solute exchange, and protein function. [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16] Thus, it is imperative to elucidate the structure and dynamics of water molecules at the membrane interface, which is expected to differ from those in the bulk owing to the interaction with hydrophilic groups on the lipid head.

Molecular dynamics (MD) simulation is a valuable tool for investigating lipid bilayers, providing molecular-level insights into not only lipid properties but also their interactions with other molecules.  [17, 18, 19, 20, 21, 22, 23, 24, 25, 26] Numerous investigations have also been conducted for water in the interface region, encompassing the distribution of water molecules, reorientation dynamics, mean square displacement, and hydrogen-bond (H-bond) dynamics. [27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40] Specifically, the slowdown of the H-bond dynamics from the bulk to the center of the membrane has been demonstrated. [39, 40] In addition, MD simulations have been used to study the structure and dynamics of lipid bilayers containing phospholipids and Chol. [41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54] These simulations have provided insights into the influence of Chol on a variety of phospholipids differing in the headgroup and tail.

Despite the numerous studies mentioned above, comprehending the water state proximal to lipid membranes, as well as understanding its connection with the membrane state in the presence of Chol, remains a significant problem. Interestingly, experimental observations have indicated that Chol was found to accelerate the water dynamics in the dipalmitoyl phosphatidylcholine (DPPC) membrane interface. [55, 56] Conversely, it has been found that water dynamics decelerate within the interior region of lipid bilayers with increasing Chol concentration. [55]

MD simulations have elucidated that the acceleration of water dynamics at the interface, particularly notable at high Chol concentrations up to 50%, arises from the inhibition of H-bonds between two oxygen atoms of lipid molecules. [49] A more recent MD study conducted a detailed analysis of the H-bond network of water within the DPPC membrane in the presence of Chol. The results unveiled that Chol fosters more bulk-like water at the membrane interface, leading to increased local water density and accelerated water dynamics. [57]

In this study, we conducted MD simulations of two types of lipid bilayers comprising of DPPC and palmitoyl sphingomyelin (PSM) with the presence of Chol. While the Chol concentrations investigated were 0 and 10%, the temperature effect was examined at 303 K and 323 K. The DPPC and PSM membranes are at the gel and liquid-crystalline phases at 303 K and 323 K, respectively. We investigate the microscopic hydration structure and dynamics by considering acceptor sites of lipid molecules, which form H-bonds with the hydrogen atoms of water molecules. Between DPPC and PSM, the choline and phosphate groups are identical, as shown in Fig. 1. There are differences in the degree of carbon chain saturation and the functional group acting as the H-bonding site. Thus, our MD investigations provide insights into H-bonds influenced by Chol, taking into account the molecular structures of the lipids and the environmental effects from the membrane composition and temperature.

II Simulation details

Figure 1: Structures of the lipid and Chol molecules studied in this paper.

Table 1: Numbers of lipid (DPPC or PSM), Chol, and water molecules in mixture and pure membrane systems.

	mixture	pure
DPPC / PSM	200	200
Chol	22	-
Water	22000	20000

Figure 2: Number density distributions of lipid carbon chain (tail), water molecule oxygen ($\mathrm{O}^{\mathrm{w}}$), and Chol along the $z$-direction. Solid lines represent pure membrane systems, while dashed lines represent systems containing Chol. (a) and (b) correspond to DPPC, and (c) and (d) to PSM.

The structures of the lipid molecules, DPPC, PSM, and Chol are depicted in Fig. 1. The hydrophilic moiety is common between DPPC and PSM, while they are different for the hydrophobic portion and the distributions of oxygen and nitrogen atoms.

The lipid bilayer system was constructed using CHARMM-GUI, [58, 59, 60, 61, 62] incorporating 200 lipid molecules with 10% of Chol if present, as listed in Table 1. For each lipid and Chol composition, 100 water molecules per molecule of lipid and Chol were added to create the lipid bilayer system. Three different initial configurations were prepared for each composition, employing the CHARMM36 force field for DPPC, PSM, and Chol, [63] and the CHARMM-compatible TIP3P model for water molecules. [64] All the MD simulations were performed using Gromacs 2022.4.  [65]

Figure 3: (a) Schematic illustration of lipid phosphorus atom position relative to the instantaneous average in the upper and lower leaflet, denoted as $z^{\mathrm{upper}}_{G}(t)$ and $z^{\mathrm{lower}}_{G}(t)$, respectively. $z^{\mathrm{P}}_{j}(t)-z^{\mathrm{upper}}_{G}(t)$ and $z^{\mathrm{P}}_{j}(t)-z^{\mathrm{lower}}_{G}(t)$ represent the distances of the P atom of lipid $j$, $z^{\mathrm{P}}_{j}(t)$, from $z^{\mathrm{upper}}_{G}(t)$ and $z^{\mathrm{lower}}_{G}(t)$, respectively, along the $z$-direction. (b) and (c) depict the distributions of $P(|\Delta z|)$ for DPPC and PSM, respectively. Snapshots in each panel are taken at 303 K (Left: with Chol, Right: pure).

The equilibration process is described in Table S1 of the supplementary material. In accordance with CHARMM-GUI guidelines, the process involved gradually relaxing restraints imposed on the phosphorus atom and the chiral carbon center of the lipid molecule (Nos. 1-6). The constants of the restraining forces on the $z$-coordinate of the phosphorus atom and on dihedral angle concerning the asymmetric center and double bond are denoted as $k_{z}$ and $k_{\mathrm{dih}}$, respectively, in Table S1. Subsequently, further equilibration steps were carried out (Nos. 7-11); the computational stability was checked in each of nos. 7-11 and the MD length was gradually increased from no. 7 to 11. Finally, three production runs under $NPT$ conditions for 10 ns each were performed (No. 12). To examine how the effect of Chol depends on the phase of the lipid membrane (gel vs liquid-crystalline), MD simulations were conducted at 303 K and 323 K for each system. The coordinate system was set so that the $z$-axis is normal to the membrane surface, which spans over the $x$- and $y$-directions.

To confirm the adequacy of the equilibration process, we examined the time evolution of surface area $S$ in the $x$-$y$ plane. Figures S1 and S2 of the supplementary material illustrate these results during the 3 $\mu$s equilibration at 303 K and 323 K, respectively. While noticeable fluctuations are observed around 1.5 $\mu$s in some systems, the area $S$ converges to a stable value at approximately 3 $\mu$s across all systems. Consequently, equilibration for 3 $\mu$s is considered adequate, and a production run was carried out after this equilibration period.

Figure 4: (a) Schematic illustration of calculations of water molecule distribution $\rho(z^{\prime})$. (b) and (c) show the ratios of the water molecule distribution $\rho(z^{\prime})$, to the number density of bulk water, $\rho_{\mathrm{bulk}}$, for DPPC and PSM, respectively. The values of $\rho_{\mathrm{bulk}}$, determined from $\rho(z)$ as illustrated in Fig. 2, are 33.50 nm^-3 and 32.95 nm^-3 for DPPC and PSM, respectively. (d) and (e) depict the average numbers of H-bonds using the same $z^{\prime}$-axis as (b) and (c), for DPPC and PSM, respectively. The dashed lines at -0.25 nm and 0.8 nm define the boundaries between region 1 and 2 and between region 2 and 3, respectively.

III Results and discussion

III.1 Density distributions of lipid, water, and Chol

Initially, we analyzed the structure of the lipid bilayer, as well as the configurations of Chol and water. Figure 2 illustrates the number density distributions $\rho(z)$ of the lipid carbon chain (tail), water molecular oxygen ($\mathrm{O}^{\mathrm{w}}$), and center of mass of Chol along the $z$-direction. Here, $z$ denotes the distance from the bilayer center, and the center of the bilayer is the center of mass of the lipid molecules. Note that the number density distributions of the upper and lower leaflets were found to be coincident with each other within the margin of errors (data not shown). $\rho(z)$ in Fig. 2 represents the average of the profiles of both leaflets. As depicted in Fig. 2, Chol is predominantly situated near the carbon chain of the lipid, leading to a broadening of the tail distribution along the $z$-direction, particularly evident at 303 K in DPPC. To examine this in detail, the number density distributions of $\mathrm{O}^{\mathrm{D7}}$ in DPPC and $\mathrm{O}^{\mathrm{S5}}$ in PSM are plotted in Fig. S3 of the supplementary material. These oxygen atoms are part of the hydrophilic functional groups within each lipid molecule, which are located in the innermost part of the lipid membranes (see Fig. 1). Figure S3 of the supplementary material demonstrates that the tail of water density distribution overlaps with those of $\mathrm{O}^{\mathrm{D7}}$ in DPPC and $\mathrm{O}^{\mathrm{S5}}$ in PSM. In the presence of Chol, particularly for DPPC at 303 K, the distribution of $\mathrm{O}^{\mathrm{D7}}$ becomes narrower, hindering the penetration of water molecules into the membrane’s inner regions, as observed in Fig. 2(a). Additionally, Fig. S4 of the supplementary material displays the number density distributions of nitrogen (N⁺) in the choline group and phosphorus (P) atoms along the $z$-direction. The peak intensities of distributions for N⁺ and P atoms were enhanced, particularly at 303 K in DPPC.

III.2 Fluctuation of the membrane interface

The surfaces of lipid membranes are soft and fluctuate with time. To elucidate the structure of the membrane interface, our focus was directed towards the lipid head, where we examined the distribution of lipid phosphorus atom position relative to the instantaneous interface defined below. The fluctuation of the interface between the membrane and water is seen evidently by employing the instantaneous interface method. [66]

We assign each lipid to either the upper or lower leaflet at each time. Here, $N^{\mathrm{upper}}_{\mathrm{p}}$ denotes the number of lipid molecules in the upper leaflet, $z_{j}^{\mathrm{P}}\left(t\right)$ is the $z$-coordinate at time $t$ of the $j$th lipid molecule, and $z^{\mathrm{upper}}_{G}(t)=(1/{N^{\mathrm{upper}}_{\mathrm{p}}})\sum_{j\in{\mathrm{upper}}}z^{\mathrm{P}}_{j}(t)$ represents the average of the $z$-coordinates of phosphorus atoms in the upper leaflet of the lipid bilayer at time $t$. Similarly, $N^{\mathrm{lower}}_{\mathrm{p}}$ and $z^{\mathrm{lower}}_{G}(t)$ can be computed for the lower leaflet.

As shown in Fig. 3(a), the deviation of the $z$-coordinate of the phosphorus atom of lipid $j$ from $z^{\mathrm{upper}}_{G}(t)$ is expressed as $\Delta z^{\mathrm{upper}}_{j}(t)=z^{\mathrm{P}}_{j}(t)-z^{\mathrm{upper}}_{G}(t)$ in the upper leaflet. Similarly, for lipids in the lower leaflet, $\Delta z^{\mathrm{lower}}_{j}(t)$ is defined. Then, the time-averaged distribution function of the absolute values $\left|\Delta z\right|$ of $\Delta z^{\mathrm{upper}}_{j}(t)$ and $\Delta z^{\mathrm{lower}}_{j}(t)$ can be assessed and is denoted as $P(|\Delta z|)$.

Figures 3(b) and 3(c) illustrate the results of $P(|\Delta z|)$ for DPPC and PSM, respectively. In both DPPC and PSM, Chol does not exert discernible effects on $P(|\Delta z|)$ at 323 K. Snapshots captured at 323 K are depicted in Fig. S5 of the supplementary material, revealing a disordered orientation of carbon chains within lipids. This observation signifies a high degree of membrane fluidity, with weak influence of Chol. In contrast, at 303 K, the distribution of $P(|\Delta z|)$ is broader in the pure lipid membrane systems, suggesting that Chol enhances membrane stability and maintains the interface position. The effect of Chol to suppress the interface fluctuations is particularly evident in DPPC, as illustrated in Fig. 3(b) (see also snapshots captured at 303 K in the insets of Figs. 3(b) and 3(c)).

III.3 Classification of water molecules

The analysis of water molecule distribution near the rugged membrane interface was conducted. A precise description of the local water distribution relative to the lipid interface was proposed [22, 23] The location of a water molecule is provided by $z^{\prime}$, which is defined as the distance from the interface using Voronoi tessellation. Unlike $z$, $z^{\prime}$ takes into account the effects of the fluctuation of the lipid/water interfaces, and the distribution $\rho(z^{\prime})$ characterizes the layered structures of water molecules.

We propose a simpler method closely resembling Voronoi tessellation. The schematic illustration of the method is described in Fig. 4(a) and the detail is provided as follows: Step 1: Project the oxygen atoms of water molecules and the lipid phosphorus atoms onto the $x$-$y$ plane. Calculate the distance $d_{ij}=|\bm{r}^{\mathrm{O}}_{i}(t)-\bm{r}^{\mathrm{P}}_{j}(t)|$ between the oxygen atom of water molecule $i$ and the phosphorus atom $j$ in the $x$-$y$ plane. Identify the lipid phosphorus atom $i^{\prime}$ that gives the smallest $d_{ij}$ for water molecule $i$. Step 2: Determine $z^{\mathrm{O}}_{i}(t)$ as the $z$-coordinate of oxygen atom of water molecule $i$ and $z^{\mathrm{P}}_{i^{\prime}}(t)$ as the $z$-coordinate of the phosphorus atom $i^{\prime}$ nearest to $i$ in the $x$-$y$ plane. Define the water molecule density of $z^{\prime}=z^{\mathrm{O}}_{i}-z^{\mathrm{P}}_{i^{\prime}}$ as $\rho(z^{\prime})$. Note that the negative $z^{\prime}$ indicates that the water molecule is located in more inner positions of the membrane than the phosphorus atom.

Figure 5: (a) Schematic illustration of $r_{\mathrm{oo}}$ and $\beta$ for H-bond formed between the oxygen atom ($\mathrm{O}^{\mathrm{D5}}$) of a functional group within DPPC and a water molecule. (b)-(i) 2D PMF, $W(r_{\mathrm{oo}},\beta)$, between water oxygen ($\mathrm{O}^{\mathrm{w}}$) and each acceptor oxygen in DPPC with Chol at 303 K. Black points represent saddle points.

Figures 4(b) and 4(c) illustrates the ratio of $\rho(z^{\prime})$ to the number density of bulk water, $\rho_{\mathrm{bulk}}$, for DPPC and PSM, respectively. In contrast to the water molecule density profile $\rho(z)$ in Fig. 2, which is influenced by the instantaneous fluctuations of lipid membranes, $\rho(z^{\prime})$ can be a more faithful representation of the water distribution near the rugged surface.

From the profile of $\rho(z^{\prime})$, the water molecules can be categorized into three regions (regions 1-3), as depicted in Figs. 4(b) and 4(c). Specifically, region 1 represents the region inside the membrane at $z^{\prime}<$ -0.25 nm, region 2 denotes the interface region at -0.25 nm $<z^{\prime}<$ 0.8 nm, and region 3 encompasses the bulk region at $z^{\prime}>$ 0.8 nm. These classifications align with previous studies. [22, 23] Note that the number of DPPC molecules was 128 in Ref. 22, which is slightly smaller than that of our system. Furthermore, similar results of $\rho(z^{\prime})$ were reported by Elola et al, where 968 DPPC molecules with the united-atom model were simulated. [49]

The water content in region 1 decreases progressively towards the center of the membrane. In the interface region (region 2), a minimum was observed near $z^{\prime}=0$ nm, with a peak occurring around $z^{\prime}=0.4$ nm, for both DPPC and PSM. This peak stems from the tendency for H-bond formation around phosphate groups. A further elucidation on the H-bond rearrangement will be provided in subsequent Sec. III.4. Remarkably, as shown in Fig. 4(b), at 303 K in DPPC, the water content in the interface region (region 2) is larger in the pure lipid membrane system than in the presence of Chol. This observation aligns with the variation in $\rho(z)$ of water molecules due to Chol, as shown in Fig. 2(a). The pronounced stabilization of the DPPC membrane by Chol at 303 K highlights the significant impact on the hydration structure near the interface. On the contrary, at 323 K in DPPC, but the effect of Chol is less significant. Moreover, for PSM, the impacts of both temperature variations and Chol on $\rho(z^{\prime})$ are not appreciable, as demonstrated in Fig. 4(c).

Figure 6: (a) Schematic illustration of $r_{\mathrm{oo}}$ and $\beta$ for H-bond formed between the oxygen atom ($\mathrm{O}^{\mathrm{S5}}$) of a functional group within PSM and a water molecule. (b)-(i) 2D PMF, $W(r_{\mathrm{oo}},\beta)$, between water oxygen ($\mathrm{O}^{\mathrm{w}}$) and each acceptor oxygen in PSM with Chol at 303 K. Black points represent saddle points.

III.4 H-bond arrangement

The H-bonding states at each of the donor and acceptor sites were analyzed. When investigating H-bond state in MD simulations, a commonly employed approach involves applying a geometric criterion to identify an H-bond between two water molecules. The predominant definition often adopts the distance between oxygen atoms (referred to as $r_{\mathrm{oo}}$) and the angle formed by the oxygen atom and the oxygen-hydrogen bond (referred to as $\beta$) within a water dimer. [67, 68, 69]

A more comprehensive understanding of the H-bond state can be obtained by analyzing the distribution function of $r_{\mathrm{oo}}$ and $\beta$, denoted as $g(r_{\mathrm{oo}},\beta)$. [70, 71, 72, 73] In this context, $2\pi\rho r_{\mathrm{oo}}^{2}\sin\beta g(r_{\mathrm{oo}},\beta)\differential{r_{\mathrm{oo}}}\differential{\beta}$ represents the average number of oxygen atoms acting as H-bond acceptors within the partial spherical shell volume characterized by $\mathrm{d}r_{\mathrm{oo}}$ and $\mathrm{d}\beta$ at the position $(r_{\mathrm{oo}},\beta)$, with the average number density of water molecules, $\rho$. The logarithm form $W(r_{\mathrm{oo}},\beta)=-k_{\mathrm{B}}T\ln{g(r_{\mathrm{oo}},\beta)}$ can be interpreted as the two-dimensional potential of mean force (2D PMF). For reference, the 2D PMF $W(r_{\mathrm{oo}},\beta)$ of bulk water at 303 K and 323 K with a density of 1 g/cm³ is depicted in Fig. S6 of the supplementary material. The temperature-independent energetically stable state is characterized by $r_{\mathrm{oo}}<0.35$ nm and $0^{\circ}<\beta<30^{\circ}$, which can be considered indicative of H-bond state.

Figures 5(a) and 6(a) provide schematic illustrations of $r_{\mathrm{oo}}$ and $\beta$ for H-bond formed between a water molecule and the oxygen atom of a functional group within DPPC and PSM, respectively. Given the presence of H-bonds between water molecules and those with acceptors within lipid molecules, the analysis of 2D PMF was conducted for water molecules and potential acceptors, including the oxygen atoms in DPPC and the oxygen and nitrogen atoms in PSM. Additionally, in the presence of Chol, the analysis was also performed for H-bonds between water molecule and oxygen atom of hydroxy group, denoted as $\mathrm{O}^{\mathrm{Chol}}$. A similar 2D PMF analysis was previously done in polymer-water mixtures. [74]

Figure 5 illustrates the 2D PMF, $W(r_{\mathrm{oo}},\beta)$, representing the interaction between water molecule as donors and DPPC oxygen atoms as acceptors at 303 K in the presence of Chol. Other results at 303 K in the absence of Chol and at 323 K both in the presence and absence of Chol are displayed in Figs. S7-S9 of the supplementary material. Similarly, in Fig. 6 and Figs. S10-S12 of the supplementary material provide the 2D PMF, $W(r_{\mathrm{oo}},\beta)$, for PSM systems. Note that the 2D PMF between water molecules is omitted since the overall profile remains unchanged for both DPPC and PSM, when compared to that of bulk water (see Fig. S6 of the supplementary material). Based on the 2D PMF analysis, we identify potential H-bond acceptors as $\mathrm{O}^{\mathrm{1}}$, $\mathrm{O}^{\mathrm{2}}$, $\mathrm{O}^{\mathrm{3}}$, $\mathrm{O}^{\mathrm{D5}}$, $\mathrm{O}^{\mathrm{D7}}$ and $\mathrm{O}^{\mathrm{Chol}}$ for DPPC, and $\mathrm{O}^{\mathrm{1}}$, $\mathrm{O}^{\mathrm{3}}$, $\mathrm{O}^{\mathrm{S4}}$, $\mathrm{O}^{\mathrm{S5}}$, and $\mathrm{O}^{\mathrm{Chol}}$ for PSM, respectively. See Fig. 1 for the notations of the acceptor sites in the lipids. As illustrated in Figs. 5 and 6 and Figs. S7-S12 of the supplementary material, the H-bond region is characterized by $r_{\mathrm{oo}}<3.5$ nm and $0^{\circ}<\beta<30^{\circ}$ in the 2D PMF, irrespective of the acceptor. Furthermore, the H-bond regions remain unchanged regardless of the presence of Chol or variations in temperature for both DPPC and PSM. Nevertheless, specific oxygen atoms such as $\mathrm{O}^{\mathrm{1}}$, $\mathrm{O}^{\mathrm{D5}}$, and $\mathrm{O}^{\mathrm{D7}}$ in DPPC, and $\mathrm{O}^{\mathrm{1}}$ and $\mathrm{O}^{\mathrm{S5}}$ in PSM, which possess higher negative charges, form more energetically stable states compared to other acceptors. In contrast, oxygen atoms $\mathrm{O}^{\mathrm{D4}}$ and $\mathrm{O}^{\mathrm{D6}}$ in DPPC, and oxygen atom $\mathrm{O}^{\mathrm{2}}$ and nitrogen (N) in PSM, form a second coordination region outside the defined H-bond region. Consequently, these are excluded from further H-bond analysis due to their indeterminate bond characteristics.

Figures 4(d) and 4(e) illustrate the distributions of the average number of H-bonds formed by water molecules at each position, using the $z^{\prime}$-axis corresponding to $\rho(z^{\prime})$, for DPPC and PSM, respectively. The average number of H-bonds in region 3 converged to 3.33 at 303 K and 3.23 at 323 K, respectively, corresponding to those observed in bulk water at each temperature. In region 2, the average number of H-bonds reaches a maximum value higher than the average observed in the bulk, gradually decreasing towards the interior of the bilayer. This peak position corresponds to that in $\rho(z^{\prime})$, as observed in Figs. 4(b) and 4(c). These observations suggest that around phosphate groups, the oxygen atoms within the lipid head group, such as $\mathrm{O}^{\mathrm{1}}$, $\mathrm{O}^{\mathrm{2}}$, and $\mathrm{O}^{\mathrm{3}}$ in DPPC, and $\mathrm{O}^{\mathrm{1}}$ and $\mathrm{O}^{\mathrm{3}}$ in PSM, act as acceptors, promoting the formation of H-bonds with water molecules. Remarkably, at 303 K, the average number of H-bonds in region 1 increases with removal of Chol for both DPPC and PSM, with this trend being particularly notable in the DPPC system. However, at 323 K, the average number of H-bonds remains unchanged regardless of the presence or absence of Chol, for both DPPC and PSM. These findings indicate variations in membrane structure induced by Chol impact the propensity for H-bond formation within the membrane.

Figure 7: Conditional probability $C_{2,j}(t)$, representing transition dynamics from region $2$ at the initial time $t=0$ to either regions 1 or 3 at time $t$ or remaining within the same region 2 in the time interval $t$ [(a) DPPC at 303 K, (b) DPPC at 323 K, (c) PSM at 303 K, and (d) PSM at 323 K].

III.5 Water molecule rearrangement dynamics

We explore the transition dynamics of water molecules among the three regions. We address the dynamics of transition by defining $C_{i,j}(t)$ ($i\neq j$) as the conditional probability that when a water molecule is in region $i$ at time 0, it visits region $j$ with the number of passing the boundaries of the regions being unity by time $t$. Further, $C_{i,i}(t)$ is the probability that the water molecule stays in region $i$ without visiting the other regions during the time interval between 0 and $t$. Note that the summation over all possible $j$ states ensures $\sum_{j}C_{i,j}(t)=1$, conserving the number of water molecules. In practice, the trajectories of water molecules are continuously monitored from $t=0$, tracking the subsequent transitions from region $i$ to region $j$ at each time $t$.

Figure 7 presents the results of $C_{2,j}(t)$, illustrating the transition dynamics of water molecules originating from region 2 at $t=0$. Note that the sum $C_{2,1}(t)+C_{2,2}(t)+C_{2,3}(t)=1$ holds at all times $t$, as explained in the definition of $C_{i,j}(t)$. Except for DPPC at 303 K, $C_{2,j}(t)$ is not affected by the presence or absence of Chol and the transition rates from region 2 to regions 1 and 3 are common between the systems with and without Chol. In contrast, for DPPC at 303 K, the decay of $C_{2,2}(t)$ exhibits a slower rate in the presence of Chol. Furthermore, Chol alters the fraction of water molecules transitioning to their respective destination. Specifically, the population of water molecules transitioning from region 2 to region 1 decreases by approximately 5% in the presence of Chol, while the transition to region 3 increases by a similar proportion. At 323 K, the saturated values of $C_{2,1}(t)$ and $C_{2,3}(t)$ resemble those of DPPC without Chol at 303 K. These observations indicate the significant impact of membrane structure variations on water molecule dynamics. For DPPC at 303 K, in particular, Chol enhance the tendency of keeping water molecules in the interface region (region 2) and relocating them towards the bulk region (region 3). This suggests the reduction of water molecule exchanges between the interface region (region 2) and the inner side of the membrane (region 1).

Figure S13 of the supplementary material illustrates the results of $C_{1,j}(t)$, which represent the transition dynamics from region $1$ at the initial time $t=0$ to regions 2 or remaining within the same region $1$ at subsequent time $t$. Here, $C_{1,3}(t)$ is excluded given that the transition from region 1 to region 3 inevitably passes through region 2, ensuring the relationship, $C_{1,1}(t)+C_{1,2}(t)=1$. Interestingly, the transition from region 1 to region 2 exhibits slower dynamics in PSM compared to in DPPC at 303 K and 323 K. This observation may be linked to the shorter H-bond lifetime $\tau_{\mathrm{HB}}$ of water molecules in the PSM system than with DPPC, as elucidated in the subsequent Sec. III.6. In addition, Chol further retards these dynamics, particularly evident at 303 K for both DPPC and PSM. This observation suggests that Chol, situated within the membrane interior, exerts a notable influence on the dynamics of water molecules within region 1.

Figure 8: Dependency of the ratio of H-bond lifetimes, $\tau_{\mathrm{HB}}$, with and without Chol, denoted as $\tau_{\mathrm{HB(Chol)}}/\tau_{\mathrm{HB(pure)}}$ for (a) DPPC and (b) PSM.

III.6 Chol influence on H-bond lifetime

Finally, we conducted an analysis of the H-bonding dynamics involving lipid molecules, Chol, and water molecules to elucidate the timescale of H-bond lifetime, by focusing on the acceptor oxygen atoms, such as $\mathrm{O}^{\mathrm{1}}$, $\mathrm{O}^{\mathrm{2}}$, $\mathrm{O}^{\mathrm{3}}$, $\mathrm{O}^{\mathrm{D5}}$, and $\mathrm{O}^{\mathrm{D7}}$ in DPPC, and $\mathrm{O}^{\mathrm{1}}$, $\mathrm{O}^{\mathrm{3}}$, $\mathrm{O}^{\mathrm{S4}}$, and $\mathrm{O}^{\mathrm{S5}}$ in PSM, and $\mathrm{O}^{\mathrm{Chol}}$ within Chol. The H-bond time correlation function $P_{\mathrm{HB}}(t)$ is defined as

$$P_{\mathrm{HB}}(t)=\dfrac{\left\langle h_{i,j}(t)h_{i,j}(0)\right\rangle}{\left\langle h_{i,j}(0)\right\rangle},$$ (1)

where $h_{i,j}(t)$ equals 1 if water molecule $i$ is H-bonded with acceptor oxygen $j$ at time $t$, otherwise 0. [75, 68, 67] We computed $P_{\mathrm{HB}}(t)$ using the Monte-Carlo bootstrap method, which employs a non-parametric approach to statistical inference. [76]

Figures S14 and S15 of the supplementary material show the $P_{\mathrm{HB}}(t)$ results for acceptor oxygen atoms in the DPPC and PSM systems, respectively. Notably, as the temperature decreases, the decay of $P_{\mathrm{HB}}(t)$ slows down for each acceptor oxygen atom, with a significant effect observed for $\mathrm{O}^{\mathrm{Chol}}$ in Chol. However, the impact of temperature variation on H-bond breakages between water molecules is negligible for both DPPC and PSM, owing to the abundance of H-bonding partners of $\mathrm{O}^{\mathrm{w}}$ in the bulk. Furthermore, the influence of Chol is more pronounced in DPPC compared to PSM, particularly at 303 K. As illustrated in Figs. S15 and S16, the H-bond correlation between $\mathrm{O}^{\mathrm{w}}$ and $\mathrm{O}^{\mathrm{Chol}}$ exhibits a slower dynamics in PSM compared to DPPC.

The H-bond time correlation function, $P_{\mathrm{HB}}(t)$, is approximated using the Kohlrausch–Williams–Watts (KWW) function, $P_{\mathrm{HB}}(t)\simeq\exp[-(t/\tau_{\mathrm{KWW}})^{\beta_{\mathrm{KWW}}}]$. The fitting results are depicted in Figs. S14 and S15 of the supplementary material for DPPC and PSM, respectively. The H-bond lifetime $\tau_{\mathrm{HB}}$ is evaluated by integrating $P_{\mathrm{HB}}(t)$, yielding

$$\tau_{\mathrm{HB}}=\int_{0}^{\infty}P_{\mathrm{HB}}(t)\differential{t}=\frac{\tau_{\mathrm{KWW}}}{\beta_{\mathrm{KWW}}}\Gamma\left(\frac{1}{\beta_{\mathrm{KWW}}}\right),$$ (2)

where $\Gamma(\cdots)$ denotes the Gamma function. The raw data of $\tau_{\mathrm{KWW}}$, $\beta_{\mathrm{KWW}}$, and $\tau_{\mathrm{HB}}$ are summarized in Tables S2-S7 of the supplementary material for both DPPC and PSM. To highlight the influence of Chol on H-bond lifetime $\tau_{\mathrm{HB}}$, the dependency of the ratio between $\tau_{\mathrm{HB}}$ with and without Chol, denoted as $\tau_{\mathrm{HB(Chol)}}/\tau_{\mathrm{HB(pure)}}$, on acceptor oxygen at 303 K and 323 K is shown in Fig. 8.

Figure 8(a) illustrates the ratio $\tau_{\mathrm{HB(Chol)}}/\tau_{\mathrm{HB(pure)}}$, ranging from 1.5 to 2.5, excluding $\mathrm{O}^{\mathrm{w}}$, at 303 K for DPPC. Moreover, $\tau_{\mathrm{HB}}$ becomes large with the internal oxygen atoms within the membrane, such as $\mathrm{O}^{\mathrm{D5}}$ and $\mathrm{O}^{\mathrm{D7}}$. In contrast, $\mathrm{O}^{\mathrm{1}}$ and $\mathrm{O}^{\mathrm{3}}$ exhibits a relatively faster H-bond lifetime, showing that the H-bond dynamics is less susceptible to the presence or absence of Chol near the aqueous region. However, in the case of PSM, the influence of cholesterol on H-bond lifetime is limited, except for $\mathrm{O}^{\mathrm{S5}}$, at both 303 K and 323 K, as observed in Fig. 8(b). The slower dynamics observed for $\mathrm{O}^{\mathrm{S5}}$ in PSM can be attributed to its position as the innermost oxygen atom within the membrane, rendering it more susceptible to Chol than other oxygen atoms.

Finally, we examine the energetic aspect of the H-bond lifetime $\tau_{\mathrm{HB}}$. Specifically, the activation energy $\Delta G$ of the H-bond is estimated by the free energy difference between the most stable and saddle points on the 2D PMF, $W(r_{\mathrm{OO}},\beta)$ (see Figs. 5 and  6, and Figs. S7-S12 of the supplementary material). We plotted the relationship between $\tau_{\mathrm{HB}}$ and $\Delta G/k_{\mathrm{B}}T$ in Fig. S16 of the supplementary material. Assuming the Arrhenius equation, $\tau_{\text{HB}}=A\exp({\Delta G}/{k_{\text{B}}T})$, we determined the constant $A$ by fitting. The straight lines with $A=0.183$ ps and $A=0.385$ ps are shown for DPPC and PSM systems, respectively, in Figs. S16(a) and S16(b) of the supplementary material. The positive correlation between $\tau_{\text{HB}}$ and $\Delta G/k_{\mathrm{B}}T$ is demonstrated for both DPPC and PSM systems, although the deviation from the Arrhenius equation is noticeable. This result suggests that a more accurate energetic description of H-bond breakage involving a large number of molecules necessitates additional detailed variables beyond the two variables of distance $r_{\mathrm{OO}}$ and angle $\beta$.

IV Conclusions

In this study, we employed MD simulations to investigate the influence of Chol on water molecule behavior within lipid membranes, with a specific focus on systems comprising DPPC and PSM. While lipid membrane structures are not susceptible to the presence of Chol at 323 K, Chol at 303 K serves to stabilize carbon chains, thereby reducing structural fluctuations at the membrane interface, particularly for DPPC.

The spatial distribution of water molecules surrounding the membrane can be classified into three distinct regions: the membrane interior, the interface, and the bulk. The transition of water from the interface to the bulk is facilitated by Chol for the DPPC system at 303 K. In contrast, the presence of Chol induces entrapment of water molecules within the membrane, leading to reduced rates of transition to the interface region from the interior region.

Our exploration into the dynamic attributes of water molecules in lipid-membrane systems, considering the influence of Chol and temperature variations has yielded insights into the intricate interplay at the membrane interface. DPPC is more susceptible to modifications at 303 K, significantly influencing H-bond dynamics within the membrane. Specifically, at 303 K in DPPC, Chol was found to markedly increase the H-bonding lifetime, particularly impacting internal oxygen atoms.

It is important to note the discrepancy in Chol density between our study, which utilized up to 10%, and 50% employed by Elola et al. [49] This emphasizes the need to further investigate H-bonding dynamics in lipid-membrane systems under conditions compatible to real Chol contents in future research. While Elola et al. reported accelerated water dynamics near the interface region, our findings have not directly corroborated these observations. However, we observed a notable migration tendency of water molecules from region 2 (interface region) to region 3 (bulk region) in DPPC in the presence of Chol, as illustrated in Fig. 7(a). Therefore, future studies should focus on systematically varying both temperature and Chol content to provide a comprehensive understanding of H-bonding dynamics.

Supplementary material

The supplementary material include equilibration scheme of MD simulations (Table S1), time evolution of surface area $S$ in the $x$-$y$ plane during the equilibration (Fig. S1 and S2), number density distribution of $\mathrm{O}^{\mathrm{D7}}$ in DPPC and $\mathrm{O}^{\mathrm{S5}}$ in PSM along the $z$-direction (Fig. S3), number density distribution of nitrogen and phosphorus atoms along the $z$-direction (Fig. S4), MD snapshots taken at 323 K (Fig. S5), 2D PMF between water molecules in bulk water (Fig. S6), 2D PMF between water oxygen and acceptor oxygen atoms in lipid molecules (Fig. S7-S12), conditional probability $C_{1,j}(t)$ (Fig. S13), H-bond time correlation function $P_{\mathrm{HB}}(t)$ (Figs. S14 and S15), raw data of $\tau_{\mathrm{KWW}}$, $\beta_{\mathrm{KWW}}$, and $\tau_{\mathrm{HB}}$ (Table S2-S7), and the relationship between the H-bond lifetime $\tau_{\mathrm{HB}}$ and activation energy normalized by thermal energy, $\Delta G/k_{\mathrm{B}}T$ (Fig. S16).

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Data availability statement

The program codes and data that support the findings of this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.11273701. Further data are available from the corresponding author upon reasonable request.

Supplementary Material

Influence of cholesterol on hydrogen-bond dynamics of water molecules in lipid-bilayer systems at varying temperatures

Kokoro Shikata, Kento Kasahara, Nozomi Morishita Watanabe, Hiroshi Umakoshi, Kang Kim, and Nobuyuki Matubayasi

Division of Chemical Engineering, Department of Materials Engineering Science, Graduate School of Engineering Science, Osaka University, Toyonaka, Osaka 560-8531, Japan

Table S1: Equilibration scheme.

No.	process	$dt$ [fs]	time [ns]	integrater	$k_{z}~{}[\mathrm{kJ/mol\,nm^{2}}]$	$k_{\mathrm{dih}}~{}\mathrm{[kJ/mol\,rad^{2}]}$
1	Energy Minimization	-	5000 (steps)	steepest descent	1000	1000
2	$NVT$	1	0.125	Leap-Flog	1000	1000
3	$NVT$	1	0.125	Leap-Flog	400	400
4	$NPT$	1	0.125	Leap-Flog	400	200
5	$NPT$	2	0.5	Leap-Flog	200	200
6	$NPT$	2	0.5	Leap-Flog	40	100
7	$NPT$	2	0.5	Leap-Flog	0	0
8	$NPT$	2	10	Leap-Flog	0	0
9	$NPT$	2	100	Leap-Flog	0	0
10	$NPT$	2	500	Leap-Flog	0	0
11	$NPT$	2	3000	Leap-Flog	0	0
12	Production ($NPT$)	2	10	Leap-Flog	0	0

Figure S1: Time evolution of surface area $S$ in the $x$-$y$ plane during the equilibration at 303K. (a) and (b): DPPC; (c) and (d): PSM. (a) and (b) refer to the systems with Chol, while (c) and (d) are for those without Chol. Systems 01, 02, and 03 correspond to three distinct runs with different initial configurations.

Figure S2: Time evolution of surface area $S$ in the $x$-$y$ plane during the equilibration at 323K. (a) and (b): DPPC; (c) and (d): PSM. (a) and (b) refer to the systems with Chol, while (c) and (d) are for those without Chol. Systems 01, 02, and 03 correspond to three distinct runs with different initial configurations.

Figure S3: Number density distributions of O${}^{\text{D}7}$ in DPPC and O${}^{\text{S}5}$ in PSM along the $z$-direction. For comparison, number density distributions of the water oxygen atoms (O${}^{\text{w}}$) are plotted, which are same as those in Fig. 2. Solid lines represent pure membrane systems, while dashed lines represent systems containing Chol. (a) and (b) refer to DPPC, and (c) and (d) refer to PSM.

Figure S4: Number density distributions of nitrogen (N⁺), and phosphorus (P) atoms along the $z$-direction. Solid lines represent pure membrane systems, while dashed lines represent systems containing Chol. (a) and (b) refer to DPPC, and (c) and (d) refer to PSM.

Figure S5: Snapshots at 323 K for (a) DPPC with Chol, (b) DPPC without Chol, (c) PSM with Chol, and (d) PSM without Chol.

Figure S6: 2D PMF, $W(r_{\mathrm{oo}},\beta)$, between water molecules in bulk water at 1 g/cm³ at (a) 303 K and (b) 323 K.

Figure S7: 2D PMF, $W(r_{\mathrm{oo}},\beta)$, between water oxygen ($\mathrm{O}^{\mathrm{w}}$) and each acceptor oxygen in DPPC with Chol at 323 K. Black points represent saddle points.

Figure S8: 2D PMF, $W(r_{\mathrm{oo}},\beta)$, between water oxygen ($\mathrm{O}^{\mathrm{w}}$) and each acceptor oxygen in DPPC without Chol at 303 K. Black points represent saddle points.

Figure S9: 2D PMF, $W(r_{\mathrm{oo}},\beta)$, between water oxygen ($\mathrm{O}^{\mathrm{w}}$) and each acceptor oxygen in DPPC without Chol at 323 K. Black points represent saddle points.

Figure S10: 2D PMF, $W(r_{\mathrm{oo}},\beta)$, between water oxygen ($\mathrm{O}^{\mathrm{w}}$) and each acceptor oxygen in PSM with Chol at 323 K. Black points represent saddle points.

Figure S11: 2D PMF, $W(r_{\mathrm{oo}},\beta)$, between water oxygen ($\mathrm{O}^{\mathrm{w}}$) and each acceptor oxygen in PSM without Chol at 303 K.

Figure S12: 2D PMF, $W(r_{\mathrm{oo}},\beta)$, between water oxygen ($\mathrm{O}^{\mathrm{w}}$) and each acceptor oxygen in PSM without Chol at 323 K. Black points represent saddle points.

Figure S13: Conditional probability $C_{1,j}(t)$, representing transition dynamics from region $1$ at the initial time $t=0$ to region 2 or remaining within the same region 1 during the time interval $t$ [(a) DPPC at 303 K, (b) DPPC at 323 K, (c) PSM at 303 K, and (d) PSM at 323 K].

Figure S14: H-bond time correlation function $P_{\mathrm{HB}}(t)$ with respect to acceptor oxygens in the DPPC systems. The solid line represents the result of fitting with the KWW function, $P_{\mathrm{HB}}(t)\simeq\exp[-(t/\tau_{\mathrm{KWW}})^{\beta_{\mathrm{KWW}}}]$.

Figure S15: H-bond time correlation function $P_{\mathrm{HB}}(t)$ with respect to acceptor oxygens in the PSM systems. The solid line represents the result of fitting with the KWW function, $P_{\mathrm{HB}}(t)\simeq\exp[-(t/\tau_{\mathrm{KWW}})^{\beta_{\mathrm{KWW}}}]$.

Table S2: $\tau_{\mathrm{KWW}}$ of acceptor oxygen atoms in the DPPC systems with and without Chol. The error is provided at the standard error and is not shown when it is smaller than $0.01$ ps.

	303 K [ps]		323 K [ps]
	pure	Chol 10 %	pure	Chol 10 %
$\mathrm{O}^{\mathrm{w}}$	$1.21$	$1.18$	$0.94$	$0.93$
$\mathrm{O}^{\mathrm{1}}$	$29.28\pm 0.02$	$41.7\pm 0.2$	${14.27}$	${15.07\pm 0.01}$
$\mathrm{O}^{\mathrm{2}}$	$24.97\pm 0.04$	$40.0\pm 0.3$	$8.81$	$9.35$
$\mathrm{O}^{\mathrm{3}}$	$22.26\pm 0.02$	$29.6\pm 0.1$	${10.32}$	${10.70}$
$\mathrm{O}^{\mathrm{D5}}$	$51.1\pm 0.1$	$88.7\pm 0.5$	${18.20\pm 0.02}$	${19.74\pm 0.02}$
$\mathrm{O}^{\mathrm{D7}}$	$40.2\pm 0.1$	$81.1\pm 0.4$	${15.01\pm 0.01}$	${16.74\pm 0.02}$
$\mathrm{O}^{\mathrm{Chol}}$	-	$110\pm 1$	-	$18.3\pm 0.1$

Table S3: $\beta_{\mathrm{KWW}}$ of acceptor oxygen atoms in the DPPC systems with and without Chol. The error is not shown since it is smaller than $0.01$.

	303 K		323 K
	pure	Chol 10 %	pure	Chol 10 %
$\mathrm{O}^{\mathrm{w}}$	$0.60$	$0.61$	$0.60$	$0.60$
$\mathrm{O}^{\mathrm{1}}$	$0.53$	$0.51$	$0.56$	$0.55$
$\mathrm{O}^{\mathrm{2}}$	$0.34$	$0.32$	$0.37$	$0.37$
$\mathrm{O}^{\mathrm{3}}$	$0.51$	$0.48$	$0.52$	$0.52$
$\mathrm{O}^{\mathrm{D5}}$	$0.42$	$0.39$	$0.45$	$0.44$
$\mathrm{O}^{\mathrm{D7}}$	$0.41$	$0.38$	$0.44$	$0.44$
$\mathrm{O}^{\mathrm{Chol}}$	-	$0.37$	-	$0.42$

Table S4: $\tau_{\mathrm{HB}}$ of acceptor oxygen atoms in the DPPC systems with and without Chol. The error is provided at the standard error and is not shown when it is smaller than $0.01$ ps.

	303 K [ps]		323 K [ps]
	pure	Chol 10 %	pure	Chol 10 %
$\mathrm{O}^{\mathrm{w}}$	$1.81$	$1.76$	$1.42$	$1.40$
$\mathrm{O}^{\mathrm{1}}$	$52.5\pm 0.5$	$81\pm 4$	$23.8\pm 0.2$	$25.4\pm 0.3$
$\mathrm{O}^{\mathrm{2}}$	$135\pm 3$	$280\pm 30$	$37.3\pm 0.5$	$40.9\pm 0.5$
$\mathrm{O}^{\mathrm{3}}$	$43.0\pm 0.7$	$64\pm 3$	$19.2\pm 0.2$	$20.1\pm 0.2$
$\mathrm{O}^{\mathrm{D5}}$	$151\pm 3$	$330\pm 30$	$46.1\pm 0.5$	$51.4\pm 0.7$
$\mathrm{O}^{\mathrm{D7}}$	$125\pm 3$	$310\pm 20$	$38.4\pm 0.4$	$44.4\pm 0.5$
$\mathrm{O}^{\mathrm{Chol}}$	-	$450\pm 60$	-	$53\pm 3$

Table S5: $\tau_{\mathrm{KWW}}$ of acceptor oxygen atoms in the PSM systems with and without Chol. The error is provided at the standard error and is not shown when it is smaller than $0.01$ ps.

	303 K [ps]		323 K [ps]
	pure	Chol 10 %	pure	Chol 10 %
$\mathrm{O}^{\mathrm{w}}$	$1.19$	$1.18$	$0.92$	$0.91$
$\mathrm{O}^{\mathrm{1}}$	$53.2\pm 0.1$	$65.8\pm 0.2$	$27.1\pm 0.1$	${28.41\pm 0.02}$
$\mathrm{O}^{\mathrm{3}}$	$39.9\pm 0.1$	$45.3\pm 0.1$	${17.68\pm 0.03}$	${19.14\pm 0.03}$
$\mathrm{O}^{\mathrm{S4}}$	$55.1\pm 0.3$	$62.9\pm 0.4$	$20.3\pm 0.1$	${21.08\pm 0.03}$
$\mathrm{O}^{\mathrm{S5}}$	$236\pm 1$	$400\pm 2$	$88.5\pm 0.4$	${109.4\pm 0.4}$
$\mathrm{O}^{\mathrm{Chol}}$	-	$670\pm 20$	-	$63.9\pm 0.7$

Table S6: $\beta_{\mathrm{KWW}}$ of acceptor oxygen atoms in the PSM systems with and without Chol. The error is not shown since it is smaller than $0.01$.

	303 K		323 K
	pure	Chol 10 %	pure	Chol 10 %
$\mathrm{O}^{\mathrm{w}}$	$0.61$	$0.61$	$0.60$	$0.60$
$\mathrm{O}^{\mathrm{1}}$	$0.48$	$0.46$	$0.50$	$0.50$
$\mathrm{O}^{\mathrm{3}}$	$0.44$	$0.42$	$0.46$	$0.45$
$\mathrm{O}^{\mathrm{S4}}$	$0.34$	$0.33$	$0.36$	$0.36$
$\mathrm{O}^{\mathrm{S5}}$	$0.39$	$0.37$	$0.41$	$0.41$
$\mathrm{O}^{\mathrm{Chol}}$	-	$0.31$	-	$0.36$

Table S7: $\tau_{\mathrm{HB}}$ of acceptor oxygen atoms in the PSM systems with and without Chol. The error is provided at the standard error and is not shown when it is smaller than $0.01$ ps.

	303 K [ps]		323 K [ps]
	pure	Chol 10 %	pure	Chol 10 %
$\mathrm{O}^{\mathrm{w}}$	$1.77$	$1.74$	$1.38$	$1.36$
$\mathrm{O}^{\mathrm{1}}$	$113\pm 4$	$153\pm 7$	$54\pm 2$	$57.3\pm 0.7$
$\mathrm{O}^{\mathrm{3}}$	$104\pm 2$	$129\pm 6$	$43\pm 1$	$47\pm 1$
$\mathrm{O}^{\mathrm{S4}}$	$320\pm 30$	$390\pm 30$	$95\pm 6$	$98\pm 2$
$\mathrm{O}^{\mathrm{S5}}$	$860\pm 80$	$1600\pm 100$	$280\pm 20$	$350\pm 20$
$\mathrm{O}^{\mathrm{Chol}}$	-	$6000\pm 1000$	-	$290\pm 30$

Figure S16: $\tau_{\mathrm{HB}}$ vs. $\Delta G/k_{\mathrm{B}}T$ for DPPC (a) and PSM (b) systems. The activation energy $\Delta G$ is determined by the free energy difference between the most stable and saddle points on the 2D PMF, $W(r_{\mathrm{OO}},\beta)$ (see Figs. 5 and 6, and Figs. S7-S12). The straight line represents the Arrhenius equation, $\tau_{\mathrm{HB}}=A\exp(\Delta G/k_{\mathrm{B}}T)$, where $A$ is determined by fitting. The values of $A$ are 0.183 ps and 0.385 ps for DPPC and PSM systems, respectively.