Translational Diffusion in Lipid Bilayers: Dynamic Free-Volume Theory

Dec 29, 1998 - those obtained from free-flight-length distributions of solute, a novel measure of ... Molecular dynamics simulation can be used to cal...
0 downloads 0 Views 168KB Size
Download PDF
J. Phys. Chem. B 1999, 103, 385-394

385

Translational Diffusion in Lipid Bilayers: Dynamic Free-Volume Theory and Molecular Dynamics Simulation Tian-xiang Xiang Department of Pharmaceutics and Pharmaceutical Chemistry, UniVersity of Utah, Salt Lake City, Utah 84112 ReceiVed: September 10, 1998

Position- and direction-dependent diffusion coefficients of spherical, nonpolar solutes as represented by noble gases in a model lipid bilayer are investigated by nonequilibrium and equilibrium molecular dynamics (MD) simulations. The nonequilibrium MD method, which is based on linear response theory, is developed by determining an optimal range of external forces in which the response of solute mean velocities is linear and the perturbation to bilayer properties is minimal and by comparing with equilibrium MD simulation and experiment. Solute diffusivity is found to be smaller in a more ordered acyl-chain region and in the plane normal to the bilayer interface and exhibits markedly different dependences on solute size, depending on location in the bilayer interior and diffusion direction. The diffusion coefficients calculated closely follow those obtained from free-flight-length distributions of solute, a novel measure of free-volume anisotropy, based on the Einstein relation, suggesting that free volume is largely responsible for molecular diffusion in lipid bilayers. Theoretical formulas for solute diffusivity in lipid bilayers where the free-volume redistribution is governed primarily by overall rotation and local isomerization of lipid molecules is proposed in the framework of dynamic free-volume theory. This theory predicts variable degrees of changes in molecular diffusivity with solute size and membrane structure depending on relative contributions of overall rotation and local isomerization of lipid molecules to free-volume redistribution and can explain the observed size dependences provided that lipid overall rotation (local isomerization) plays a major role in determining transverse (lateral) diffusivity of solute molecules.

Introduction One of the major functions of biomembranes is to regulate bioactivities and transport of various molecular species in living cells. Thus, a thorough understanding of molecular mechanisms for solute diffusion in lipid bilayers, one of the major components of biomembranes along with proteins, is crucial to elucidate structure-activity relationships for many biological processes (e.g., anesthesia, drug potency and toxicity, etc.). Unfortunately, there are still tremendous difficulties to directly observe translational motion of small molecules within various regions of biomembranes, as demonstrated by a wide variability of membrane microviscosity observed by different experimental methods.1-3 As in polymers and liquid crystals,4,5 our current analyses of molecular diffusivity in biomembranes are mostly based on free-volume theory.6-10 Since the pioneer work by Cohen and Turnbull in the late 1950s,11 free-volume theory has undergone remarkable improvement.12-14 In fact, attempt has been made to render free-volume theory for self- and mutualdiffusion coefficients in polymers completely predictive, so that all model parameters can be estimated without knowledge of any diffusion data.15 However, previous free-volume models fail to account for the influence of free-volume mobility on solute diffusion and the contributions of different motional modes of solvent to free-volume redistribution. This is in contrast with the kinetic theory of Enskog16 and an empirical model proposed by Bueche17 for diffusion in polymers in which the effects of relative motions of solute and solvent molecules and polymer relaxation rates are considered, respectively. In fact, a diffusion Deborah number has been proposed by Vrentas et al. in amorphous polymers for anticipating conditions under

which the viscoelastic relaxation dynamics of polymers becomes important,18,19 though no general theory has been formulated. The importance of solvent relaxation dynamics in determining solute diffusivity is highlighted in recent molecular dynamics (MD) studies,20,21 in which the increase of transition energy barriers for chain isomerization considerably reduces solute diffusivity in polymers and membranes while the static free volume remains unchanged. A novel dynamic free volume was recently developed by Xiang21 to interpret these findings. Molecular dynamics simulation can be used to calculate various transport properties which are difficult to obtain experimentally and to explore molecular mechanisms underlying these transport properties. Despite these advantages, only recently have certain aspects of molecular diffusion in lipid membranes been investigated by means of MD simulation21-27 because of the enormous amount of computation resources required. A more severe obstacle to calculating molecular diffusivity is the heterogeneous nature of lipid membranes. Recent statistical mechanics and MD studies from this and other laboratories revealed a strong preference for hydrophobic solutes to move into the more disordered region near the center of the bilayer.28-30 This is responsible for the failure of the meansquared displacement of solute at a long time to follow the wellknown Einstein relation,31 as demonstrated by Huertas et al.32 using a Brownian dynamics method with a simple potential model for solute motion in lipid bilayers. An attempt has recently been made by Marrink and Berendsen25,26 to minimize this bias of the chemical potential gradient to diffusivity calculations by constraining solutes to a given position within the membranes and calculating the

10.1021/jp9836731 CCC: $18.00 © 1999 American Chemical Society Published on Web 12/29/1998

386 J. Phys. Chem. B, Vol. 103, No. 2, 1999 diffusion coefficients by a method based on force autocorrelation functions. However, the force autocorrelation function decays rapidly and the long-time behavior, which is known to play an important role in determining diffusivity in a condensed fluid33 (especially for larger solutes for which the limit of overdamped Markovian dynamics may be in reach26), may not be welldetermined because of the substantial noise along the tail of the correlation decay. Nonequilibrium molecular dynamics (NEMD) simulation, which is based on linear response theory,34-36 appears to be suited to overcome the poor response inherent in equilibrium MD methods by imposing a much larger fluctuation via an external force. By measuring the steady-state and instantaneous response to this perturbation, problems associated with the longtime tailing of correlation functions and the heterogeneous nature of lipid bilayers are avoided. This method has grown in popularity over the past decade37 and has recently been applied to calculating gas transport properties in polymers,38 but it has not been used to calculate diffusion coefficients of solute in lipid bilayers and other biomembranes. In this study, theoretical formulas capable of predicting changes of molecular diffusivity in lipid membranes and other interphases with various molecular properties of solute as well as membrane are developed in the framework of dynamic freevolume theory,21 in which contributions of lipid chain overall rotation and local isomerization to free-volume redistribution in lipid membranes are considered explicitly. Following these determinations, equilibrium MD and NEMD simulations are conducted to calculate diffusion coefficients for nonpolar spherical solutes as represented by noble gases into the alkylchain region of lipid bilayers as this region is known by both experiment39,40 and MD simulation25 to be the transport barrier domain across lipid bilayers. The applicability of the NEMD method to lipid membranes is carefully evaluated. The position-, direction-, and size-dependent diffusion coefficients obtained are then compared with solute diffusion behaviors in simple liquids of different molecular structure and analyzed in the framework of free-volume theory. The use of these simple solutes is crucial to reveal basic diffusion mechanisms as the extra complexity associated with nonsphericity, internal rotation, or/and specific interface interactions for a polar solute can be evaluated only if the diffusion patterns for spherical, nonpolar solutes are established. In addition, these calculations have practical implication as noble gases, in particular, argon, are known to be “ideal” anesthesia,41 though their diffusivities in lipid membranes are yet to be determined by direct experiment. Computational Section Molecular Dynamics Simulation. The main features of the MD model used in this study have been described elsewhere.21,42 In brief, the bilayer assembly consisted of 2 × 100 lipid chain molecules and six solutes and was contained in a cubic cell (55, 55, 36 Å) with periodic boundary conditions. The lipid surface area chosen is close to that of dipalmitoylphosphatidylcholine (DPPC) bilayers at 323 K (31 Å2/lipid chain).43 Each lipid molecule had one headgroup, 14 methylene groups, and one methyl group, all modeled as anisotropic united atoms (AUA).42,44 The headgroup had the same size as a carboxylate group and was anchored by an harmonic force (k ) 0.84 kcal/ molÅ2) to a mobile interface plane.42 The headgroups can still move along normal to the bilayer interface and laterally. This treatment is consistent with experimental analyses in which Gaussian functions were found to be appropriate to describe headgroup distributions.45 The same type of force field has been

