Properties of water hydrating the galactolipid and phospholipid bilayers: a molecular dynamics simulation study*

Molecular dynamics simulations of 1,2-di-O-acyl-3-O-β-D- galactopyranosyl- sn -glycerol (MGDG) and 1,2-dioleoyl- sn -glycero-3-phosphatidylcholine (DOPC) bilayers were carried out to compare the effect of the lipid head group’s chemical structure on the dynamics and orientational order of the water molecules hydrating the bilayer. The effect of the bilayers on the diffusion of water is strong for the neighbouring water molecules i.e., those located not further than 4 Å from any bilayer atom. This is because the neighbouring water molecules are predomi- nantly hydrogen bonded to the lipid oxygen atoms and their mobility is limited to a confined spatial volume. The choline group of DOPC and the galactose group of MGDG affect water diffusion less than the polar groups located deeper in the bilayer interface, and similarly. The latter is an unexpected result since interactions of water with these groups have a vastly different origin. The least affected by the bilayer lipids is the lateral dif- fusion of unbound water in the bilayer plane ( x,y -plane) — it is because the diffusion is not confined by the pe- riodic boundary conditions, whereas that perpendicular to the plane is. Interactions of water molecules with lipid groups also enforce certain orientations of water dipole moments. The profile of an average water orientation along the bilayer normal for the MGDG bilayer differs from that for the DOPC bilayer. In the DOPC bilayer, the ordering effect of the lipid head groups extends further into the water phase than in the MGDG bilayer, whereas inside the bilayer/water interface, ordering of the water dipoles in the MGDG bilayer is higher. It is possible that differences in the profiles of an average water orientation across the bilayer in the DOPC and MGDG bilayers are responsible for differences in the lateral pressure profiles of these bilayers.


INTRODUCTION
The classical dogma that every protein function determined by its structure is encoded in the respective structural gene is, at least, inaccurate. In fact, to perform their biological functions, proteins often require ions but most of all water. Neither the presence nor the structure of water is encoded in genome but water is an indispensable element of any living organism in which it fulfils important biological functions. It is also an essential structural component of biological membranes. On the molecular level, water molecules readily interact with biomolecules and biosurfaces (Fogarty & Laage, 2014). These interactions modify the properties of water molecules in close vicinity of the interacting partner, e.g., (Gawrisch et al., 1992;Marrink et al., 1993;Cheng et al., 2003). It has been shown that the dynamics of water near surfaces of membranes (König et al., 1994;Balasubramanian et al., 2002), proteins (Bagchi, 2005), nucleic acids (Michalarias et al., 2005) and sugars (Köper et al., 2005) is hindered. Our previous MD simulation study on the dynamics of water hydrating the 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine (POPC) bilayer demonstrated that both translational and rotational motion of water molecules was affected by the bilayer (Rog et al., 2002). That study was later extended on other phospholipid head groups, phosphatidylethanolamine (PE), and phosphatidylglycerol (PG), and showed that despite different head group structures and hydrogen bonding capacities, the effect of PC, PE, and PG on the water dynamics was similar (Murzyn et al., 2006). In the present study, the properties of the water molecules hydrating the bilayer composed of 1,2-di-O-acyl-3-O-β-dgalactopyranosyl-sn-glycerol (MGDG, Fig. 1a) are compared to those hydrating the bilayer composed of 1,2-di- Only the polar hydrogen atoms (empty circles) of MGDG are displayed; the oxygen (O), nitrogen (N), and phosphorus (P) atoms are dark and the carbon atoms are light grey, respectively; the chemical symbol for carbon atoms, C, is omitted. Vol. 62, No 3/2015 475-481 http://dx.doi.org/10.18388/abp.2015_1077 oleoyl-sn-glycero-3-phosphatidylcholine (DOPC, Fig.  1b). MGDG is the main lipid component of thylakoid membranes of higher plants and algae, and plays an active role in the process of photosynthesis however; molecular level studies on model MGDG membranes are not numerous so the effect of MGDG on the hydrating water has not been described so far. In contrast, DOPC bilayers have been extensively studied but the effect of DOPC on the hydrating water has not been analysed in details either. Comparison of the influence of the two lipids with structurally diverse head groups -galactose and phosphatidylcholine -on the dynamics of water is essential to fully explain our results on structural organisation of the bilayer/water interface and elastic properties of the galactolipid and PC bilayers obtained by us previously (Baczynski et al., 2015) 1 .

