Molecular dynamics of water transport through membranes: Water from solvent to solute

An application of Molecular Dynamics computer simulation (MD) to the process of transport of water through a lipid bilayer membrane is described. The permeation process is far too slow to be modeled by straightforward MD. In stead the inverse of the permeability coefficient is expressed as an integral over a local permeation resistance, which itself is inversely proportional both to the local density and the local diffusion constant. These quantities are recovered from MD simulations. The local density relates to a free energy profile, which is constructed by a combination of density determination, of the mean force on constrained molecules, and particle insertion. Thus a slow process can be accurately predicted from relatively short MD simulations. 1. I N T R O D U C T I O N This presentation concerns a computer simulation study, using molecular dynamics (MD), of the transport of water through a bilayer membrane, in this case consisting of dipalmitoyl phosphatidylcholine (DPPC, lecithin). On the one hand this process is of considerable biological interest and on the other hand it poses a methodological challenge. The challenge is twofold. The first challenge is to perform a reliable MD simulation of a lipid bilayer in its liquid crystalline state using atomic detail of lipid molecules and solvent. Section 2 reviews such simulations as they have been performed in our laboratory over the past decade.The second challenge is to model a process (water penetration) that is very slow on the time scale of the simulation: during a simulation of typically one hundred picoseconds, it may happen once or twice, or not a t all, that a water molecule will wander through the bilayer.Thus no statistics a t all can be obtained on the permeation process. In section 3 the permeation process is described in the diffusion limit in terms of the thermodynamics of irreversible processes; this yields the key to computation of the permeability coefficient by combining the depth profiles over the membrane of free energy and diffusion. Sections 4 and 5 explain how these profiles can be determined from MD simulations. In Section 6 the results of both free energy and diffusion simulation are combined to obtain the permeation coefficient. A full account of this research is given in ref. [l]. 2. M D S I M U L A T I O N OF LIPID BILAYERS Simulation of lamellar liquid-crystalline systems has been a subject of research in our laboratory for over a decade, starting with a simplified bilayer of decane chains without detailed modeling of head groups and solvent [2,3], via detailed dynamics of multilamellar mesophases of sodium decanoate/decanol/water systems [4-61 to atomic simulations of DPPC in water (6,7]. The nature of the 'hydration force' between DPPC bilayers has also been investigated [8] and found to be related to the ordering of water by the interfacial dipolar charges. 2513 2514 H. J. C. BERENDSEN AND S:J. MARRINK The stability of the bilayers and their aggregation phase depend critically on details of the force field used. Lecithin bilayers can exist in (at least) two phases: a low-temperature gel phase in which the hydrocarbon chains are ordered but the head groups are still fluid, and a liquid-crystalline or phase in which the hydrocarbon chains are also fluid. With the standard GROMOS force field [9], DPPC was found to equilibrate towards the gel state, even at a temperature of 350 K, well above the phase transition temperature between gel and liquid-crystalline state. It waa found necessary to make two modifications in order to obtain the liquid-crystalline phase at 350K: use a more accurate description of the dihedral angle function in the hydrocarbon chains, and reduce the ionic charges (on phosphate and choline) by a factor of two. The latter has the same effect aa using a dielectric constant of 4 for the interactions of charges; this compensates partly for the lack of explicit polarisability in the force field. The DPPC simulations and its force field aspects have been recently reviewed [lo]. The simulations were carried out in a rectangular periodic box using an algorithm that effects a coupling to a 'bath' of constant pressure and constant temperature [ll]. The dimensions of the box adjust themselves to be consistent with the set pressure (1 atmosphere). While the overall density quickly equilibrates, the depth of the box, related to the thickness of the bilayer, takes several hundred picoseconds to equilibrate. It is in fact the slowest and most sensitive indicator of equilibrium, requiring the equilibration of rotational chain isomerism which takes place on the time scale of tens of picoseconds. It was found that either a gel or a liquid-crystalline phase could be obtained, depending on force field and temperature. The twodimensional radial distribution function of chains showed a much stronger evidence of hexagonal packing in the gel than in the liquid-crystalline phase. The equilibrium surface areas for the two phases are in good agreement with experimental data: 0.485 nm3 for the gel phase (experimental value 0.48) and 0.575 nmz for the liuid-crystalline phase (experimental value 0.60). Another sensitive comparison with experiment can be made by computing the deuterium resonance order parameters which can be obtained from deuteron magnetic resonance line splittings. The order parameter SCD that is measured is an average over the second spherical harmonic function of the angle B between the direction of the CD-bond and the membrane normal: Such order parameters can be obtained for aJl carbon atoms of the chain except the terminal ones. For a fully ordered, extended, non-tilted chain, the order parameter attains the value of -0.5; for an isotropically disordered chain its value is zero. Order parameters typically show a plateau value for the first half of the chain, tapering off to lower values towards the tail end of the chain. Well-ordered model systems like sodium decanoate/decanol/water have plateau values of about -0.3; in the less well-ordered lecithin bilayers the order parameters are around -0.2 in the liquid-crystalline phase We have found almost exact agreement with experiment both for the model system simulations [5] and the lecithin simulations (6,7]. This gives confidence in those results that cannot be tested by experiment. Fig. 1 shows a projection of a typical snapshot of the lecithin bilayer. There is more disorder than biochemistry text books suggest. The head group region spreads over about 0.7 nm and the interface has a fluctuating corrugated structure. A film conveys the fluctuations in the structure much better than a picture of this kind can do1. Four different zones, indicated in Fig. 1, can be distinguished: Region 1 Here the head groups start to appear, but still a t low density. The water density decreases to about half the bulk value. Region 2 This is the region where the head group density is high and the water density decreases to almost zero. Hydrocarbon chains begin to appear. Region 3 This is the well-ordered and rather dense region of the hydrocarbon chains. The density here drops from 1.1 to about 0.75 g/cm3, typical for liquid decane. Region 4 In this tail region with its high proportion of methyl end groups, the density drops to about 0.6, typical for liquid ethane. 'A copy of L 7 minute video film can be obtained for US S 40,from BIOMOS bv, Nijenborgh 4, 9747 AG GRONINGEN, the Netherlrndi Molecular dynamics of water transport through membranes 2515 Figure 1: Snapshot of a lecithin bilayer. Choline and phosphate groups are highlighted, Water molecules are dashed V-shaped particles The density profile in terms of electrons per unit volume agrees with the experimental profile determined from x-ray diffraction measurements on multilamellar phases. The density in the tail region is quite low and allows for a few percent free space for a probe molecule of the size of a water molecule [l]. 3. TRANSPORT THROUGH MEMBRANES: A THERMODYNAMIC DESCRIPTION During a simulation of typically 100 ps, it occurs once or twice that an isolated water molecule 'dissolves' into the membrane hydrocarbon phase, where it can diffuse quite easily. We never observed a cluster of two or more water molecules in the interior of the membrane, detached from the bulk aqueous medium. Thus it seems that the major transport of water takes place by permeation of single water molecules, although a very infrequent occurrence of permeating clusters or even of a complete pore by thermal fluctuation cannot be excluded on the basis of these simulations alone. It is clear that water permeation cannot be studied quantitatively by observing this process in an equilibrium simulation. Therefore we resort to a different approach. The permeation process is considered as a diffusion process in an external field, described as a free energy profile or potential of mean force. Our first concern is to describe the diffusional theory of transport in terms of the thermodynamics of irreversible processes. This enables us to link the permeation coefficient to experiments as well as to an integral over local properties in the inhomogeneous membrane. 2516 H. J. C. BERENDSEN AND S:J. MARRINK 3.1 General diffusional theory of transport We consider the motion of particles of the i-th species (in this case water, but the theory is applicable to any solute as well) in the diffusional limit, where the average velocity u, = (a,) is proportional to the thermodynamic dn'uing force, which is the negative gradient of the thermodynamic potential p;: where C; is the frictional coefficient of the particles. [Note that, if the thermodynamic potential is expressed in J/mol, the frictional coefficient is expressed per mole of particles as well]. The flux J, in mole m-' s-' is given by (2) J , . . ' , c,u, = vpi ci The friction coefficient C, is related to the diffusion constant D, by as can be easily seen when a concentration gradient in an ideal solution is considered for which ~1 = + R T l n c and eq. 2 is compared to Fick's law J = -DVc. The linear flux relations for the case that material properties depend on one coordinate z can now be written as Together with a continuity equation, eq. 4 predicts the spatial and temporal evolution of the local density distribution. We are, however, interested in the steady state solution of the flux in the linear regime, i,e., under the influence of a small deviation from equilibrium. Steady state means that J; is not a function of z, and after rearranging we can integrate eq. 4 over the membrane from 21 in the bulk phase on one side to zz in the bulk phase on the other side: Here c,(z) is the concentration of component i in the presence of the imposed gradient. Under the assumption of small gradients we can replace this concentration by the equilibrium Concentration c y in the absence of an imposed gradient. If we define the permeation resistance RP as where cT is the concentration in the bulk solutions on either side of the membrane in the absence of an imposed gradient, the linear response relation, eq. 5 , becomes The permeation resistance is directly related to the experimental permeability coeficient (section 3.2), and is also amenable to computation on the basis of detailed simulation. For the latter the calculation of the local equilibrium density C?(Z) (section 4) and the local diffusion constant (section 5 ) are separate tasks that must be accomplished first before eq. 6 can be integrated. The result of the integration is given in section 6. Before proceeding, let us summarise the assumptions that have been made either explicitly or implicitly in the derivations thus far: (1) The whole system is isothermal at absolute temperature T (2) The membrane component is stationary in the frame of reference (3) The local diffusional model is valid: the thermodynamic gradient can be considered constant over the correlation distance of the particle (the distance given by the displacement of a particle during the time over which its velocity correlation function differes from zero) (4) The fluxes are proportional to the gradients in thermodynamic potential. This means that the limit of small gradients is considered where this is needed. ( 5 ) The permeation process for water is dominated by single molecules that only feel friction with the stationary membrane component. Molecular dynamics of water transport through membranes 2517 Assumption (3) is the most questionable one, because the concentration gradients in a membrane are very large and the diffusion is relatively fast. It is possible to refine the barrier-crossing dynamics by including details of the velocity (or force) autocorrelation of the particles [12], but we restrict our considerations to the simple diffusional limit. 3.2 Experimental quantities The driving force for permeation processes can be imposed by the following causes: hydrostatic pressure difference Ap (water), osmotic pressure difference All (water), or concentration difference Ac (solute or water isotope). Hydrostatic and osmotic differences are equivalent in their influence on the thermodynamic potential of water: Comparing eq. 7 with eq. 8, the flux can be expressed as A/lw = (Ap AII)/cL (8 1 J 1 (Ap-AII) RP, RT w (9) For the flux J, of an isotope i of water, we consider the z-dependent mole fraction zi(z) of the isotope. Its thermodynamic potential is given by /li(z) = / l w ( ~ ) + RTln xi(.) (10) Assuming water to be in equilibrium over the membrane, pw is constant and equal to its bulk value &,.Integration of eq. 4 using eq. 10 and equating ci(z) with zi(z)c,(z), we find 1 J; = --Aq, RP, where Ac; = ctAx; is the concentration difference of the isotope across the membrane. For the flux of a solute resulting from a (small) concentration difference Ac; over the membrane, for which it is easily derived that