Analysis on the Water Retaining Capacity of Membrane by Molecular Dynamics Simulations∗

Lipid membrane has a significant role to retain waters on its surface due to high hydrophilicity of head groups of lipid molecules. Although many physical chemical properties such as thickness, phase state, fluidity etc. can be measured by the recent experimental equipments, the water retaining capacity of membrane is still difficult to be analyzed. Molecular dynamics simulation was performed to theoretically estimate the retaining time of waters on membrane surface. Different kinds of lipid bilayers were modeled using sphingomyelin molecules, which is one of the typical components of biomembrane. Our simulation successfully generated the stable and equilibrated bilayer structures of sphingomyelin with different carbon chain length. Thickness, area per lipid, order parameter, and gauche ratio were deduced from the simulation data. In this work, we propose a technique to estimate the distribution function of waters on membrane surface, which enables us to evaluate the water retaining time on membrane. Our calculation technique will be helpful for design an artificial membrane whose moisture property is controlled. [DOI: 10.1380/ejssnt.2009.591]


I. INTRODUCTION
Biomembranes are composed of several kinds of lipid molecules which are classified into three families; phospholipids, sphingolipids, and cholesterols.Typical molecular species of them are phosphatidylcholine, sphingomyelin and cholesterol, respectively.The properties of biomembrane substantially depend on the composition of lipid molecules.For the purpose of measuring the physical-chemical properties of membrane, many approaches have been developed so far.The critical temperature for the phase transition between gel and liquidcrystal or the fluidity of lipid molecules at a defined temperature can be monitored by the recent thermal or spectroscopic equipments.In spite of such technical developments, the water retention capacity of membrane has not been well analyzed.This analysis requires the understanding of the movement of water molecules at the atomistic scale, which is difficult to be detected by the current in vitro experimental techniques.
Molecular dynamics (MD) simulation is one of the most popular theoretical approaches to make the motion in biological phenomena visible at an atomistic or molecular level.This method enables us to imitate an actual in vivo condition in a computer.For better structural insight into lipids, many simulation studies have been performed on the biomembranes with lipid bilayer [1,2] and the channels inside the lipids [3].For example, Scott et al. performed MD simulation of ceramide bilayer for the first time and they introduced order parameter, electron density profile and so on [4].Other groups carried out the simulations on phospholipids [5] or mixed membrane [6].Recently, a computational technique to analyze the interaction between protein and lipid membrane has been improved because more complex computational system such as a protein embedded in a mixed membrane and a protein attached to a membrane surface can be modeled [7,8].Therefore, the atomistic scale understanding of the action of membrane proteins will rapidly advance with the assistance of computer simulation.A standard method for describing the physical properties of biomembrane, however, has not been established sufficiently.Hence, a variety of analyzing approaches should be developed for interpreting experimental data precisely and/or for predicting biological phenomena reliably.
In this study, we have figured out a novel analyzing technique that can measure the retaining time of waters on different kinds of membranes.Our analysis successfully deduces the time for a water molecule to be kept on a membrane by statistically dealing with MD simulation results.This retaining time is difficult to be measured experimentally.Indeed, it is impossible to monitor the motion of the respective water molecules, even though the time-resolution X-ray crystallography is developed to analyze the structural change of molecules in the hundreds pico-second order [9].In the field of computational study, some reports introduced the radial distribution function (RDF) of waters with respect to a particular atom [10].Unfortunately, the conventional RDF cannot be applied to the waters on a membrane because water molecules on a membrane are not distributed spherically.Therefore, we will propose a modified distribution function (DF) of waters with respect to the membrane surface in order to calculate retaining time of waters; i.e., how long water molecules stay on the membrane.This function is different from the conventional RDF that represents particle distribution against a particle.Instead, our present approach is intended to describe particle distribution against a plane.
The retaining time of waters will be affected by the lipid structure, so we modeled 6 different kinds of membranes with sphingomyelin bilayer, changing the length of acyl chain.That is, 12:0, 14:0, 16:0, 18:0, 20:0, and 18:1 (cis ∆9) sphingomyelins are prepared as computational membrane systems.The backbones of these acyl chains are derived from lauric, myristic, palmitic, stearic, arachidic, and oleic acid, respectively, which are free fatty acid mainly included in animal body.In this study, we carried out modeling of sphingomyelin membranes, execution of MD simulations, subsequent analysis of the simulation results, observation of a common feature of sphingomyelins, and estimation of the retaining time of waters.

