Structural properties of water confined by phospholipid membranes

Biological membranes are essential for the cell life and hydration water provides the driving force for their assembly and stability. Here we study the structural properties of water in a phospholipid membrane. We characterize local structures inspecting the intermediate range order (IRO) adopting a sensitive local order metric, recently proposed by Martelli et al., which measures and grades the degree of overlap of local environments with structures of perfect ice. Close to the membrane, water acquires high IRO and changes its dynamical properties, e.g., slowing down its translational and rotational degrees of freedom in a region that extends over $\simeq 1$ nm from the membrane interface. Surprisingly, we show that at a distance as far as $\simeq 2.5$ nm from the interface, although the bulk-like dynamics is recovered, water's IRO is still slightly higher than in bulk at the same thermodynamic conditions. Therefore, the water-membrane interface has a structural effect at ambient conditions that propagates further than the often-invoked $1$ nm-length scale, a results that should be taken carefully into account when analyzing experimental data of water confined by membranes and could help us understanding the role of water in biological systems.


Introduction
Biological membranes provide a limiting structure that separates the interior and exterior of cells and organelles. Being selectively permeable, membranes control the flow of substances in and out of the cell, which permits the regulation of the cell composition and communication between cells through signaling. Membranes are also involved in the capture and release of energy.
Biological membranes are composed of many biomolecules, including proteins, sugars, cholesterol, and phospholipids. Among these components, phospholipids provide structure to biological membranes. This is due to their spontaneous self-assembly, arising from *Special Topic: Water and Water Systems (Eds. F. Mallamace, R. Car, and Limei Xu). arXiv: 1703.07835 [cond-mat.soft] the hydrophobic effect and resulting in the formation of bilayers. For this reason, phospholipid membranes are used as models to investigate the fundamental properties of biological membranes, both experimentally and theoretically. There is a wide variety of phospholipids, but, in this study, we investigate dimyristoylphosphatidylcholine (DMPC), a phospholipid with a choline headgroup and a tailgroup formed of two myristoyl chains. Choline-based phospholipids are ubiquitous in cell membranes and are also used in drugs targeting liposomes [1]. The hydration water surrounding the membrane is fundamental for the proper functioning of biological membranes. Indeed, the presence of water strongly influences the stability, fluidity, and phase behavior of phospholipid membranes, determining their function. In addition, hydration water mediates the interactions of biological membranes with other biomolecules and with ions. Because of this importance, hydration water at phospholipid membranes has been extensively studied experimentally [2][3][4][5][6][7][8] and using computer simulations [9][10][11][12][13][14][15]. The results show that water is abundant in the interfacial region of bilayers (lipid headgroups) and can even penetrate into deeper regions of the membrane [10,15]. As a result, fully hydrated DMPC can take ≈ 30 hydration water molecules per lipid [15]. Furthermore, computer simulations show that water molecules form strong hydrogen bonds with the lipid phosphate groups, as well as with the carbonyl oxygens of the phosphatidylcholine lipids [9,10], slowing down both the orientational [5][6][7][13][14][15] and translational [4,15] dynamics of water. At low hydration, water is maintained inside the membrane, revealing its essential role in guaranteeing the structural stability of the membrane [15]. In highly hydrated systems, i.e., when lipid surfaces are in contact with many water molecules, The dynamical properties of bulk water, such as diffusivity, rotational dynamics, and density are recovered at relatively short distances from the lipid surface (see Calero et al. [15]). On the other hand, the structural characterization of water in such systems is lacking, mostly because structural properties are taken for granted once the density is known.
In this article, we characterize the local structure at the intermediate range order (IRO), i.e., the second neighbor shell, of water molecules confined between two membranes, each made of a lipid bilayer. We employ a sensitive local order metric (LOM) recently introduced by Martelli et al. [16]. The LOM characterizes the molecular order in the neighborhood of a site in a liquid, an amorphous compound, or a crystal. We compare the IRO of each water molecule with the IRO of perfect cubic Ic ice, and we grade the degree of deviation from the perfect Ic structure. 1) We measure the extent to which the lipid surfaces affect the structure of water at different distances from the surface, and we relate the structural changes to dynamical changes. We found that the membrane perturbs the IRO of water at a distance of at least ≃ 2.5 nm, while the bulk dynamical properties are recovered closer to the membrane (∼ 1.5 nm). This result should be considered when analyzing experimental and theoretical results of confined water in general. Furthermore, the results provide further insights into understanding the properties of water under such conditions [17,18]. Water molecules acquire increasing local order in regions close to the fluctuating lipid surfaces. This effect is small but not negligible, and we expect local order to be present (to some extent) in water in contact with hydrophobic surfaces, where rotational order in water molecules emerges [18]. In the present case, the increase of local order, as measured by the LOM, is ac-1) Similar results hold when examining the second neighbor shell on hexagonal ice as a reference structure.
companied by a slowing down in the translational and rotational degrees of freedom, in agreement with previous experimental and numerical observations [4][5][6][7][13][14][15]. Interestingly, on approaching the interface, the water molecules show a sudden drop in the number of translational degrees of freedom, reaching values comparable to water in deeply undercooled conditions. We also observed that, on approaching the lipid interface, the rotations around the water dipole are less affected than the rotations round the −→ OH vector, indicating that the interactions between P-groups and water hydrogen atoms are stronger -and more probable -than the interactions between N-headgroups and water oxygen atoms.
This article is organized as follows. In Section 2, we present the details of the LOM and the details of the numerical simulations. In Section 3, we show our findings. Our conclusions and final remarks are presented in Section 4.