Xiang TABLE 1: Selected Force-Field Parameters Used in the MD Simulation interaction pair

σ (Å)

 (kcal/mol)

CH2‚‚‚CH2 CH3‚‚‚CH3 X‚‚‚Xa He‚‚‚He Ne‚‚‚Ne Ar‚‚‚Ar Kr‚‚‚Kr Xe‚‚‚Xe Rn‚‚‚Rn

3.53b 3.67c 4.22d 2.63g 2.79g 3.41g 3.67g 3.96g 4.23g

0.159b 0.238b 0.220d 0.012g 0.071g 0.249g 0.336g 0.431g 0.555g

AUA parameterse δ1 ) 0.42 Å, δ2 ) 0.31 Å torsional potential (kcal/mol)f U(φ) ) 1.7704 + 4.3160 cos φ + 1.1314 cos2 φ - 7.3805 cos3 φ + 0.1981 cos4 φ - 0.03539 cos5 φ a X denotes a carboxylate headgroup. b From ref 62. c From the AMBER model.77 d From ref 78. e In the AUA model, the anisotropy is realized by taking the geometric mean position of the valence electrons (whose displacement from the carbon center is denoted as δ1 for methylene and δ2 for methyl) as the origin for intermolecular interactions of the methylene subunits. δ1 was obtained from the refinements of the unit-cell parameters in crystalline polyethylene that optimize the lattice energy,20 and δ2 was determined from δ1 and the ratio δ2/δ1 given in Padilla and Toxvaerd.62 f From ref 42. g From ref 48.

used in MD and statistical mechanics studies of various lipid bilayer properties.46,47 The bond lengths and angles were allowed to vibrate near their equilibrium positions and a torsional potential developed previously was used for torsional motions of lipid chains.42 The Lennard-Jones (L-J) potential was employed for intermolecular interactions and for interactions between segments separated by more than three bonds in the same molecule. The combination rules, ij ) (iijj)1/2 and σij ) (σii + σjj)/2, were used for cross interactions between particles i and j. A cutoff radius of 13 Å was used for all the L-J interactions. The molecular interaction parameters used in the simulation are listed in Table 1. Six nonpolar spherical solutes with size and energy parameters48 equal to those for noble gases (He, Ne, Ar, Kr, Xe, Rn) were investigated. The solutes were assumed to interact with lipid molecules via the same type of the L-J potential as that used for interactions between the lipid chain segments themselves but with different interaction parameters (, σ), which are listed in Table 1. The effects of water-bilayer interfaces on solute motion were described by a continuous medium model in which the solute is subjected to an additional potential V(z) in the interfacial regions49,50

V(z) )

β 1 + (1 - (z ( d)/R)τ

(1)

where d is defined as the average position of the headgroups in the same monolayer and the upper (lower) sign is for the interface at -d (+d). In the bilayer interior and at a distance of >R away from the interfaces, V(z) is assumed to be 0. β is the energy of solvation at the interface. A value of 5.5 kJ/mol was used for β, representing the free energy of transfer of solute from hydrocarbon to water. The steepness of the potential was controlled by R ) 21.2 Å and τ ) -32, similar to the approach employed previously for modeling the energy barrier imposed on methylene groups upon traveling from the air to the water side of a monolayer.49 This interface potential is essentially zero in the bilayer interior, increasing rapidly when the solute crosses

Translational Diffusion in Lipid Membranes

J. Phys. Chem. B, Vol. 103, No. 2, 1999 387

the interface plane. The solutes were subjected to random and frictional forces, both determined by the viscosity of water (0.55 cp), in the regions outside the bilayer interfaces. MD simulations were performed on several Silicon-Graphics Workstations. The equations of motion were integrated by using the “velocity Verlet” algorithm37 and a time step of 2 fs. The system was maintained at 323 K by coupling to a thermal bath.51 The simulation consisted of two phases. In the equilibration phase, the bilayer assembly was equilibrated via a 300 ps run while various bilayer properties were monitored. The presence of solutes were found to induce no significant change of bilayer thickness and lipid chain order parameters. The production phase consisted of several 2-4 ns trajectory runs. The Cartesian coordinates for solutes and lipid molecules were stored every 0.04 and 100 ps, respectively, in files to be used for subsequent calculations. Diffusion Coefficient: Equilibrium MD Calculation. In a plane parallel to the bilayer interface, the bilayer properties are independent of position in the plane. As a result, one can apply the Einstein relation to the calculation of the lateral diffusion coefficient (D⊥) of solute constrained within a given position normal to the bilayer interface

DR )

1 〈|r(t) - r(0)| 〉 D⊥ ) lim 4tf∞ t

(2)

In the above equation, r(t) and r(0) are the position vectors of solute at times t and 0, respectively, in a plane parallel to the bilayer interface. The angled brackets denote an average over all the time origins separated by 20 time steps along the trajectory of solute. A linear segment of the mean-squaredisplacement curve was used to calculate the diffusion coefficient. Mean-square-displacements were also obtained for the movement of unconstrained solutes normal to the bilayer interface to demonstrate the effects of bilayer heterogeneity on solute diffusive motion. Diffusion Coefficient: Nonequilibrium MD Calculation. Various NEMD methods to measure diffusion coefficients which are based on linear response theory34-36 have been described previously.37 In this study, we developed the following NEMD technique to obtained position-dependent diffusion coefficients. The system (bilayer + solute) Hamiltonian H was modified to include an external field on solutes