METHODS
The MGDG and DOPC bilayers were previously equilibrated in 440-ns molecular dynamics simulations (Baczynski et al., 2015). In this study, the simulations were extended for additional 10 ns using Gromacs 5.0.4 package (Pronk et al., 2013) to investigate in detail the properties of the bilayer water. Each bilayer consisted of 8 × 8 × 2 (128) lipid molecules. The MGDG bilayer was hydrated with 3840 and the DOPC bilayer with 6747 water molecules (Fig. 2). For lipids, the all-atom optimized potentials for liquid simulations (OPLS-AA) force field (Jorgensen et al., 1996;Kaminski et al., 2001) was used; for water the transferable intermolecular potential three point model (TIP3P) was used (Jorgensen et al., 1983). A 5366-molecule water box simulated for 20 ns constituted the reference system. The van der Waals interactions were cut-off at 1.0 nm. The long-range electrostatic interactions were evaluated using the particle-mesh Ewald summation method with the β-spline interpolation order of 5, and a direct sum tolerance of 10 −6 (Essmann et al., 1995). For the real space, a cut-off of 1.0 nm, threedimensional periodic boundary conditions, and the usual minimum image convention, were used (Essmann et al., 1995). The linear constraint solver (LINCS) algorithm (Hess et al., 1997) was used to preserve the length of any covalent bond with a hydrogen atom, and the time step was set to 2 fs. The list of non-bonded pairs was updated every 5 steps.
The MD simulations were carried out in the NPT ensemble, under a constant pressure of 1 atm, and at constant temperature of 295 K (22 °C), which is above the main phase transition temperature for the pure DOPC bilayer of 253 K (-20°C) (Kučerka et al., 2006); at 295 K temperature, the MGDG bilayer is in the liquid-crystalline phase (Baczynski et al., 2015). Temperatures of the solute and solvent were controlled independently by the Nosé-Hoover method (Hoover, 1985), with the relaxation time of 0.6 ps. The pressure was controlled anisotropically by the Parrinello-Rahman method (Parrinello & Rahman, 1981), with the relaxation time of 1.0 ps. For analyses, 10-ns trajectories (the last 10 ns for the reference water system) sampled every 50 fs, were used. To analyse water diffusion and orientation of bound water, each trajectory was divided into 200 50-ps fragments; the results presented below are averaged over the fragments.

RESULTS
To study the effect of the bilayer on the dynamics of the hydrating water, the water molecules were classified into seven groups, following the criteria used in Refs (Rog et al., 2002;Murzyn et al., 2006). The first group consisted of water molecules that were not further than 4 Å from any bilayer atom (neighbouring water), the second group consisted of water molecules that remained within a layer between 4 and 12 Å from any bilayer atom (intermediate water), and the third group contained water molecules that were not closer than 7 Å to any bilayer atom (far water). Next four groups consisted, respectively, of water molecules hydrogen bonded (H-bonded) to the phosphate oxygen atoms of DOPC (Op water); the carbonyl oxygen atoms of DOPC and MGDG (Oc water), and the galactose oxygen atoms of MGDG (Og water) as well as those clathrating the choline group of DOPC (Nch water) -H-bonding and clathrating criteria are given in (Pasenkiewicz-Gierula et al., 1999;Murzyn et al., 2001). A water molecule was qualified to one of the groups if it fulfilled the given criteria for 100% 1 The manuscript was submitted to Colloids Surf. B: Biointerfaces and is in the reviewing process.  of the analysis time (50 ps). Water selections were performed independently for each trajectory fragment. Each group contained a different number of water molecules: neighbour water, ~700 MGDG and ~1300 DOPC; intermediate water, ~500 MGDG and ~200 DOPC; far water, ~35 MGDG and ~1100 DOPC; Oc water, ~33 MGDG and ~40 DOPC; Op water, ~320; Og water, ~650; Nch water, ~1050. To test how the chosen criteria affect the results, the requirement and the observation time were modified and a water molecule was qualified to one of the groups if it belonged to the group for at least 140 ps (70%) or 190 ps (95%) of the 200-ps trajectory.
Mean square displacement (MSD) curves of water molecules belonging to selected groups are shown in Fig. 3 and Fig. 4. MSD curves for neighbouring, intermediate, far and bulk water were calculated in three dimensions ( Fig. 3a, b), two dimensions -in the x,y-plane ( Fig. 3c, d), and in one dimension -along the normal (z-axis) (Fig. 3e, f). Whenever it was possible, the diffusion coefficient was obtained from the linear part of the MSD curve, and is given in Table 1.