The local order metric
Recently Martelli et al. [16] introduced a LOM that measures the degree of order in the neighborhood of an atomic or molecular site in a condensed medium.
The local environment of a water molecule i (i = 1, . . . , N ) in a configuration defines a local pattern formed by M neighboring sites. Here, we include only the oxygen atom second neighbors of the oxygen site i.
i.e., the average equilibrium ij distance in the pattern. The pattern-i centroid and the reference-i centroid are set to coincide, while the reference orientation is arbitrary. The local order metric, S(i), at site i is the maximum of the overlap function when the reference orientation is varied and the pattern indices are permuted, i.e., where θ, ϕ, and ψ are the Euler angles for a given orientation of the reference, i P are the permuted indices of the pattern sites corresponding to a permutation P, and σ ≡ d/4.4 is a parameter that controls the spread of the Gaussian functions. This value is chosen such that the tails of the Gaussian functions spread to half of the interatomic distance between oxygen atoms in the second neighbor shell in perfect cubic ice under ambient conditions. As a consequence of the point symmetry of the reference, the overlap function defined in Eq. (1) has multiple equivalent maxima. To compute S(i), it is sufficient to locate only one of these maxima. The point symmetry of the reference allows us to explore only a fraction (1/L) of the Euler angle domain Ω, which we call Ω/L. This is the irreducible domain of the Euler angles, L being the number of proper point symmetry operations of the reference. Inside Ω/L, we pick 15 orientations at random and with uniform probability, and we select permutations by assigning the nearest patternreference pairs. Using the best alignment out of the 15 orientations, we fix the permutation and optimize the orientation via conjugate gradients. The LOM is an intrinsic property of the local environment at variance with the straightforward overlap function that depends on the orientation of the reference and on the ordering of the sites in the pattern. The LOM satisfies the inequality 0 S(i) ≤ 1. The two limits correspond, respectively, to a completely disordered local pattern [S(i) → 0] and to an ordered local pattern matching perfectly the reference [S(i) → 1]. Consequently, this grades each local environment on an increasing scale of local order from zero to one.
We define order parameters based on S(i) as the average score, S, or site averaged LOM: and we refer to S as the score.

Simulations details
We performed molecular dynamics (MD) simulations of 7040 water molecules confined between two layers of 128 DMPC lipids.
2) We used periodic boundary conditions in all directions in such a way that our systems corresponds to water confined between the two different sides of the same lipid bilayer. The bilayer has an average width of 3 nm and the extension of the entire system along the direction perpendicular to the bilayer is 8.5 nm.
We used the simulation package NAMD 2.9 [19] at a temperature of 303 K and an average pressure of 1 atm. We set the simulation time step to 2 fs. We described the structure of the phospholipids and their mutual interactions by the recently parameterized force field CHARMM36 [20,21], which reproduces the area per lipid in excellent agreement with experimental data. The water model employed in our simulations, consis-2) We used a higher number of water molecules with respect to that considered in Ref. [15] for a similar system. tent with the parametrization of CHARMM36, is the modified TIP3P model [22,23]. The Van der Waals interactions were cut off at 12 Å with a smooth switching function starting at 10 Å. We computed the long ranged electrostatic forces using the particle mesh Ewald method [24], with a grid space of about 1 Å. After energy minimization, we equilibrated the hydrated phospholipid bilayers for 10 ns, followed by a production run of 2 ns in the N pT ensemble at 1 atm. In the simulations, we controlled the temperature using a Langevin thermostat [25] with a damping coefficient of 0.1 ps −1 , and we controlled the pressure using a Nosé-Hoover Langevin barostat [26] with a piston oscillation time of 200 fs and a damping time of 100 fs.