∑i ni‚ri

H)K+V-F

(3)

where K is the kinetic energy, V the potential energy in the absence of external forces, F the amplitude of the external forces, ni the unit vector of the external force on solute i, and ri the position vector of solute i. The forces on one-half of the solutes acted in the opposite direction to the forces on the other half in order to minimize a translation of the system as a whole. The response of the system at any given time t was represented by the mean velocity VR(t, z) (R ) | and ⊥) of solute as a function of position z normal to the bilayer interface

VR(t,z) )

1

∑j (-1)jVR(t,z)

(4)

n(t, z)

where n(t, z) is the number of solutes whose positions fall within a region centered at z at time t. The diffusion coefficient DR at position z was then given by

(5)

Here k is the Boltzmann constant and T is temperature. Error Estimation. To evaluate errors in the diffusion coefficients calculated from the NEMD method, the total time points along solute trajectories at which the solute resides in a given region centered at z in the bilayer interior were divided into several sets and the diffusion coefficient was calculated from the solute mean velocity in each sample set, from which the mean and error were obtained. Results and Discussion Dynamic Free-Volume Theory for Membrane Diffusivity. In a condensed fluid, solute diffusive motion can be considered as a random-walk process consisting of a large number of uncorrelated elementary displacements. In accord with the Einstein relation and for a continuous distribution of the elementary displacement r, p(r), the diffusion coefficient of the solute can be expressed as52

D) 2

kT 〈V (t,z)〉 F R

2

r p(r) dr ∫0∞D(r)p(r) dr ) 61∫0∞t(r)

(6)

where t(r) is the time a diffusing solute spends to make an elementary displacement r. As a principal hypothesis, the diffusion process is considered as consisting of three sequential steps of opening up a free volume, moving the diffusing solute into this free volume, and filling up the position previously occupied by the solute. The overall time for an elementary diffusive displacement is thus the summation of the characteristic times for the above three steps (τf1, τd, τf2)

t ) τd + 2τf

(7)

where τf [(τf1 + τf2)/2] is the average time for free-volume relaxation. The average rate for free-volume redistribution depends on the probability distribution of free volume V17 or

1 1 ) e-γV/Vf τf τo

(8)

where τo is the relaxation time of the solvent in the limit of low solvent density or that for isolated solvent molecules, Vf is the mean free volume, and γ is a factor in the range of 1/2-1 to account for free-volume overlap.11 This result is consistent with the Williams-Landel-Ferry model in which the temperature dependences of various viscoelastic and dielectric relaxation times of a large number of glass-forming substances are described through their dependences on free volume.53 In general, free-volume redistribution in a condensed solvent may involve several motional modes,54 depending on solvent molecular structure. For lipid bilayers, it is governed primarily by overall rotation and local isomerization motions of lipid molecules. If the characteristic times for these motions in an isolated lipid molecule and the relative contributions of these motions to free-volume redistribution are denoted by (τo,r, τo,i) and (wr, wi), respectively, where wr + wi ) 1, the overall relaxation time for free volume, τf, can be expressed as

τf ) (wrτo,r + wiτo,i)eγV/Vf

(9)

388 J. Phys. Chem. B, Vol. 103, No. 2, 1999

Xiang

Figure 2. Dependence of the mass function, ξ((2ml/md)1/2/Rc)/md1/2, on molecular weight (Da) for noble gases. The molecular weight of the lipid chain, ml, is assumed to be 211 Da: (b) Rc ) 1.0; (+) Rc ) 0.5; and (O) Rc ) 0.1.

Figure 1. Schematic depiction of solute elementary displacement in the acyl-chain region of a lipid bilayer (panel A), and the effects of lipid overall rotation (panel B) and local isomerization (panel C) on solute elementary displacement.

Referring to Figure 1A, an elementary displacement r can occur only if a free volume with a length l g r and a crosssectional area S g that of the diffusing solute (Sd) is created adjacent to the solute. Similar to the treatments by Turnbull and Cohen,55 the elementary displacement r is assumed to be proportional to the free-volume length l

r ) Rl

(10)

where R is a proportionality constant of less than 1. We now consider two limiting cases of wr ) 1 and wi ) 1. If wr ) 1, lipid overall rotation is solely responsible for freevolume redistribution. As shown in Figure 1B, the length l of a free volume created by lipid rotation is dependent on the rotation angle of the lipid or l/τo,r ) cuθL, where uθ is the mean angular velocity for lipid rotation, L is the total length of a lipid chain, and c is a geometric constant depending on position and diffusion direction of solute in the membrane. If the lipid molecule is considered as a classical rotor with a moment of inertia, Is ) 1/3mlL2, where ml is the total mass of the lipid chain, uθ can be expressed as uθ ) (2kT/Is)1/2. Considering that r/τd is the mean translation velocity of solute, or (3kT/md)1/2, we can rewrite eq 6 as

D)

γRx3kT/md 6Vf

le-γSl/Vf

∫0∞∫S∞ d

1 + x2m1/mdeγSl/Vf/Rc

dl dS

(11)

Rx3kT ξ(x2m1/md/Rc) Vf D) 6 γSd m

x

(12)

d

where ξ(x) is

∫0∞ln(1 + e-t/x)dt

D)

γ 12Vfτo,i

∫∫r2e-γSl/V dl dS f

(14)

Associated with lipid isomerization (trans/gauche transitions) is the discrete movement of one or more chain segments. Various collective torsional motions (e.g., jogs and kinks) were previously proposed as one of the primary routes for diffusive motion of water molecules in biological membranes.58 Thus, unlike classical rotation by which the size of a free volume can be varied continuously from zero, a critical free volume, V*, due to chain isomerization for solute diffusion may exist. If V* is approximated as a cylinder, V* ) l*Sd (where l* and Sd are, respectively, the length of the cylinder and the cross-sectional area of the diffusing solute), one can integrate eq 14 from the low limits of (Sd, l*) over the free-volume distribution to obtain the diffusion coefficients of solute

D)

Vf R2 [l*Sd + Vf/2γ]e-2γl*Sd/Vf 2 48τo,i γS

(15)

d

Integrating the above equation, we have

ξ(x) ) 1 - x

For a given lipid membrane, ξ((2ml/md)1/2/Rc) is a function of solute mass md and Rc with both R and c smaller than 1. The mass function, ξ((2ml/md)1/2/Rc)/md1/2, in eq 12 for noble gases at different Rc values is shown in Figure 2. It is noted that the mass function increases nearly linearly with Rc and decreases only slightly with solute mass md, especially at a smaller Rc value. This weak dependence may be due to the fact that lipid overall rotation is the rate-limiting process for diffusive motion of a small solute as lipid rotation is generally on the order of 50-5000 ps,56 which is 2-5 orders of magnitude longer than the time a small solute spends to travel a distance of 1-5 Å. If wi ) 1, lipid isomerization motion is solely responsible for free-volume redistribution. Since the relaxation rate for a barrier-crossing isomerization process is usually much slower than the rate for translation motion of small solutes,20,57 or τd , τo,i, we can write eq 6 as

