Quantifying dynamical heterogeneity length scales of interface water across model membrane phase transitions

Lipid bilayers form a continuous semipermeable barrier surrounding cells and compartmentalize the intracellular space. They are essential for controlling the dynamics, arrangement, and function of membrane proteins¹ at the fluid (or liquid-crystalline, L_α) phase at physiological temperature. The bilayers at the fluid phase comprise heterogeneous functional microdomains that aid in activities including cell signaling, cell adhesion, and membrane trafficking. Thus, the lateral heterogeneity is accepted as a criterion for the function of a complex biological membrane, specifically in the context of a raft² relevant for signaling cascades.³ Quantifying the heterogeneity in a bio-membrane using the current biochemical or biophysical tools is an incredibly challenging task due to the presence of unusual liquid-like domains with microseconds to seconds time scale and nanometers to micrometers length scales having both short- and long-range orders. Molecular-level characterization of such dynamical domains through experiments is still fragmentary due to the lack of relevant atomistic trajectories and the unknown spatiotemporal correlations.

At the physiological temperature interface, water molecules (IW) near a single component fluid membrane exhibit signatures of dynamical heterogeneity⁴ and can probe regional lipid dynamics that is difficult to access otherwise.⁵ This is the stepping stone of the current study that aims to measure the spatiotemporal scales of different phases of functionally relevant bilayers through the IW dynamics. Using sum-frequency generation, the spectroscopy heterogeneity of IW at the water-charged lipid interface has been identified.⁶ The fluid or disordered membrane undergoes a phase transition to a ripple (P_β) phase, which quickly changes to a gel (L_β) phase⁷ upon cooling. The membrane remains in the gel phase upon supercooling. Similar disorder-to-order transitions of the membrane occur from dehydration-induced drastic slow-down in structural relaxation originating from the dynamical heterogeneity of the interface water molecules (IW).⁸ Understanding cellular membranes at supercooled temperatures are extremely crucial for cryo-preservation methods, which have a wide range of applications in food science, pharmacology, clinical medicine, and the preservation of human embryos for in vitro fertilization.⁹ Despite such innumerable applications, investigations on the lipid dynamics at the gel phase are rarely found.¹⁰ Using quasielastic neutron scattering (QENS), a higher contribution of rotational water dynamics is found near multilamellar 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) bilayers compared with the translational one at supercooled temperature.¹¹ Although it is widely accepted that IW aids in different activities of cellular membranes, very few works attempt to establish the correlations of hydration dynamics to membrane functions.¹² Using coarse-grained (CG) and all-atom (AA) molecular simulations and different experimental techniques, the phase behavior of model membranes has been investigated extensively.^13–17 However, none of the work mentioned above focuses on the underlying role of dynamical heterogeneity across phase transitions. The thermodynamics and dynamics of the IW exhibit sharp cross-overs across the bilayer melting transition, which demonstrates their correlations.¹⁸ Thus, the effect of membrane phase transitions on the dynamical heterogeneity of the IW is fundamentally important and long overdue. Although many theoretical models¹⁹ and numerical methods are developed to study glass transition phenomena, their applicabilities on the membrane phase transitions and water dynamics are not understood. The current research finds out the relevance of different theories²⁰ and recent numerical approaches²¹ of glass dynamics on model membranes across phase transitions. The objective is manyfold: (a) to gain a fundamental understanding of the dynamical heterogeneity of the IW and the lipid membranes near phase transitions, (b) to establish relations of the spatiotemporal scales of the IW to that of the membranes, and (c) to build up an alternate method through IW dynamics to probe membrane domains under physiological and stressed conditions. Furthermore, this strategy will help in comprehending the domain length scale of a raft, which is rarely found in theory and extremely difficult to gain from an experiment. Few theoretical analyses and simulations are reported^22–24 to measure the dynamical length of persistent domains with distributions of the most mobile and immobile water and its correlation with the structure as a function of distance from the interface.