Structure
We characterized the local structures by using, as a reference, a Cuboctahedron,C ( Fig. 1), belonging to the class of Archimedean solids enriched with edgetransitivity [27] and describing the oxygen lattice in the second neighbor shell as perfect cubic Ic ice. Analogous results can be drawn using an anticuboctahedron (C), which describes the oxygen lattice in the second shell of neighbors as perfect hexagonal Ih ice. In Fig. 1, the six-folded ring on the σ-plane is a plane of inversion. In the case of hexagonal ice, one of the two three-folded rings above or below the σ-plane is rotated by 60 • .
We divided the entire system along the direction perpendicular to the bilayer in 10 bins and calculated, for the water molecules within each bin centered at a distance z from the center of the bilayer, the average score, SC [ Fig. 2(a)]. We found that SC is approximately constant in the central region within 2.8 nm < z < 5.8 nm [vertical dashed orange lines in Fig. 2(a)] of the lipid layers. The score slightly increased moving from the center toward the lipid interfaces at z 2.8 nm and z 5.8 nm. The sudden drops in SC at z ≃ 1.5 nm and z ≃ 7 nm correspond to the average position of the fluctuating water-membrane interface (vertical dashed-dotted blue lines). The behavior of SC suggests that the interface has a significant effect on the structure of water at distances of < 1.5 nm from the interface, while, at larger distances, water might be bulk-like, at least concerning some properties.
To better understand how, and to what extent, confined water becomes bulk-like, we calculated the number of water molecules as function of the distance z (or z ′ ≡ 8.5−z) from the center of the membrane [ Fig. 3(a)]. We found that, on average, water penetrates the membrane by up to z ≃ 1.5 nm (or z ′ ≃ 1.5 nm, i.e., z ≃ 7 nm) and populates the confined region between the membranes with a density profile that saturates at a distance ≃ 1.5 nm from the membrane at z ≃ 3 nm (or z ′ ≃ 3 nm, i.e., z = 5.5 nm).
These results are consistent with the results of a previous study [15], where Calero et al. introduced the intrinsic distance (ξ) from the membrane by performing a twodimensional Voronoi tessellation of the average plane of the membrane (the xy-plane) using the phosphorous and nitrogen atoms of the phospholipid heads as centers of the Voronoi cells. The introduction of ξ ≡ z − z Voronoy , where z Voronoy is the z coordinate of the center of the Voronoi cell, allows one to characterize the water density profile while filtering out the noise induced by the fluctuations of the water-membrane interface. In particular, the use of ξ emphasizes the penetration of water within the membrane, the layering of water near the membrane, and the existence of a central region with an approximately constant density profile ≃ 0.5 nm away from the interface when the hydration level is comparable to that considered here [15]. Calero et al. showed that the density of water is the same as that of bulk water under the same thermodynamic conditions at distances 0.5 nm from the membrane. This observation, therefore, suggests that water in the central region has almost bulk-like properties.
A further comparison with the results of Calero et al. [15] shows that the increment of SC towards the lipid interface resembles the increase of water density, ρ, in the same direction. This concomitant increment of ρ and SC might seem counterintuitive, but we rationalize that it is a consequence of the isotropic reduction of the distances between all second shell molecules (not shown). Furthermore, the increase in SC (i.e., the IRO) is also related to the strengthening of the interactions of water with the lipid membrane, as shown by our analysis of the hydration water dynamics in the following subsection.
Based on the previous work [15] and from the current analysis [ Fig. 3(a)], we focused our structural analysis in a central region at a distance of 2.5 nm from the membrane, i.e., 4.0 nm < z < 4.4 nm. We compared the IRO in this region with the IRO of bulk water as ob- tained from simulations of 6000 water molecules under the same thermodynamic conditions and with the same interaction potentials. In particular, we computed SC only for those water molecules with all 12 second neighbors within the region with a constant density profile (2.5 nm < z < 5.5 nm). We found that SC for the central region follows a P (SC) distribution that differs slightly, but neatly, from the bulk case [ Fig. 3(b)]. In particular, bulk water is slightly less ordered than confined water. This result suggests that the lipid membrane affects the IRO of water at distances as far as ≃ 2.5 nm from the interface, corresponding to approximately seven water diameters. We infer that a possible source of such longrange effects is the fluctuations of the membrane, which create dynamical density heterogeneities.