II. METHOD
A. Modeling of membranes VMD [11] was used to make a model of phosphatidylcholine (POPC) bilayer that has 30 POPC molecules per half leaflet; i.e., totally 60 molecules in bilayer.POPC molecules were then converted into sphingomyelin (SM) by modifying POPC head and tail groups.Thus we acquired 6 SM bilayer models composed of single molecular species in which the length of acyl chain varies from, 12:0,14:0, 16:0, 18:0, 20:0 to 18:1.An unsaturated bond in acyl chain has cis conformation and is located between 9 C and 10 C. TIP3P explicit water model [12] was used to solvate each membrane system, generated above and below the bilayer with the thickness of 40 Å.Sodium and chloride ions were added into the water layers in the concentration of 0.15 mol/L to reproduce the biological condition.As a result, the number of total atoms in each model system was 23629 for 12:0 SM, 23989 for 14:0 SM, 24349 for 16:0 SM, 24349 for 18:0 SM, 25069 for 20:0 SM, 24589 for 18:1 SM.

B. MD simulation
NAMD package [13] was employed to perform all processes of simulation.Energy minimization was initially carried out for 10,000 steps to release the steric hindrances.Next, the model system was heated to 310 K for 6,000 steps.Simulations for equilibration were then carried out under the NPT ensemble condition; 310 K and 1 atm.The head group of lipids was initially constrained with the harmonic constraint method in order to maintain the structure of the whole system.The simulation was executed for 2 ns where the constraints were weakened gradually and subsequently for 2 ns where the constraints were all free.Finally, 10 ns simulations were executed for 6 SM models as a production run.CHARMM 27 parameters [14] were used for the force field of lipids and CHARMM 22 [15] was used for the others.The long-range electrostatic interaction was calculated at every step by the particle mesh Ewald (PME) method under the periodic boundary condition [16].The integration time step was set to 2 fs.All covalent bonds connecting hydrogen atom were constrained by the SHAKE algorithm or the SETTLE method for water [17,18].Langevin thermostat method was used to control the system temperature [19].Short-range non-bonded interactions were cut off smoothly between 10 and 12 Å.

C. Analysis
All the results shown in Tables I and II and Figures 3  and 4 are produced by the analysis on the production run for the last 5 ns of each MD simulation.Thickness was extracted from the distances among the phosphorus atoms of phosphate group of SM.Area per lipid was obtained by dividing the size of unit cell by the number of molecules; i.e., 30.Deuterium order parameter indicating the degree of ordering of lipids was calculated from the Eq. ( 1); where θ represents the average angle between the normal vector of membrane (z axis) and each vector of the C-H bonds of CH 2 groups of lipids.The brackets denote the averaging over time and over all lipids.
The radial distribution function (RDF) is usually applied to describe how many objective particles exist around a particular particle.RDF of objective particles in a space is derived from the Eq.(2); where r represents the distance from a particular particle, n(r) represents the number of objective particles in ∆r between r and r + ∆r, and ρ represents the density of objective particles at the infinite separation from the particular particle.We modified this equation to be applicable to objective particle on a surface, in which r is changed to represent the distance along the surface normal measured from the bilayer surface and n(r) represents the number of objective particles per unit area.Figure 1 shows the scheme of algorithm we figured out for calculating n(r) of waters that are distributed from membrane surface to bulk-water layer.The distance between each water molecule and the lipid atom closest to the water molecule was calculated and was presented as a histogram.This histogram shows the distribution of water molecules when all lipid atoms closest to any water molecules were settled to the position at the top surface (x-y plane) and all water molecules were rearranged at the distance r along the direction normal to the membrane (z axis) (Fig. 1(d)).Thus the acquired plot (Fig. 1(e)) represents how many water molecules exist at a particular distance from the surface, for example, the plot usually has a peak at a distance of hydrogen bonding.
The distribution function (DF); g(r), proposed here is useful to understand the interaction between biomembrane and waters.The time for water molecules to stay on the membrane interface or surface can be calculated by solving the Eqs.(3)-(5) [20]; where k B is Boltzmann constant, T is temperature (310 K), W is the potential energy with respect to the mean force on a water, ∆E is the energy necessary for transferring the water from membrane to bulk, h is Planck's http://www.sssj.org/ejssnt(J-Stage: http://www.jstage.jst.go.jp/browse/ejssnt/) e-Journal of Surface Science and Nanotechnology constant, and τ is the time for passing through the energetic barrier.From the static viewpoint, the capacity for water retaining might be determined by the number of water molecules distributed around membrane.Instead, we define the capacity for water retaining from the dynamic aspect.Namely, based on the molecular theory of solvation, the water retaining is deduced from the kinetic value such as the time necessary for one water molecule to transfer across the interface between membrane and bulk water layer.This membrane-bulk water layer interface will correspond to the barrier at the first maximum in Fig. 1(e).