Therefore, 11.55 µs long all-atom molecular dynamics (MD) simulations of DMPC bilayers are performed at various temperatures and two different system sizes (6 × 6 and 18 × 18 nm²) to cover the three primary phases. Simulation details are mentioned in Sec. I of the supplementary material, and the run lengths of all bilayers are mentioned in Table S1 of the supplementary material. CHARMM36 force fields are used for the 128 DMPC lipids solvated in 5743 TIP4P/2005²⁵ (6 × 6 nm² box size) water molecules. This particular pairing of the lipid force fields and the water model accurately reproduces the fluid to ripple to gel phase transitions.²⁶ The area per head groups (a_h) of membranes from our simulations at all phases match well with the available literature^27,14,28 (Sec. II and Table S2 of the supplementary material). Snapshots of the equilibrated bilayers at three phases are shown in Figs. 1(a)–1(c), and the bilayers at the remaining temperatures are shown in Fig. S1 of the supplementary material (considered for the calculations of structural relaxations). The a_h drops abruptly from 308 to 292 and from 284 to 268 K, demonstrating the fluid (L_α) to the ripple (P_β) to the gel (L_β) phase transition temperatures (Fig. S2 of the supplementary material). The ripple phase is identified from the existence of periodic interdigitated and non-interdigitated regions of lipid tails.²⁹ The homogeneous thickness and interdigitation of the fluid phase (324–303 K) [Figs. 1(d) and 1(g)] become periodically heterogeneous for the ripple phase (292–284 K) [(e) and (h)] that disappears for the gel phase (273–218 K). The region with a smaller thickness is correlated with higher interdigitation at all phases. The thickness and interdigitation of the remaining temperatures are shown in Fig. S4 in the supplementary material (interdigitation). At the ripple phase, there is a coexistence of a highly interdigitated ordered domain followed by a fully expanded disordered domain (shown in S4 in the supplementary material). The ripple wavelength of ∼4.2 nm from Figs. 1(e) and 1(h) matches well with that obtained from the spectral intensity of thickness and interdigitation fluctuations discussed in Sec. II of the supplementary material.

The homogeneous thickness and interdigitation of the fluid phase (324–303 K) [Figs. 1(d) and 1(g)] become periodically heterogeneous for the ripple phase (292–284 K) [(e) and (h)] that disappears for the gel phase (273–218 K). The region with a smaller thickness is correlated with the higher interdigitation at all phases. The thickness and interdigitation of the remaining temperatures are shown in Fig. S4 in the supplementary material (interdigitation). At the ripple phase, there is a coexistence of a highly interdigitated ordered domain followed by a fully expanded disordered domain (shown in S4 in the supplementary material). The ripple wavelength of ∼4.2 nm from Figs. 1(e) and 1(h) matches well with that obtained from the spectral intensity of thickness and interdigitation fluctuations discussed in Sec. II of the supplementary material. As the largest survival time of the water molecule near membranes at the highest temperature is found to be ∼100 ps,³⁰ it is considered as the minimum residence time of the IW at all temperatures. Thus, if a water molecule is present within ±0.35 nm of the lipid P head atoms continuously for 100 ps, it is classified as an IW³¹ (Sec. III of the supplementary material). To find out the structural relaxation, self-intermediate scattering functions (SISF) of the bulk water (BW), IW, and lipid heads are calculated (Sec. IV of the supplementary material). The initial decay of the SISF [Figs. 2(a)–2(c)] is due to the ballistic or cage motion followed by the emergence of Boson peaks at lower temperatures (228–218 K), as observed commonly in supercooled liquids. The time at which F_s(q, t) = 1/e is referred to as the α relaxation time-scale, τ_α. Note, τ_α extracted by this method matches well (data not shown) with the values of τ_α obtained from the fitting of the modified Kohlrausch-Williams-Watts (KWW) equation^5,32 with an additional time scale of τ_l. With lowering the temperature, τ_α starts growing for the BW, IW, and P heads. The growths are significant for the IW and P, once the temperature goes below 273 K and the membrane is at the gel phase. The growing τ_α of the IW and P heads are strongly correlated at all three phases of the membranes even at super-cooled temperatures (Fig. S5 in the supplementary material). Block analysis has been employed on the membrane interface with an 18 × 18 nm² surface area (Table S1 of the supplementary material). The non-Gaussian van Hove distributions, G_s(X, t), of the IW on the membrane surface become more and more Gaussian as the L_B increases [see Figs. 3(a), 3(d), and 3(g) for the fluid, ripple, and gel phases, respectively, and Fig. S6 of the supplementary material for remaining temperatures]. Binder cumulant, B(L_B, T) = 0, for higher temperatures and larger blocks as the underlying distribution is Gaussian. With the decreasing block size, the values of the B(L_B, T) increase at all temperatures due to larger deviations from Gaussianity seen in Figs. 3(b), 3(e), and 3(h) for three phases of the membrane. As the binder cumulant is a scaling function of the underlying correlation length, it is assumed to follow the relation of B(L_B, T) = f[L_B/ξ(T)]. Thus, B(L_B, T) at different temperatures are subjected to a data collapse and plotted in Figs. 3(c), 3(f), and 3(i) for the fluid, ripple, and gel phases, respectively. The insets of Fig. 3 depict the temperature dependence of the scaling length scale, ξ, known as the heterogeneity length scale. For the fluid L_α phase, ξ increases with a decrease in T. The dependence disappears for the ripple P_β phase, unlike the glassy liquids.