Dynamics
The interface affects the density and structure of the water, and this is accompanied by changes in the dynamical behavior of confined water. For example, Calero et al. [15] found that the coefficient of diffusion parallel to the membrane and the rotational relaxation time of water confined between lipid membranes approach within 25% of the bulk values when water moves further than 0.5 nm from the interface.
Here, to study the dynamics of the membrane-confined water, we collected data over 1 ns trajectories and calculated the standard displacement, defined as the distance traveled by water molecule i normalized to the ith O-O distance, BU d i , averaged over all the molecules i. We performed the calculation bin by bin [ Fig. 2(b)].
Our data show that the water molecules in the central region between the lipid layers travel up to ≃ 10 BU, while they significantly slow down approaching the interface, where SC increases [ Fig. 2(b)]. For the water within 1 nm of the interface, i.e., at z 2.5 nm and z 6 nm, the displacement drops by 60% to ≃ 4 BU. For water penetrating the membrane, i.e., z 1 nm, the translational motion reduces by a further 50% to ≃ 2 BU. A standard displacement of < 1 BU would correspond to water molecules rattling in the cage formed by their nearest neighbors. This case would represent a liquid in which the translational degrees of freedom are frozen. Therefore, our data show that the interface has a significant effect on the water dynamics, at least within a distance of 1 nm from the membrane.
To understand the magnitude of this effect, we made a comparison with the bulk case. The definition of standard displacement depends on the total time of the simulation; therefore, we made comparisons with benchmark cases calculated for the same simulation time. As a reference case, we used 1 ns trajectories for bulk TIP4P/2005water. Although we simulate here TIP3P-water, whose properties are fairly different from those of TIP4P/2005-water, a qualitative comparison of the results with the two models is still possible. In particular, we observed that water close to the lipid interfaces reaches values similar to those found for supercooled TIP4P/2005-water after 1 ns simulation time [16].
This suggests that the increment in the IRO discussed in the previous subsection can be ascribed to the significant slowing down of the translational degrees of freedom. Indeed, a decrease in the translational diffusion can be interpreted as a reduction in the magnitude of the thermal noise. Hence, the configurational entropy contribution to the free energy minimization becomes less important, while the potential energy term, arising from the electrostatic interactions, acquires more weight, leading to a more structured configuration.
Further insights into the effects of the dynamical properties on structural properties can be achieved by calculating the correlation functions: where A is either (i) the dipole vector µ or (ii) the −→ OH vector. In particular, Calero et al. showed [15] that C µ (t) yields information concerning the dynamical processes involving molecular rotations.
On the other hand, DMPC lipids contain N-and Pheadgroups, where the N atoms interact mostly with water oxygen atoms and the P atoms interact mostly with water hydrogen atoms, denoted HO. It is reasonable to assume that the N-O and the P-HO interactions have different strengths. This assumption can be directly tested by calculating the time correlation function C(t) ≡ ⟨δ(t) · δ(0)⟩, where δ is the N-O vector or (black) and for the P-HO vector (red) at the water-membrane interface. Exponential fits of the calculations yield estimates of the correlation time that are larger for the P-HO vector. We calculated C(t) by averaging over distances up to the first minimum in the N-O (black) and the P-HO (red) radial distribution functions (inset). the P-HO vector (Fig. 4). We found that the P-HO vector has a longer lifetime compared to that of the N-O vector, indicating that the interactions between P and water hydrogen atoms are stronger than the interactions between N and O. This conclusion is also consistent with the observation that the P-HO radial distribution function (r.d.f.) has its first peak at shorter distance than the N-O r.d.f. (inset in Fig. 4).
Based upon the observation that the N-O and the P-HO interactions have different strengths, we assume that rotations around µ are different from those around −→ OH. In bulk water under these thermodynamic conditions, the two rotations are, instead, indistinguishable. Here, we expect them to play a different role in the structure of water at the interface.
We, therefore, calculated C µ (t) and C−→ OH (t) for water molecules belonging to bins centered at different z (not shown) and fitted them with a double exponential function.
3) From the fits, we extracted two characteristic times for each z, one for each exponential function: τ 1 (z), associated with fast relaxation modes, and τ 2 (z), associated with the slow relaxation modes (Fig. 5). We interpreted the two families of relaxation processes as resulting from the dynamical heterogeneities induced by the interactions of water molecules with the lipid layers.
We found that the rotation of water slows as molecules approach the membrane. In particular, the distances at which the change is appreciable are the same as those for for the OH vector (red squares) as a function of z. Vertical (orange) lines indicate the region where the variations in SC become appreciable and are at the same z as in Fig. 2. (b) Same as panel (a) but for the slow relaxation times τ2.
3) Calero et al. [15] correctly observed that it is formally more correct to integrate the correlation function to obtain the rotational time τ rot ≡ ∫ ∞ 0 C A (t)dt. Although the analysis of τ rot yields qualitatively similar conclusions to that presented here, we have adopted the common biexponential fitting because it intuitively reveals the effects of the electrostatic interactions on the slow relaxation, which play an important role in our analysis. the translations, at z 2.8 nm and z 5.8 nm.
We observed that, On approaching the lipid membrane, both τ 1 and τ 2 increase. Interestingly, the relaxation times for −→ OH increase at a pace higher than the relaxation times for µ. This observation is consistent with our finding that the P-HO interaction is stronger than the N-O interaction, and this can be rationalized by observing that the lipids have different (delocalized) charges on the N-heads and on the P-functional groups and that these charges affect the rotation of water around the two vectors in different way.
We found that the rotational relaxation times of water penetrating the membrane, at z < 1.5 nm, increased significantly, as it would be expected for degrees of freedom near their freezing point. In addition, in this case, the slowing down is stronger near the P-groups than near the N-heads.
The slowing down of the rotational dynamics decreases as water moves from the interface toward the central region, i.e., 2.8 nm z 5.8 nm. At these distances, we found that the relaxation times are almost indistinguishable, and τ µ 1 ≃ τ µ 2 ≃ τ OH 1 ≃ τ OH 2 , as one would expect in bulk water.
Therefore, the water recovers the bulk dynamic behavior ≃ 1.3 nm away from the membrane interface, corresponding, approximately, to four water diameters. However, this result contrasts with the structural analysis described in the previous section, where we found that the IRO structure is affected as far as ≃ 2.5 nm away from the interface. We have interpreted this difference as a consequence of the fact that the structural parameter SC has a higher sensitivity, with respect to the correlation functions, to small effects arising from the interface. In particular, because we have adopted a definition for the second water hydration shell, the effect on SC is observed up to three water diameters further than those for the correlation times, reaching a total of seven water diameters, consistent with the conclusion of the previous subsection.