Diffusion of bulk water
The experimentally derived value of the self-diffusion coefficient of water, D xyz , at 25°C is ~23 ± 0.05 × 10 -6 cm 2 /s (Holz et al., 2000). The value of D xyz for bulk water at 22°C obtained in this study of 55.8 ± 2.3 × 10 -6 cm 2 /s is over twice larger than the experimental value. The discrepancy is due to the choice of simulation parameters, the water model, and the temperature control algorithm used in this study. It is known that values of D xyz for bulk water obtained from MD simulations depend on the simulation details and water models, and for the most popular water models they range from 27 to 59 × 10 -6 cm 2 /s at temperatures between 24 and 28°C (Mark & Nilsson, 2001), and from 24.9 to 51.9 × 10 -6 cm 2 /s at 25°C (Chaplin). In this study, the Nosé-Hoover thermostat (Hoover, 1985) and TIP3P water model were used and D xyz is 55.8 ± 2.3 × 10 -6 cm 2 /s. When using the Berendsen thermostat (Berendsen et al., 1984) and TIP3P water, D xyz at 37°C was 38 ± 1 × 10 -6 cm 2 /s (Rog et al., 2002) and was close to the experimental values of D xyz at 35 and 40°C of 28.59 and 32.22 × 10 -6 cm 2 /s, respectively. The Berendsen thermostat has proven to be an efficient algorithm to thermally equilibrate the simulated system. However, once the system is at equilibrium, the Nosé-Hoover thermostat should be used as it produces a correct kinetic ensemble. For this reason, in this and in our previous MD simulation studies of lipid bilayers, we have used the Nosé-Hoover thermostat to control tem-perature of the system at equilibrium. Additionally, one of the aims of this study was to explain differences in the lateral pressure profiles for the DOPC and MGDG bilayers and to establish correlation between pressure at a given depth of the bilayer and lipid-water interactions at the same depth. Therefore, in this study we were particularly careful to maintain the same simulation conditions as in the time-consuming productive simulations of MGDG and DOPC bilayers (a half microsecond long), used to calculate the lateral pressure profiles (Baczynski et al., 2015).