(13)

The above equation can be solved by a numerical calculation.

As noted, the diffusion coefficient in this limit is independent of solute mass. It is thus concluded that solute mass has a minimal influence on diffusivity regardless of which free-volume redistribution mechanism (overall rotation or local isomerization) dominates. This is supported by our previous MD studies.21 In contrast, as predicted in eqs 12 and 15, solute diffusivity may exhibit markedly different dependences on solute size and membrane mean free volume depending on which free-volume redistribution mechanism (overall rotation or local isomerization) dominates (Vida supra). The model equations derived above

Translational Diffusion in Lipid Membranes

Figure 3. Average velocity for Ar as a function of the external force constant F exerted on the diffusing solute in the ordered peripheral chain region of the lipid bilayer: (b) 〈V|〉 and (O) 〈V⊥〉.

may, therefore, provide a theoretical base on which various diffusion behaviors of solute in lipid bilayers can be rationalized. Since lipid chain ordering in the bilayer interior is known to exhibit a plateau value in the peripheral-chain region and decrease abruptly near the center of the bilayer,42 various diffusion properties in these two distinctive regions will be presented in the following sections. Determination of Diffusion Coefficients by Nonequilibrium MD Method. In this study, a large lipid bilayer with no explicit headgroups and water phase was used as our interest is in the diffusion properties in the alkyl-chain region where solute permeation barrier resides.26,40 Although the chemical nature of headgroups and the degree of hydration can affect membrane structural properties, these factors influence the transport properties of the solute indirectly through their effects on alkylchain order as there is a direct correspondence between acylchain order and surface density, regardless of whether the surface density changes are made by changes of temperature or lipid composition.43,59,60 Selecting a model lipid membrane at a fixed surface density of 31 Å2 per lipid chain, at which the order parameters calculated42 were found to agree closely with those obtained experimentally,61 enabled us to avoid a complete, detailed account of the interfacial interactions, which is extremely time-consuming. The modeling of lipid-chain segments as anisotropic united atoms substantially reduces the number of centers simulated and is widely used in MD simulations of lipid membranes and other thermodynamic systems. This model has been shown to reproduce the experimental diffusivity of small solutes in hydrocarbon solvents and amorphous polymers.20,62 The above arguments, thus, support the model lipid bilayer developed in this study as an adequate model of real lipid membranes. A successful NEMD method must comply with at least three criterions. First, as the NEMD method is based on linear response theory, it works only if the response as represented here by the mean velocity of solute varies linearly with the external force. Second, the applied perturbation should have a minimal impact on bilayer thermodynamic properties such as temperature and bilayer thickness. Third, the diffusion coefficients calculated from the NEMD method should be consistent with those calculated from other well-established MD methods and with available experimental results. The range of the external force F over which the response of solute mean velocity is linear is evaluated for noble gases varying in size. Representative plots for Ar diffusing in the ordered peripheral region (the | and ⊥ components) are shown in Figure 3. As noted, the response increases approximately

J. Phys. Chem. B, Vol. 103, No. 2, 1999 389

Figure 4. Average kinetic energy for Ar (O) and for the whole system (b) as a function of the external force constant F exerted on the solutes in the lipid bilayer.

Figure 5. Bilayer thickness as a function of the external force constant F.

linearly in the plane normal to the bilayer interface, but this linearity appears to break down at an external force of g8 kJ/ mol/Å normal to the bilayer interface. Similar behaviors are also found for other noble gases. The presence of an external force will accelerate solute motion, leading to an elevation of solute kinetic velocity. The dissipation of this kinetic energy to surrounding lipid molecules may increase the overall system temperature and affect other thermodynamic properties such as bilayer thickness, but this increase of system temperature may be curtailed by the use of the loose-coupling thermostat25 as in the present MD simulation. Figure 4 presents representative plots of the kinetic energies (K) for Ar and the whole bilayer system as a function of the external force F. The dependence of bilayer thickness on F is shown in Figure 5. As noted, the system temperature remains approximately constant ((0.5 K) for F e 2 kJ/mol/Å and shows a small increase (ca. +3 K) at F ) 5 kJ/mol/Å while the solute kinetic energy increases more substantially than the overall system temperature (+29 K at F ) 5 kJ/mol/Å). This can be attributed to the larger number of degrees of freedom for the bilayer system as a whole. The external force has no significant influence on bilayer thickness, suggesting that the effects of the external force on bilayer structure are small. A smaller solute is subjected to a larger increase of the kinetic velocity than a larger solute (data not shown) due to the larger acceleration of a lighter solute at a given external force. To deliver a maximum response while leaving the system in the linear-response region, the external forces of 1-5 kJ/mol/Å are used for the NEMD calculations of diffusion coefficients for noble gases. Figure 6 shows natural log-log plots of mean-squaredisplacements (MSD) for Ar along (〈r|2〉) and perpendicular to

390 J. Phys. Chem. B, Vol. 103, No. 2, 1999

Xiang

TABLE 2: Diffusion Coefficients (cm2/s) in the Bilayer Interior Obtained from NEMD Simulation peripheral region

center of bilayer

diffusant

D|

D⊥

D|

D⊥

He Ne Ar Kr Xe Rn