The values of ξ at the ripple phase are close to half of the ripple wavelengths as obtained from the spectra of thickness or interdigitation at all temperatures (Fig. S7 of the supplementary material). This signifies that the periodic patch or domain of the ripple of the membrane can be captured by the heterogeneity length scale of the IW that does not vary with temperature within the same phase. The T dependency of the ξ is very weak at the gel phase. The heterogeneity length scale follows a power law behavior with the amplitude of the first peak of the radial distribution functions, g₁³⁰ of the IW at each of the three phases of membranes [Figs. 4(a)–4(c)], suggesting that the correlation length scale is structure dominated. Figure 4(d) shows that the change in ξ is minor within a given phase but the length scale significantly grows once there is a phase transition. However, within a given phase, there is a drastic slow-down in the τ_α with decreasing temperature without a significant change in ξ [Fig. 4(e)]. On the contrary, there are sharp changes in ξ, across the membrane phase transitions. Thus, the heterogeneity time scales are slowed down mainly by supercooling, which is not accompanied by a similar growth in the heterogeneity length scales. The large growth in the length scales is governed by the membrane phase transitions. The relaxation time scale, τ_α, of the IW obeys the relation, τ_α ∼ e^aξ/T, a = constant, solely near the gel phase [Fig. 4(f)], which is not valid at the fluid or the ripple phase (Fig. S8 of the supplementary material).

To correlate the dynamics of the IW with the structure of the bilayers and to visualize dynamical heterogeneity, IW is classified into two groups, caged and mobile, at three phases of the bilayers using two-dimensional (2D) van Hove correlation functions. Similar to the iso-configurational analysis (ISOCA) that shows the most mobile and the most immobile water calculated from local diffusivity,²² the caged and mobile IW are calculated based upon two cut-off distances, r₁ and r₂,^8,24 obtained [Figs. 5(a)–5(c) for three phases] from the deviations of the 2D simulated van Hove correlation functions, $Gssimulation(r,t)$⁠, of the IW from that of the theoretical van Hove correlation function, G^theory(r, t) as mentioned in the supplementary material. If the displacement of an IW r < r₁, the IW is considered caged. If r > r₂, it is a mobile IW. The caged and mobile IW on the bilayer surface are color coded and shown in Figs. 5(d)–5(f) for the fluid, ripple, and gel phases. The number of caged IW or the most immobile IW increases from the fluid to the ripple to the gel phase. A correlation length can be estimated from the probability of two molecules belonging to the same mobile cluster separated by the distance r as developed in Ref. 24 and compared with ξ extracted from the block analysis. Using this approach, the correlation between the number of caged IW and the change in the bilayer structure from the fluid to the ripple to the gel phase can be estimated. Such a correlation between the structure and dynamics at the atomic scale can be useful for both modeling and experimental investigations in the future.

In conclusion, the study uses 11.55 µs all-atom molecular dynamics simulations to correlate and quantify the spatiotemporal heterogeneities of the IW to the DMPC lipid membranes upon fluid to ripple to gel phase transitions. Although there is an α-relaxation time lag between the IW and the lipids, the lipid heads and the IW are strongly correlated with each other at all three phases with a strong enhancement in dynamical heterogeneity upon supercooling. The length scales of spatially heterogeneous dynamics of the IW are quantified using a block analysis technique and then correlated with the lateral organizations of the membranes across phase transitions. One-dimensional van Hove correlation functions of the IW show a shift from the non-Gaussianity to Gaussianity once the block size increases at each phase of the membrane. The deviations from the Gaussianity are estimated by Binder cumulant that increases with lower block size at all temperatures. Remarkably, the size of the rippling of the membrane can be obtained from the heterogeneity length scale of the IW. The length scale does not change significantly upon supercooling at the gel phase and follows an activated dynamical scaling with the relaxation time scale. This observation is in a similar line to what is predicted by the random-field Ising model.³³ The drastic growth in the ξ is dominated by the membrane phase transitions and the drastic slow-down in the relaxation time scale is dictated by the supercooling. Thus, the dynamic cross-over in the heterogeneity length scale of the IW will provide a new avenue to identify the micro-domain that can be critical for membrane functionalities. Since probing nano-domains in biological membranes are extremely challenging and no current biophysical tool is available yet to this end,² our method opens up a new direction to characterize functionally relevant nano-domains of membranes through the heterogeneity length scale of the IW. As per our knowledge, for the first time, a systematic analysis is performed to estimate the spatiotemporal correlations of the hydration water to the model lipid membranes across phases relevant to physiological and low-temperature stressed conditions. The findings can help us to better understand bio-protection processes under extreme conditions and will aid future research on domain-associated transport and signaling in biomembranes at low temperatures.

See the supplementary material for the simulation details, bilayer properties, identification of the interface water (IW), and dynamics including tables and figures.

A.D. acknowledges the financial support of the Grant No. SERB CRG/2019/000106. S.M. was thankful to Dr. Indrajit Tah for useful discussions and to Professor Jeffery B. Klauda for providing the pdb of the equilibrated DMPC bilayer (with 72 lipids) at 273 K.