Diffusion of neighbouring, intermediate and far as well as bulk water
The entries in Table 1 for bulk water indicate that its translational diffusion is fully isotropic, whereas MSD curves in Fig. 3, particularly in Fig. 3e and f, show that the diffusion of water hydrating the bilayer is not isotropic and also that it is not a free Brownian diffusion. This is due to periodic boundary conditions (PBC) employed in the simulations and inhomogeneity of the bilayer system. PBC make the system continuous in the x,y-plane, whereas along the z-axis make the system periodic, that is, multi-lamellar. Water layers are 'locked' between bilayers, so water diffusion undergoes in the space confined by two, to a larger extent, impermeable walls. In effect, water diffusion in the x,y-plane is not or is only a little restricted but along the z-axis it is restricted. A detailed analysis of the diffusion of inter-lamellar water using MD simulation was performed by Sega et al. (Sega et al., 2005). They demonstrated that for the perpendicular diffusion (along the bilayer normal) MSD became non-linear after a relatively short time and later on reached a plateau. As can be seen in Fig. 3e and f, 1D MSD curves for intermediate and far water, after initial linear increase, start to level off. In consequence, the corresponding MSD curves for 3D diffusion ( Fig. 3a and b) become also non-linear. This prevented us from calculating D xyz and D z directly from the slopes of the MSD curves. Thus, for intermediate and far water, only numerical values for D xy were calculated directly from the respective MSD curves and are given in Table 1. Sega et al. (Sega et al., 2005) derived a formula enabling calculation of MSD for the perpendicular diffusion where D z was a parameter (eqn. 9; Sega et al., 2005), however, their analysis concerned the diffusion of all bilayer water molecules in the long-time asymptotic regime, thus, the formula is not applicable in the present study.
The non-linear behaviour of the 3D and 1D MSD curves for intermediate and far water might seem to be due to the restrictive requirement of 100% presence of a water molecule in the given group during the analysis time as this criterion certainly favours slower water molecules. However, the opposite applies as totally excluding faster water molecules from the analysis and including only slower ones cannot result in non-linearity of the MSD curves but merely in their smaller tilting. Thus, levelling-off of the MSD curves is due to the confinement of the diffusing water in the inter-lamellar space. The effect of including both faster and slower molecules was apparent when two less restrictive criteria, 70% and 95% presence and a longer observation time, were applied -indeed, the MSD curves were additionally levelled-off (not shown). Figure 3c and d, and the values of D xy for intermediate and far water in both MGDG and DOPC bilayers and for bulk water (Table 1) demonstrate that the lateral diffusion (in the x,y-plane) of the water is weakly affected by the bilayer. An apparent non-linearity in the upper fragment of the MGDG far water MSD (Fig. 3c) is mostly due to a small number of water molecules in this group, so 'environmental' effects are not properly averaged out during the analysis time. When a longer time scale analysis was attempted (200 ps instead of 50 ps, and 100% belonging to the group), it turned out that merely three water molecules fulfilled the criteria for belonging to the group. The lateral diffusion of intermediate water is very similar in the DOPC and MGDG bilayers. A large standard deviation of the determined D xy for far water in the MGDG bilayer, for reasons stated above, makes it difficult to compare the values of D xy in the MGDG and DOPC bilayers. However, one would expect that far water is the least affected by the bilayer and its D xy should be close to that of bulk water, which is the case for the DOPC bilayer.
Neighbouring water in both bilayers is over one order of magnitude slower than bulk water and its diffusion is almost isotropic (Fig. 3 and Table 1). This is in consequence of the fact that the neighbouring water molecules are mainly H-bonded to lipid oxygen atoms. During the bonding lifetimes, their motion is limited to a confined spatial volume within which they wobble.

Diffusion of H-bonded and clathrating water
MSD curves for water molecules H-bonded to the carbonyl oxygen atoms, Oc, in the MGDG and DOPC bilayers, to the galactose oxygen atoms, Og, in the MGDG bilayer, and to the phosphate oxygen atoms, Op in the DOPC bilayer as well as for water molecules clathrating the choline group, Nch, of DOPC are shown in Fig. 4 for diffusion in 3D (a, b), 2D (in the x,y-plane) (c, d) and 1D (along the membrane normal, z-axis) (e, f). Among these groups, the slowest translational diffusion is for Oc water. Op water diffuses two times faster than DOPC Oc water, whereas Og water diffuses almost one order of magnitude faster than MGDG Oc water. The diffusion of Nch water is similar to that of Og water but for both Nch and Oc water the diffusion is over one order of magnitude slower than that of bulk water.