(1.8 ( 0.4) × 10-4 (6.9 ( 0.9) × 10-5 (5.9 ( 1.2) × 10-5 (4.7 ( 2.0) × 10-5 (4.2 ( 1.3) × 10-5 (3.2 ( 0.9) × 10-5

(1.5 ( 0.1) × 10-4 (8.5 ( 1.2) × 10-5 (5.0 ( 0.6) × 10-5 (2.9 ( 0.7) × 10-5 (2.4 ( 1.0) × 10-5 (1.2 ( 0.3) × 10-5

(1.8 ( 0.5) × 10-4 (1.5 ( 0.2) × 10-4 (9.6 ( 1.2) × 10-5 (6.9 ( 1.5) × 10-5 (7.0 ( 2.3) × 10-5 (4.0 ( 1.6) × 10-5

(1.7 ( 0.4) × 10-4 (1.0 ( 0.1) × 10-4 (6.0 ( 0.5) × 10-4 (4.9 ( 1.8) × 10-5 (3.7 ( 1.1) × 10-5 (3.2 ( 1.3) × 10-5

Figure 6. Natural log-log plots of mean-square-displacements (Å2) vs time (ps) for Ar in the bilayer interior: (×) 〈r|2〉; and (+) 〈r⊥2〉.

(〈r⊥2〉) the bilayer interface versus time. In these equilibrium MD simulations, 〈r⊥2〉 is obtained when the solute is constrained to the ordered chain region in the bilayer while 〈r|2〉 is obtained for unconstrained solute. The slope in the log-log plot of 〈r⊥2〉 vs t is 1.1, suggesting that the lateral diffusion for Ar approximately follows the Einstein relation. In contrast, 〈r|2〉 is nonlinear with a slope substantially smaller than 1, especially at long times. This leveling of 〈r|2〉 may result from the fact that the membrane interior is not homogeneous normal to the bilayer interface. Lipid chains are more ordered in the peripheral region than near the center of the bilayer, and noble gases also suffer unfavorable hydrophobic interaction when they move from the acyl-chain region toward the hydrophilic interface, which is modeled by eq 1. As a result, these solutes tend to move into the more disordered region near the center of the bilayer.28-30 Consequently, the mean-squared displacement of solute normal to the bilayer interface, 〈r|2〉, is nonlinear and levels off at long times, as shown in Figure 6. Similar longtime behaviors are also found in a simple Brownian dynamics simulation of solute enclosed within membrane walls.32 This violates the assumption for the Einstein relation, namely, the mean-squared displacement increases linearly with time in the limit of t f ∞ as usually found in a homogeneous medium. It is then clear that the equilibrium MD simulation on the basis of the Einstein relation is unable to obtained accurate D|. The lateral diffusion coefficients (D⊥) for Ar obtained from the MSD, as shown in Figure 6, is 4.3 × 10-5 cm2/s in the ordered chain region, which compares favorably with D⊥ obtained from NEMD simulation (5.0 × 10-5 cm2/s; cf., Table 2). Smaller diffusion coefficients on the order of (12-40) × 10-6 cm2/s are obtained for the largest noble gas, Rn. In an MD study of diffusion of a larger nonspherical solute, benzene, in DMPC bilayers by Bassolino-Klimas et al.,23 the diffusion coefficient was found to be (1-3) × 10-6 cm2/s in the alkylchain region. In an earlier 13C NMR spin-relaxation experiment by Dix et al.,63 the diffusion coefficient for an even larger solute, di-tert-butyl nitroxide, in DPPC bilayers was found to be (0.52.0) × 10-6 cm2/s depending on position in the bilayer interior. Although a direct comparison of our NEMD data with these previous results is not warranted considering the differences in

experimental conditions and solute size/shape, the lower diffusion coefficients for larger solutes used in these previous studies are consistent with the general propensity of diffusivity decrease with increasing solute size. It is also interesting to compare the diffusion coefficients in the membrane to those in a hydrocarbon solvent as the latter has similar chemical composition and is commonly use as a model for biomembranes. Evans et al.64 conducted experimental studies on diffusivity of noble gases and other spherical solutes in hydrocarbon solvents using the Taylor dispersion technique. The diffusion coefficients for Ar, Kr, and Xe in tetradecane at 298 K, which has a viscosity of 2.1 cp, are 3.4 × 10-5, 2.4 × 10-5, and 1.7 × 10-5 cm2/s, respectively. These results are close to lateral diffusion coefficients in the ordered peripheral region of the lipid bilayer but somewhat smaller than diffusion coefficients near the center of the bilayer. The latter may be attributed to a smaller free-volume fraction in tetradecane (0.4365) than in the middle of the bilayer (0.4922) where more chain-ending methyl groups reside. The closeness of the diffusion coefficients in the membrane and in tetradecane appears to contradict with over 100-fold higher viscosity found for lateral diffusion of lipids and proteins in lipid membranes.2,66,67 However, a recent comparative MD study of chain segmental motions in lipid bilayers and hydrocarbon solvents found no significant difference in the corresponding reorientation times.68 It was concluded that the effective viscosity in the alkyl-chain region of a DPPC bilayer at 323 K is very similar to that (∼2 cp) of a neat alkane of approximately the same chain length and the apparently higher viscosity is more closely related to molecular interactions on the hydrated interface rather than the viscosity in the bilayer interior. This is further supported by NMR relaxation and fluorescence studies of correlation times of chain segments and translational diffusivity of small molecules in lipid bilayers.3,69 Together, the above results indicate that one can choose an optimal value for the external force that generates a linear response of solute mean velocity, minimizes the perturbation to bilayer properties, and gives solute diffusion coefficients consistent with other simulation and experimental methods. Dependence of Diffusion Coefficients on Bilayer Structure. Table 2 presents the diffusion coefficients (D| and D⊥) for noble gases in the ordered peripheral-chain region and near the center of the bilayer. It is noted that the diffusion coefficient for the smallest solute, He, is only weakly dependent on its position in the bilayer interior (1.5-1.8 × 10-4 cm2/s in the ordered peripheral region and 1.7-1.8 × 10-4 cm2/s in the center of the bilayer), whereas the diffusion coefficient for the largest solute, Rn, in the center of the bilayer (3.2-4.0 × 10-5 cm2/s) is 1.3-2.6 times larger than near the bilayer interface (1.2-3.2 × 10-5 cm2/s). Large diffusivity differences were also found for diffusion of other larger solutes such as benzene and ditert-butyl nitroxide in phospholipid bilayers in previous MD and NMR spin-relaxation studies.23,63 The diffusion coefficient for He is only weakly dependent on diffusion direction, while a significant disparity in D| and D⊥ is observed for Rn. In

Translational Diffusion in Lipid Membranes

J. Phys. Chem. B, Vol. 103, No. 2, 1999 391 TABLE 3: The Diffusion Coefficients (cm2/s) Estimated from the Distribution of Free-Flight Length According to Eq 17 peripheral region

center of bilayer

diffusant

D|

D⊥

D|

D⊥

He Ne Ar Kr Xe Rn

4.4 × 1.6 × 10-4 8.4 × 10-5 6.6 × 10-5 5.9 × 10-5 4.1 × 10-5

3.4 × 1.2 × 10-4 7.0 × 10-5 4.3 × 10-5 3.6 × 10-5 2.1 × 10-5

4.0 × 1.4 × 10-4 9.1 × 10-5 6.9 × 10-5 6.0 × 10-5 3.3 × 10-5

4.1 × 10-4 1.5 × 10-4 9.5 × 10-5 6.6 × 10-5 4.3 × 10-5 3.2 × 10-5

10-4

10-4

10-4

Figure 7. Probability distribution of the free-flight length for the motion of Ar in a plane normal to the bilayer interface and located in the ordered peripheral chain region.

general, D| is larger than D⊥, reflecting the strong effects of chain ordering on solute diffusion. Similar behaviors are also found in liquid crystals which also exhibit long-range order.5 According to eqs 12 and 15, the size of the free volume in the membrane interior has a strong effect on solute diffusivity. The free-volume distribution in lipid bilayers has been investigated by Xiang.22 It was found that the free-volume fraction is larger near the center of the bilayer than in the ordered peripheral-chain region and, on average, the free volume is elongated and prefers to align normal to the bilayer interface. These findings parallel with the present MD results of a larger diffusion coefficient near the center of bilayer and normal to the bilayer interface, though a quantitative measure of freevolume anisotropy is lacking. Here, we develop a quantitative method of describing free volume that takes into consideration free-volume anisotropy, nonsphericity, connectivity (overlapping), and the perturbation of solute to the local organization of lipid molecules, namely, the probability distribution of free-flight length, lR, of solute. It is defined as

lR ) tf,Rud,R

(16)

where tf,R is the free-flight time and ud,R the velocity of solute along the R direction (⊥ or |). The free-flight time, tf,R, is determined here by the time interval during which the direction of solute motion along the R direction remains unchanged or the velocity product between an initial time to and a later time t, ud,R(t) × ud,R(to), is positive. Figure 7 shows a representative plot of the probability distribution, P(l⊥), for Ar located within the ordered peripheral chain region. As noted, the distribution reaches a maximum at a small l⊥ (∼1 Å) and declines rapidly as l⊥ increases. Similar behaviors are also found for diffusion normal to the bilayer interface, in other regions of the bilayer, and for other noble gases. For Ar, the mean free-flight lengths, l⊥ and l|, are 1.2 and 1.5 Å, respectively, in the ordered peripheral region and 1.6 and 1.7 Å, respectively, near the center of the bilayer. Thus, a solute near the center of the bilayer can travel freely for a longer distance without reversing its direction and sees a roughly isotropic environment while a solute in the ordered peripheral region sees a more anisotropic neighborhood with a significantly shorter free-flight length in a plane normal to the bilayer interface. Given the probability distribution for the free-flight length, P(lR), and the assumption that the elementary displacement of solute is close to the free-flight length, or r ≈ l, one can estimate

Figure 8. Average R parameter as a function of solute diameter σ (Å).

the diffusion coefficient on the basis of the Einstein relation (cf., eqs 6), namely,

DR )

∫0∞lRud,RP(lR)dlR

1 2

(17)

The diffusion coefficients so obtained are presented in Table 3. Comparing results in Tables 2 and 3, we note that the diffusion coefficients calculated from eq 17 (D(eq 17)) reproduce the general trends of lower diffusion coefficients in a more ordered chain region and in a plane normal to the bilayer interface as found from the exact NEMD calculations (D(exact)). The deviation of D(eq 17) from the exact NEMD calculation can be attributed to the difference of the elementary displacement r and the free-flight length l or the assumption r ≈ l should be replaced by r ) Rl as defined in eq 10. The R parameter can then be obtained from the ratio of the diffusion coefficients D(exact) and D(eq 17)

R)

[

]

D(exact) D(eq 17)

1/2

(18)

R is close to 1 if D(exact) ≈ D(eq 17). Figure 8 shows the dependence of the R parameter averaged over different positions and diffusion directions on solute size. It is noted that the diffusion coefficients obtained from eq 17 are close to the exact results within 20% for most noble gases, though R is generally smaller for smaller noble gases, especially for He. This may be attributed to a greater degree of rattling motions for smaller noble gases which have higher kinetic velocities. These rattling motions tend to reduce the elementary displacements r and thereby the diffusion coefficients of solute. Dependence of Diffusion Coefficients on Solute Size. Molecular size is one of the primary determinants of molecular diffusivity and structure-activity relations in many biological processes.70,71 The dependence of diffusivity on solute diameter (σ) is known to be highly variable depending on the molecular structure of the medium and can be evaluated, for small

392 J. Phys. Chem. B, Vol. 103, No. 2, 1999

Xiang TABLE 5: The n Parameter Predicted from Dynamic Free-Volume Theory free-volume redistribution mechanism translation overall rotation

Vf/γ (Å3)

TABLE 4: The n Parameter Defined in Eq 19

a

diffusion component/medium

n

D| (peripheral region) D⊥ (peripheral region) D| (bilayer center) D⊥ (bilayer center) D (CCl4)a D (hexane)b D (tetradecane)b

2.9 4.7 2.2 2.8 2.1 3.0 4.5

1.0 0.5 0.1

solutes,72 by the following power law

Do σn

2.0 2.0 2.0 5.0 5.0 5.0 10.0 10.0 10.0 30.0 30.0 30.0

0.0 0.1 0.2 0.5 0.1 0.2 0.5 0.1 0.2 0.5 0.1 0.2 0.5

n 2.35b 2.42 2.23 2.05 4.00 4.85 6.28 11.2 4.19 4.60 6.28 4.06 4.19 4.28 4.01 4.03 4.14

a Independent of this parameter. b Predicted for diffusion of small solutes in carbon tetrachloride according to a formula derived in our previous study.21

From ref 79. b From ref 64.

D)

Rc

a

local isomerization

Figure 9. Log-log plots of the diffusion coefficients for noble gases in the ordered peripheral region ((b) D⊥ and (O) D|) as a function of solute diameter σ. The diffusion coefficients for three noble gases (Ar, Kr, Xe) in tetradecane64 (×) and small gaseous molecules in carbon tetrachloride79 (+) are also plotted for comparison.

l* (Å)

(19)

For example, n is 2.1 for solute diffusion in simple liquids73,74 and 3.3-11.4 in natural rubber and polymethylacrylatecan.75,76 Therefore, it is important to explore how this size dependence varies with solute position and diffusion direction in lipid bilayers where chain organization is anisotropic and varies with position. Figure 9 shows log-log plots of the diffusion coefficients (D| and D⊥) for noble gases in the ordered peripheral-chain region as a function of solute diameter, σ. For comparison, the diffusion coefficients for three noble gases (Ar, Kr, Xe) in tetradecane and small gaseous molecules in carbon tetrachloride are also plotted in Figure 9. The data in Figure 9 are fitted according to eq 19. The resulting n values are reported in Table 4 along with those for diffusion near the center of the bilayer and in hexane. It is noted that the n parameter is highly variable depending on position and diffusion direction in the bilayer interior. For D| and for diffusion near the center of the bilayer, n varies in a narrow range of 2.2-2.9, whereas for D⊥ in the ordered peripheral region, a much stronger size-dependence, n ) 4.7, is obtained. This latter result is closer to that (n ) 4.0) obtained for lateral diffusion (D⊥) of lipophilic fluorescent probes in DMPC membranes,72 which, because of their large molecular weights (223-856 Da), should reside, at least partially, within the ordered plateau chain region. Rigorously speaking, the diffusion models predicted from our dynamic free-volume theory (cf., eqs 12 and 15) are more complex than the simple power law in eq 19 as the ξ function in eq 12 depends on solute mass which varies with solute size in a complex manner and there is an exponential factor of solute size in eq 15. To compare the size dependences of solute

diffusivity obtained from the present MD simulations and previous experiments with the dynamic free-volume theory, the size dependences predicted from the free-volume theory are plotted in a format of log D vs -log σ and fitted according to eq 19 (data not shown). The slopes obtained are considered to be the n parameters and are listed in Table 5. If solvent overall rotation is responsible for free-volume redistribution, the n parameter stays in a narrow range of 2.1-2.4. In contrast, if solvent chain isomerization is the only mechanism for freevolume redistribution, the size dependence is substantially stronger (n g 4.0) and varies with the mean free volume Vf and the minimum free-flight-length, l* (cf., eq 15). From Figure 7, it is noted that there is a nonzero probability for free-flight length, lf, at lf > 0.5 Å, though lf approaches zero at a smaller value. Thus, l* should be less than 0.5 Å. In our previous MD study on free-volume distributions in lipid bilayers,22 we found that the mean free volume varies in the range of 3-10 Å3 depending on position in the bilayer interior. Thus, the n values at Vf/γ ) 2, 5, 10, and 30 Å3 are presented in Table 5. As noted, the n value increases with decreasing Vf/γ or when the medium is more densely packed. The power law in eq 19 works well for n values not considerably larger than 4.0. For large l* and small Vf/γ where n . 4.0, the log D - log σ plot becomes increasingly nonlinear (data not shown) because of the greater contribution of the exponential factor in eq 15 to solute diffusivity. The strong correlation between the n parameter and the molecular structure of solvent, as shown in Table 5, is consistent with the size dependences of molecular diffusion in carbon tetrachloride, hexane, and tetradecane (cf., Table 4). With the increase of chain length from hexane to tetradecane, chain isomerization becomes more important for free-volume redistribution leading to a stronger size dependence for diffusion in tetradecane (n ) 4.5) in comparison with that in hexane (n ) 3.0). An even weaker size dependence (n ) 2.1) is found in carbon tetrachloride where chain isomerization is absent. This is close to the n value (2.35) predicted from the dynamic freevolume theory if solvent translation motion is responsible for free-volume redistribution.21 As illustrated in Figure 1C, lipid local isomerization can open up a large free volume for lateral diffusion (D⊥) in the ordered peripheral chain region while lipid overall rotation may be more effective to create elongated free volume for solute diffusion

Translational Diffusion in Lipid Membranes normal to the bilayer interface (D|). The region near the center of the bilayer contains more chain-ending methyl groups resembling a short-chain hydrocarbon solvent such as hexane. Thus, based on our dynamic free-volume theory, lateral diffusion (D⊥) in the ordered peripheral chain is expected to exhibit a stronger size dependence than transverse diffusion (D|) and diffusion near the center of the bilayer. These predictions are in close agreement with the present MD results in Table 4, suggesting that molecular diffusivity in lipid bilayers conforms to dynamic free-volume theory. Conclusions In summary, a nonequilibrium MD method is developed to calculate various diffusion coefficient components for noble gases as a function of position in the model membrane. With a proper choice of the external force, this method is found to give lateral diffusion coefficients in close agreement with the equilibrium MD method. It avoids problems associated with the long-time behaviors of various correlation functions and the heterogeneous nature of lipid chain organization as encountered in an equilibrium MD simulation. The diffusion coefficients obtained are found to depend strongly on position and diffusion direction in the membrane and follow closely those obtained from free-flight-length distributions, a newly developed measure of free-volume anisotropy, on the basis of the Einstein relation, except for smaller solutes for which rattling motions within cavities become important. Furthermore, the diffusion coefficients vary with solute size with the degree of changes depending strongly on solute position and diffusion direction in the membrane. To explore molecular mechanisms underlying these observations, theoretical formulas for solute diffusion in lipid bilayers where overall rotation and local isomerization of lipid chain molecules are responsible for free-volume redistribution are derived in the framework of dynamic free-volume theory. This theory predicts that solute diffusion exhibits variable dependences on solute size depending on relative contributions of overall rotation and local isomerization of lipid molecules to free-volume redistribution. The observed strong (weak) dependences of lateral (transverse) diffusion coefficients on solute size can then be explained by the greater contribution of lipid chain isomerization (overall rotation) to free-volume redistribution. Acknowledgment. This work was supported by a grant from the National Institute of Health (Grant No. RO1 GM5134701A2). The authors also thank the College of Pharmacy and Dr. J. N. Herron of the University of Utah for the use of their computers in some of the simulations. References and Notes (1) Cone, R. A. Nature New Biol. 1972, 236, 39. (2) Peters, R.; Cherry, R. J. Proc. Natl. Acad. Sci. U.S.A. 1982, 79, 4317. (3) Lakowicz, J. R. Spectroscopy in Biochemistry; Chemical Rubber: Boca Raton, FL, 1981; Vol. 1. (4) Vrentas, J. S.; Vrentas, C. M. J. Polym. Sci., Polym. Lett. 1990, 28, 379. (5) Moscicki, J. K.; Shin, Y.-K.; Freed, J. H. J. Chem. Phys. 1993, 99, 634. (6) Stein, W. D. Transport and Diffusion Across Cell Membranes; Academic Press: Orlando, FL, 1986. (7) Vaz, W. L. C.; Clegg, R. M.; Hallmann, D. Biochemistry 1985, 24, 781. (8) Almeida, P. F. F.; Vaz, W. L. C.; Thompson, T. E. Biochemistry 1992, 31, 6739.

J. Phys. Chem. B, Vol. 103, No. 2, 1999 393 (9) Xiang, T.-X.; Anderson, B. D. Biophys. J. 1997, 72, 223. (10) Xiang, T.-X.; Anderson, B. D. Biophys. J. 1998, 75, 2658. (11) Cohen, M. H.; Turnbull, D. J. Chem. Phys. 1959, 31, 1164. (12) Macedo, P. B.; Litovitz, T. A. J. Chem. Phys. 1965, 42, 245. (13) Vrentas, J. S.; Duda, J. L. J. Polm. Sci. 1977, 15, 403. (14) Fujita, H. Polym. J. 1991, 23, 1499. (15) Zielinski, J. M.; Duda, J. L. AIChE J. 1992, 38, 405. (16) Chapman, S.; Cowling, T. G. The mathematical theory of nonuniform gases; Cambridge University Press: Cambridge, 1970. (17) Bueche, F. Physical properties of polymers; Interscience: New York, 1962. (18) Vrentas, J. S.; Jarzebski, C. M.; Duda, J. L. AIChE J. 1975, 21, 894. (19) Wu, J. C.; Peppas, N. A. J. Polym. Sci. Part B: Polym. Phys. 1993, 31, 1503. (20) Krishna Pant, P. V.; Boyd, R. H. Macromolecules 1993, 26, 679. (21) Xiang, T.-X. J. Chem. Phys. 1998, 109, 7876. (22) Xiang, T.-X. Biophys. J. 1993, 65, 1108. (23) Bassolino-Klimas, D.; Alper, H. E.; Stouch, T. R. Biochemistry 1993, 32, 12624. (24) Bassolino-Klimas, D.; Alper, H. E.; Stouch, T. R. J. Am. Chem. Soc. 1995, 117, 4118. (25) Marrink, S.-J.; Berendsen, H. J. C. J. Phys. Chem. 1994, 98, 4155. (26) Marrink, S. J.; Berendsen, H. J. C. J. Phys. Chem. 1996, 100, 16729. (27) Jin, B.; Hopfinger, A. J. Pharm. Res. 1996, 13, 1786. (28) Marqusee, J. A.; Dill, K. A. J. Chem. Phys. 1986, 85, 434. (29) Xiang, T.-X.; Anderson, B. D. Biophys. J. 1994, 66, 561. (30) Pohorille, A.; Cieplak, P.; Wilson, M. A. Chem. Phys. 1996, 204337, 337. (31) Vertenstein, M.; Ronis, D. J. Chem. Phys. 1987, 87, 5457. (32) Huertas, M. L.; Cruz, V.; Cascales, J. J. L.; Acuna, A. U.; de la Torre, J. G. Biophys. J. 1996, 71, 1428. (33) Tildesley, D. J.; Madden, P. A. Mol. Phys. 1983, 48, 129. (34) Kubo, R. J. Phys. Soc. Jpn. 1957, 12, 570. (35) Kubo, R. Rep. Prog. Phys. 1966, 29, 255. (36) Zwanzig, R. Annu. ReV. Phys. Chem. 1965, 16, 67. (37) Allen, M. P.; Tildesley, D. J. Computer simulation of liquids; Clarendon Press: Oxford, 1992. (38) Muller-Plathe, F.; Rogers, S. C.; Gunsteren, W. F. v. J. Chem. Phys. 1993, 98, 9895. (39) Katz, Y. Biochim. Biophys. Acta 1981, 647, 119. (40) Xiang, T.-X.; Anderson, B. D. J. Pharm. Sci. 1994, 83, 1511. (41) Kennedy, R. R.; Stokes, J. W.; Downing, P. Anaesth. IntensiVe Care 1992, 20, 66. (42) Xiang, T.-X.; Anderson, B. D. J. Chem. Phys. 1995, 103, 8666. (43) Nagle, J. F. Biophys. J. 1993, 64, 1476. (44) Toxvaerd, S. J. Chem. Phys. 1990, 93, 4290. (45) Wiener, M. C.; White, S. H. Biophys. J. 1992, 61, 434. (46) Edholm, O.; Nyberg, A. M. Biophys. J. 1992, 63, 1081. (47) Jacobs, R. E.; Hudson, B. S.; Andersen, H. C. Biochemistry 1977, 16, 4349. (48) Pierotti, R. A. Chem. ReV. 1976, 76, 717. (49) Karaborni, S.; Toxvaerd, S.; Olsen, O. H. J. Phys. Chem. 1992, 96, 4965. (50) Peters, G. H.; Toxvaerd, S.; Sverndsen, A.; Olsen, O. H. J. Chem. Phys. 1994, 100, 5996. (51) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (52) Frenkel, J. Kinetic theory of liquids; Dover: New York, 1955. (53) Williams, M. L.; Landel, R. F.; Ferry, J. D. J. Am. Chem. Soc. 1955, 77, 3701. (54) Sceats, M. G.; Dawes, J. M. J. Chem. Phys. 1985, 83, 1298. (55) Turnbull, D.; Cohen, M. H. J. Chem. Phys. 1961, 34, 120. (56) Pastor, R. W.; Venable, R. M. J. Chem. Phys. 1988, 89, 1128. (57) Choi, K.; Jo, W. H. Macromolecules 1995, 28, 8598. (58) Trauble, H. J. Membr. Biol. 1971, 4, 193. (59) DeYoung, L. R.; Dill, K. A. Biochemistry 1988, 27, 5281. (60) Boden, N.; Jones, S. A.; Sixl, F. Biochemistry 1991, 30, 2146. (61) Seelig, A.; Seelig, J. Biochemistry 1974, 13, 4839. (62) Padilla, P.; Toxvaerd, S. J. Phys. Chem. 1991, 94, 5650. (63) Dix, J. A.; Kivelson, D.; Diamond, J. A. J. Membr. Biol. 1978, 40, 315. (64) Evans, D. F.; Tominaga, T.; Chan, C. J. Solution Chem. 1979, 8, 461. (65) Bondi, A. J. Phys. Chem. 1954, 58, 929. (66) Saffman, P. G.; Delbruck, M. Proc. Natl. Acad. Sci. U.S.A. 1975, 72, 3111. (67) Hughes, B. D.; Pailthope, B. A.; White, L. R.; Sawyer, W. H. Biophys. J. 1982, 37, 673. (68) Venable, R. M.; Zhang, Y.; Hardy, B. J.; Pastor, R. W. Science 1993, 262, 223.

394 J. Phys. Chem. B, Vol. 103, No. 2, 1999 (69) Brown, M. F. J. Chem. Phys. 1984, 80, 2808. (70) Xiang, T.-X.; Anderson, B. D. J. Membr. Biol. 1994, 140, 111. (71) Potts, R. O.; Guy, R. H. Pharm. Res. 1995, 12, 1628. (72) Johnson, M. E.; Berk, D. A.; Blankschtein, D. B.; Golan, D. E.; Jain, R. K.; Langer, R. S. Biophys. J. 1996, 71, 2656. (73) Powell, R. J.; Hildebrand, J. H. J. Chem. Phys. 1971, 55, 4715. (74) Hildebrand, J. H.; Lamoreaux, R. H. Proc. Natl. Acad. Sci. U.S.A. 1974, 71, 3321.

Xiang (75) Lieb, W. R.; Stein, W. D. Nature 1969, 224, 240. (76) Coughlin, C. S.; Mauritz, K. A.; Storey, R. F. Macromolecules 1990, 23, 3187. (77) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, J.; Weiner, P. J. Am. Chem. Soc. 1984, 106, 765. (78) van der Ploeg, P.; Berendsen, H. J. C. J. Chem. Phys. 1982, 76, 3271. (79) Nakanishi, K.; Voigt, E. M.; Hildebrand, J. H. J. Chem. Phys. 1964, 42, 1860.