A. Structural Property
Top and side views of 6 kinds of sphingomyelin bilayers after 10ns MD simulation are shown in Fig. 2. The structures are different among 6 kinds of bilayers, in spite that the atom geometry in head groups were set similar to each other in their initial structures.Then, this difference is a consequence of the effect of hydrophobic interaction due to the difference in length of acyl chain in tail group.
Some major analysis on lipid membrane were performed; i.e., thickness, area per lipid, and tendency for gauche conformation are summarized in Table I.The order parameter profiles are presented in Fig. 3.All results in Table I were obtained by sampling the structure every 10 ps for the last 5ns and averaging over the values for the 500 snapshot structures.14:0 SM, whose two carbon chains are same in length, are compactly aligned and the terminal carbon atoms are well-ordered .The difference between two carbon chain lengths induces disordering at the center of membrane.18:1 SM are also loosely packed due to cis-double-bond.These factors are reflected in the results of gauche ratio and area per lipid.
When determining the phase state of membrane experimentally, the gauche ratio acquired from IR data is used in a standard technique, in which the ratio is low in the gel ordered phase whereas high in the liquid crystalline disordered phase.In this report, gauche ratio was measured, following the definition by Zaraisukaya et al [21].Koynova et al. collect the data of the transition temperature (Tm) of sphingomyelin having saturated lipid acid [22].According to the data, Tm increases with chain length.Around 310K, the lipid bilayers consisting of the molecules with 16 or shorter carbon chain length form not liquid-crystal phase but gel phase.This experimental findings is consistent with our MD simulation results   except for 18:0 SM.