Ordering of water molecules by the bilayers
H-bonding and clathrating are directional interactions that orient bonded water molecules relative to the interacting partner. To estimate the ordering effect of the bilayer lipids on water, an average orientation of the water dipoles at the given bilayer depth was calculated. To do that, the simulation box was divided into 1 Å thick slices along the bilayer normal and an average water dipole orientation in each slice was recorded. A water dipole orientation was calculated as cosine of the angle θ between the dipole moment of the water molecule and the bilayer normal. Each water molecule was assigned to a slice based on the position of its oxygen atom at each time frame (50 fs). Water dipole orientation profiles as functions of the bilayer depth were calculated for all water molecules in the MGDG and DOPC bilayers during 10-ns trajectories and are presented in Fig. 5.
The water dipole orientation profile for the MGDG bilayer differs from that for the DOPC bilayer. In the former one, the water dipoles in the water/bilayer interface have a tendency to orient along the bilayer normal (Fig. 5). In the DOPC bilayer, the water dipoles in the water/bilayer interface have a similar tendency also however, it is less pronounced. Moreover, in the upper interfacial region, there is a large fraction of water molecules oriented rather along the bilayer plane (Fig. 5). To support the general picture of water orientation in the bilayers drawn from the <cos(θ )> profiles, distributions of θ angles in three representative slices located -1.3 (deep interface), -2.2 (bilayer surface), and -3 (far  water) nm from the bilayer centre, c.f., Fig. 5, were calculated and are shown in Fig. 6. As can be seen from Fig. 6a, inside the interfacial region (slice -1.3 nm) of the MGDG bilayer, the distribution of θ angles is relatively narrow, ranges between 90 and 180°, with a tail towards smaller angles. The most probable orientation is 140° but population of angles larger than 120° occurs much more often than smaller angles. This indicates that the dipole moments in this slice make a relatively large angle with the bilayer plane. Inside the slice -1.3 nm in the DOPC bilayer water dipoles orientations have broader distribution which steeply drops towards larger angles and their most probable orientation is 120° (Fig.  6a). Population of angles smaller than 120° occurs more often than larger angles. Thus, the water dipoles in this slice are less inclined relative to the bilayer plane than the corresponding dipoles in the MGDG bilayer. At the upper edge of the interfacial region (slice -2.2 nm) of the MGDG bilayer (Fig. 6b), the distribution of θ angles is broader than in Fig. 6a, and almost symmetric about 90°, however, the number of angles larger than 90° is slightly greater than smaller angles. This indicates that most of the water dipoles in this slice are randomly oriented but a small fraction is oriented away from the plane. In the DOPC bilayer, the distribution of angles in the slice -2.2 nm has a flat maximum ranging between 50 and 90° and drops steeply towards smaller angles and smoothly towards larger angles. This results in the net orientation of water dipoles at small angle relative to the bilayer plane. In the slice -3.0 nm that is located above the bilayer surface, the distributions of orientations of water dipoles in both bilayers are uniform and symmetric about 90° (Fig. 6c). Thus, due to random orientations of the water dipoles away from the bilayer surface, the average cos(θ) in both bilayers is zero. Distributions of θ angles of water dipoles in the selected slices representing different locations along the bilayer normal (Fig. 6) agree perfectly well with the water dipole orientation profiles for the MGDG and DOPC bilayers, therefore they confirm trends apparent from mean values analysis.
To better understand the origins of the shapes of the water dipole orientation profiles, orientations of dipoles of the water molecules H-bonded to Op, Oc, and Og oxygen atoms of DOPC and MGDG and those clathrating the choline group of DOPC were calculated, separately for each interaction. Orientations were calculated for each slice from 50-ps trajectory fragments and averaged over 200 fragments (c.f. Methods). As was stated at the beginning of this section, orientation of water molecules hydrating the bilayers can be affected only by these interactions. The orientation (<cos(θ)>) profiles for bound water in the DOPC and MGDG bilayers are presented in Fig. 7a and b together with populations of orientations induced by a certain interaction ( Fig. 7c and  d). Each population is proportional to the number of water molecules bound by the given group in the given slice. The Oc oxygen atoms in both bilayers bind the smallest number of water molecules as they are buried most deeply in the bilayer interface, whereas the choline groups of DOPC bind the largest number of water molecules. A surprizing result is that the ranges of the distributions in Fig. 7c and d, to a large extent overlap indicating that deeply in the bilayer interface of the DOPC bilayer there are water molecules H-bonded to Oc and Op oxygen atoms and also clathrating the choline group. The same applies to the Oc and Og water in the MGDG bilayer. Nevertheless, Oc water reaches out much less into the upper interfacial region than Op, Og and Nch water, and Op water does not reach as deeply into the interface as Oc and Nch water. Another interesting result is that the net orientation of dipoles of water molecules bound by each group changes sign at some interface depth. This means that an average orientation of dipoles relative to the bilayer plane changes ( Fig. 7a and b); the effect is the least pronounced for Og water (Fig. 7b). In each case, the change of sign takes place near the maximum of the distribution of orienta-