Conclusions
We have employed the sensitive order metric S to characterize the structural properties of water confined between phospholipid membranes. Commonly inspected quantities, such as the O-O radial distribution function, molecular displacement, and rotational correlation functions, do not capture differences between bulk water and the farthest water layer from the lipid surface. Instead, we show that the high sensitivity of S allows the detection of small, but not negligible, structural differences in the IRO. The lipid membrane, therefore, perturbs the structure at the IRO of water as far as ≃ 2.5 nm from the lipid surface. Notably, this value is limited by the finite size of our sample. Therefore, larger simulations are required to better estimate the extent of such long-range effects.
We found that the IRO of water increases as the membrane surface is approached surfaces, and this is concurrent with the slowdown of the translational and rotational water dynamics. Our observations indicate that water translational and rotational correlation times increase significantly as water approaches the membrane. The translational times reach values comparable to that of water at supercooled conditions. On the other hand, we show that the membrane acts unevenly on the rotational degrees of freedom, slowing down the rotations related to the −→ OH vectors more than the rotations related to the water dipole vector. Our calculations reveal that this is due to the stronger interaction between Pgroups and water hydrogen atoms with respect to the interactions between N-heads and oxygen atoms.
Such effects have physical implications for the structure and local density of water, and molecules interacting with P-groups are slightly more ordered than molecules interacting with N-heads. Moreover, water molecules interacting with P-groups have closer second neighbors than those interacting with N-heads. Although only qualitative, these observations are in agreement with our considerations linking the slowing down of the translational degrees of freedom with the increment in SC. Further analysis using the intrinsic distance, ξ, from the membrane [15] and in-depth statistical analysis will be necessary to make these observations quantitative.
In conclusion, our results suggest that the effect of the phospholipid membrane on the structural properties of water extend as far as ≃ 2.5 nm from the membrane, a larger distance than previously calculated using less sensitive observables. The bulk density and dynamical properties are recovered at much shorter distances from the membrane, i.e., ∼ 1 nm, showing that the definition of a bulk-like region depends intrinsically on the property used to compare with the bulk. This effect should, therefore, be considered when studying properties of water at interfaces. Accordingly to our findings, the order metric S is an ideal tool to inspect such structural changes.