B. Water retaining property
Radial distribution function (RDF) is often used as a fundamental approach for measuring the capacity of water retaining for metals or ions.The RDF represents the relative density distribution of objective particles surrounding a particular particle (e.g.atom).In this work, instead of RDF, we figured out the distribution function (DF) of objective particles against a plane as shown in Method section.If the DF of water molecules is available with respect to membrane, it quantitatively suggests the capacity of water retaining or moisturizing.With a simple algorithm as shown in Figure 1, we have successfully obtained the DF of waters against biomembrane.In biomembranes, the transient layer from lipid to water is in a range of 2 to 5 Å (Fig. 4, black line).Although RDF of objective atoms against a particular atom usually shows the first and second hydration peaks, the DF for membrane gives a small maximum like shoulder after the first maximum, followed by a shallow minimum and no prominent second hydration peak.Most of the water molecules gathering at the first maximum are close to oxygen atom of SM and located at 2.6-2.8Å away from the closest oxygen atom of lipids.That is, this peak is corresponding to the hydro- gen bonding with oxygen atom of SM and the peak position matches the hydrogen bond distance.As for water molecules at the shoulder-like small peak (around 3 Å), approximate half of the atoms closest to the waters in lipids are carbon.From these observations, we separated the DF for membrane into two parts from the closest atom of lipids; oxygen and carbon, and presented in blue and red lines of Fig. 4, respectively.In all cases of sphingomyelin bilayers, almost all water molecules at the first maximum are close to oxygen atom of lipid molecules and those at the shoulder peak are confirmed to be close to carbon of lipids.It is surprising that DF of water molecules with respect to carbon atom have a peak in case of membrane.It should be also noted that the distance between the water molecule and carbon atom was different from that between water and oxygen.DF of 18:1 SM has a large shoulder peak as seen in Fig. 4. In DF of membrane, the density of water molecules becomes constant far away from membrane and, then, the waters beyond  We estimated the retaining time of waters on membrane from the calculated DF.One water molecule can stay on the membrane surface for some period and this value is regarded as an index for the capability of water retaining (Table II).The estimation indicated that it takes 30 to 36 fs for one water molecule on the membrane surface to transfer to bulk water layer.Although, for example, the time-resolved X-ray diffraction technique has been developed and is now expected to detect a femtosecond phenomena, the time reported here can not be provided easily even from the state-of-the-art experimental approaches of today.The difference of the length of acyl chain of lipids apparently effects the water-retaining time.As the length of the acyl chain increases, the membrane is likely to retain water longer.The effect of carbon atom to hold waters is enhanced with the increase of the unbalance in length of two acyl chain.There are many small cavity on the membrane surface.Much water molecules are trapped inside the cavity and not rapidly transfer to the bulk layer.We speculate the formation of those cavity is an origin for the secondary peak that the closest lipid atom to the trapped water is carbon.No experimental observation in the relationship between lipid carbon chain and capability of water retaining has been reported.It will be necessary to carried out the simulation with other types of lipid molecules such as phosphatidylcholine and ceramides to clarify the influence of the difference of chemical structure of the head groups.

IV. CONCLUSIONS
We carried out MD simulations for 6 kinds of membranes with sphingomyelin bilayer changing the length of acyl chain for 12:0, 14:0, 16:0, 18:0, 20:0, and 18:1.The thickness of membrane increases with the length of acyl chain, except for 20:0 and 18:1.The sphingomyelin membrane with 20:0 and 18:1 show the interdigitation of acyl chains at the center of bilayer.Owing to the interdigitation, the area per lipid becomes considerably high for membranes with 20:0 and 18:1.
In this study, we proposed a calculation technique to obtain the distribution function of waters on the membrane surface.The distribution function can be separated into two components.One component is due to the waters whose closest atom of lipid is oxygen.The corresponding peaks appear at 2.6-2.8Å and this distance is compatible with the hydrogen bond distance .Another component is due to the waters whose closest atom of lipid is carbon.The corresponding peaks appear around 3 Å.Most of the water molecules in the latter component are trapped inside the small cavity on the membrane surface.The retaining time of waters on membrane, that is the average time for one water molecule to transfer from membrane surface to bulk water layer, was estimated from the distribution function of waters.The estimation value of the retaining time of waters on lipid membrane is 30-36fs.

FIG. 1 :
FIG. 1: Algorithm for estimating the distribution function (DF) of waters against biomembrane.(a) Algorithm 1, where the number of water molecules between r and r+∆r is counted, referring the z-coordinate of water O atom.The open and solid circles denote waters and lipids atoms, respectively.(b) Typical DF profile of waters along the membrane normal (z axis) using Algorithm 1.No peak appears.(c) Algorithm 2, where the distance between a water molecule and its nearest lipid atoms is measured.(d) All the waters are rearranged along the direction normal to the membrane.(e) Typical DF of waters obtained by Algorithm 2. A critical peak appears.

FIG. 2 :
FIG. 2: Snapshot structures of several kinds of Sphingomyelin membranes (upper: top views, lower: side views).Carbon, Oxygen and Nitrogen are represented in green, red and blue stick, respectively.

4 :
Distribution function (DF) profiles of waters on the surface of sphingomyelin bilayer.DF for all water is denoted by black solid line.DF for water whose closest atom of lipid ozxygen and DF for waters whose closest atom of lipid is carbon are colored red and blue, respectively.