Figure 7. Profiles of an average dipole orientation along the bilayer normal for Oc (red bars) and Op (blue bars) H-bond water and Nch (green bars) clathrating water in the DOPC bilayer (a), and Oc (red bars) and Og (green bars) H-bond water in the MGDG bilayer (b). Distributions of numbers (Occurrence) of a given average dipole orientation along the bilayer normal for Oc (red bars) and Op (blue bars) H-bond water and Nch (green bars) clathrating water in the DOPC bilayer (c), and Oc (red bars) and Og (green bars) H-bond water in the MGDG bilayer (d).
tions ( Fig. 7c and d), which, on the other hand, corresponds to the maximum of the density distribution of the given group across the bilayer (unpublished result). As the number of Oc water molecules is small, their orientations have less impact on the properties of the bilayer than the remaining bound water molecules. For the Op, Og, and Nch water, the change of sign occurs near the distance of ~2 and ~2.2 nm from the bilayer centre of the DOPC and MGDG bilayer, respectively, which closely corresponds to the distance where the density of the head group region drops to ½, which on the other hand, is considered to specify the bilayer surface. As  Figs 6b and 5 indicate, the distributions of orientations of water dipoles, measured both as <cos(θ)> and θ, in the slice -2.2 in the DOPC bilayer are not fully random and an average dipole orientation makes a small angle with the bilayer surface (~10° or less), whereas those in the MGDG bilayer, are almost random and cos(θ) averages out to nearly zero. At distances between 2.2 and 3.0 nm, DOPC choline groups still affect orientation of water molecules, whereas Og groups of MGDG practically do not.

DISCUSSION
In this article, the effect of two bilayers, one comprising lipids with galactose head group (MGDG bilayer) and the other with phosphatidylcholine head group (DOPC bilayer), on the dynamics and orientational order of the water molecules hydrating the bilayer was analysed. The bilayers affect both translational motion of water molecules and orientation of their dipole moments relative to the bilayer normal. In general, the diffusion of the bilayer water is slower than bulk water (Rog et al., 2009) and is affected by the confinement in the inter-lamellar space (Sega et al., 2005). In detail, water molecules directly interacting with lipid groups diffuse significantly slower than unbound water molecules. But even within the group of bound water there are differences. The diffusion of Oc water, which is buried deep in the interfacial region, is the slowest. Op water is located higher in the interfacial region of the DOPC bilayer and diffuses two times faster than Oc water. The Og and Nch water molecules that interact with the lipids terminal groups diffuse several times faster than water molecules Hbonded to Op and Oc oxygen atoms. An unexpected result of this study is observation that water molecules Hbonded to Og in the MGDG bilayer diffuse similarly to those clathrating the choline group in the DOPC bilayer. This might be explained by similar mobility of the galactose and the choline moieties to which Og and Nch water molecules are bound. It is also possible that exchange of Og and Nch water molecules between neighbouring galactoses and choline groups, respectively, is faster than 50 fs, the trajectory sampling interval.
It is interesting to note that the effect of the bilayer lipids on the lateral diffusion (in-plane) of water extends not farther than 7 Å or rather not much further than 4 Å away from any atom on the bilayer surface. This result contradicts our previous results (Rog et al., 2002) as well as those of Murzyn and coworkers, (Murzyn et al., 2006) where the effect of water confinement on the lateral diffusion was overlooked.
The ordering effect of the bilayer on water dipoles extends from the bilayer centre up to ~3 nm in the DOPC bilayer and up to ~2.8 nm in the MGDG bilayer. At larger distances orientations of water dipoles are random in both bilayers. Thus, the effect of the bilayer lipids on the orientation of dipole moments of water molecules is limited predominantly to groups directly interacting with water molecules; in this study they are Oc, Op, and Og oxygen atoms of MGDG and DOPC, and choline groups of DOPC. An average dipole orientation relative to the bilayer plane of bound water changes near the bilayer surface, which indicates that near the surface, the net orientation of water dipoles is close to horizontal.
Net orientations of water dipoles at different depths of the MGDG and DOPC bilayers correlate very well with lateral pressure profiles calculated in our previous study (Baczynski et al., 2015). We suggested there that higher lateral pressure at the upper edge of the interfacial region of the DOPC than the MGDG bilayers might result from a better ordering of water molecules in this DOPC bilayer region. In the present study it is shown that indeed, ordering of water in the DOPC bilayer extends further into the water interface than in the MGDG bilayer. This might lead to higher local expansion of the bilayer and a more positive value of the lateral pressure. In effect, the lateral pressure above the surface of the DOPC bilayer is higher than it is the case in the MGDG bilayer.