## Abstract

Recently measured water permeability through bilayers of different lipids is most strongly correlated with the area per lipid *A* rather than with other structural quantities such as the thickness. This paper presents a simple three-layer theory that incorporates the area dependence in a physically realistic way and also includes the thickness as a secondary modulating parameter. The theory also includes the well-known strong correlation of permeability upon the partition coefficients of general solutes in hydrocarbon environments (Overton's rule). Two mathematical treatments of the theory are given; one model uses discrete chemical kinetics and one model uses the Nernst-Planck continuum equation. The theory is fit to the recent experiments on water permeability in the accompanying paper.

## Introduction

A highly favored theory of passive permeability through lipid bilayers and biomembranes uses the solubility- diffusion (SD) model. This supposes that, for the purpose of understanding permeability *P*, the bilayer can be modeled as a single layer of hydrocarbon of thickness *d _{C}*. This leads directly to the well-known formula,(1)where

*K*is the partition coefficient of the solute into the hydrocarbon core of the lipid bilayer and

*D*is the coefficient of diffusion of the solute in the same environment. For small solutes,

_{C}*D*is often (but not always, see Lieb and Stein, 1986) assumed to be weakly dependent upon solute. The strong dependence of

_{C}*P*, varying over nearly six orders of magnitude for different solutes for a given lipid bilayer (often egg lecithin), is interpreted as the dependence of

*K*on the solute (Overton's rule). This conclusion is reinforced by the fact that the measured partition coefficients of solutes into bulk hydrocarbon correlate fairly well with permeabilities measured over the same six orders of magnitude (Walter and Gutknecht, 1986; Finkelstein, 1987). This is a major result that any theory must account for. Nevertheless, the fact that the single layer SD theory easily accommodates Overton's rule does not prove that it is the correct model. One concern about the SD theory is that the value of

*d*obtained from calculating

_{C}*KD*exceeds 10 nm for egg lecithin bilayers (Finkelstein, 1987), but the structural thickness of the hydrocarbon core for that lipid bilayer is only 2.7 nm (Nagle and Tristram-Nagle, 2000). The theory in this paper removes this discrepancy.

_{C}/PAnother concern with the single layer SD theory regards how to incorporate the dependence of *P* for a given solute with the area per lipid *A* for different bilayers. It may be noted first that correlation with *A* is different than correlation with inverse thickness 1/*d _{C}* because, even though the product

*Ad*=

_{C}*V*is the volume of the hydrocarbon region,

_{C}*V*is considerably different for lipids with different numbers of carbons in the hydrocarbon chains. Indeed, there is no apparent experimental correlation of the water permeability with

_{C}*d*whereas there is a strong, though not perfect, correlation with

_{C}*A*(Mathai et al., 2007). The more relevant structural quantity for discussing the SD theory is the volume per methylene group

*V*in the hydrocarbon core. The partition coefficient

_{CH2}*K*should increase monotonically with

*V*, as in the “free volume” theory, so

_{CH2}*V*should be the first order structural quantity to correlate with

_{CH2}*K*. If there were a strong correlation of

*V*with

_{CH2}*A*, then the

*A*dependence of

*P*could be easily understood as a

*K*dependence within the single layer SD theory. Contrarily, all the lipid bilayers employed in the recent experimental study of water permeability have essentially the same value of

*V*. It may be emphasized that the structural values of

_{CH2}*V*were obtained from straightforward measurements of the total lipid volume that are highly accurate (Nagle and Tristram-Nagle, 2000; Koenig and Gawrisch, 2005; Greenwood et al., 2006; Heerklotz and Tsamaloukas, 2006). The largest uncertainty was how much to subtract for the volume of the headgroup, but that number should be the same for all phosphatidylcholine lipids in their fully hydrated bilayers, so any discrepancy only changes

_{CH2}*V*by essentially the same amount for all bilayers. This volumetric result precludes a simple reconciliation of the single layer SD theory with experiment, although a more complex reconciliation based on a lattice model has been proposed (DeYoung and Dill, 1990; Xiang and Anderson, 1995).

_{CH2}This paper therefore goes beyond the single layer SD theory and considers three layer theories such as have been considered by Zwolinski et al. (1949) and Diamond and Katz (1974). The main new idea is that the area dependence is quite naturally included in the interfacial headgroup layers rather than in the fluid hydrocarbon core layer. This theory will be implemented with close comparison to recent water permeability measurements that were made on five pure lipid bilayers, all with the same phosphatidylcholine headgroup and all at the same temperature and all with structures recently determined by the same high resolution X-ray method for fully hydrated lipid bilayers.

## THEORY AND RESULTS

### I. Three Layer Theory

Before deriving detailed equations from mathematical models, let us develop the major ideas in a phenomenological and intuitive manner. The underlying theory assumes three layers, an inner hydrocarbon core, as in the single layer SD theory, and two interfacial headgroup layers. Let us define the permeabilities through each part separately. Let *P _{C}* be the permeability that would apply just within the hydrocarbon core and let

*P*be the permeability through the interfacial region and including, importantly, transfer into the hydrocarbon core. Then, as is well known (Zwolinski et al., 1949; Diamond and Katz, 1974) and as will be shown in detail in the following two subsections, the permeability

_{H}*P*of the three layer composite model is given by(2)which is just the formula for addition of resistances in series where each of the three separate resistances is proportional to its inverse permeability. A recent experimental study suggested that the headgroup regions and the hydrocarbon region each offer independent and additive resistance to permeation (Krylov et al., 2001).

The most important aspect of our model is the functional form for *P _{H}*. As suggested by Fig. 1, we suppose that the headgroups sterically block the entrance of water into the hydrocarbon region. We therefore propose a structural factor of (

*A*−

*A*)/

_{0}*A*in

*P*to account for the fraction of the total area A that is not blocked. The parameter

_{H}*A*is the area at which the headgroups are packed so tightly that the permeability becomes negligible. Xiang and Anderson (1997) have measured the permeability of acetic acid in the gel phase of DPPC to be 482 times less than in the fluid phase, so a first approximation for

_{0}*A*is the area of the gel phase. The theory will not attempt to account for gel phase permeability, which appears to be qualitatively different from fluid phase permeability (Xiang et al., 1998). For phosphatidylcholine lipids the gel phase area is ∼48 Å

_{0}^{2}and the chains are tilted (Tristram-Nagle et al., 2002). As pointed out by McIntosh (1980), tilting shows that the phosphatidylcholine (PC) headgroups are tightly jammed together in the gel phase. Our permeability model essentially assumes that the headgroups comprise a partial barrier for entry of water into the hydrocarbon region, and the effect of this barrier is naturally proportional to the fractional “free” area (

*A*−

*A*)/

_{0}*A*. This is the single most important feature in our model that will account for the major area dependence found by Mathai et al. (2007).

The second part of our model assumes that the hydrocarbon core, by itself, has a permeability *P _{C}* =

*KD*/

_{C}*d*, given by the simple SD model for the hydrocarbon core. In this simplest model that we will first consider, the only parameter that will vary between different lipid bilayers is the structural parameter

_{C}*d*, the thickness of the hydrocarbon core. Of course, it might be considered that the effective hydrocarbon core thickness for permeability could be smaller than

_{C}*d*due to tight packing of the first few methylene groups in the hydrocarbon chains (Subczynski et al., 1994; Xiang et al., 1998). One might also suppose that a larger fraction of “free volume,” (

_{C}*V*−

*V*)/

_{0}*V*, would increase the space available for water and thereby increase the partition coefficient

*K*. Larger fraction of free volume would also allow for more dynamical motion that would increase the intrinsic coefficient of diffusion

*D*. However, the volume per methylene is nearly constant for all the fluid phase lipids studied, so such a factor would make no difference between the different lipids we studied.

_{C}Therefore, this theory quantitatively predicts that, for pure lipid bilayers, the dependence of *P* on structural parameters is given by(3)At this point the linear factors α and γ are just fitting parameters that are assumed only to be independent of the structural quantities *A*, *A*_{0}, and *d _{C}* whose postulated dependencies are explicitly displayed in Eq. 3. Of course, α and γ must be affected by

*K*and by the coefficients of diffusion that may be different in different parts of the bilayer, as will be seen in the following two sections. In first approximation, α and γ will be assumed to be the same for all fully fluid phase lipid bilayers. Fitting these formulae to permeability data for five lipid bilayers with different structural parameters therefore determines α and γ from which the individual permeabilities

*P*and

_{H}*P*are determined for each of the bilayers.

_{C}The first question to investigate is whether both terms on the right hand side of Eq. 3 are significant. It has already been shown (Mathai et al., 2007) that setting α = 0, which is just the single layer SD model, is not adequate because there is a poor correlation of *P* with 1/*d _{C}*. The other extreme is to set γ = 0, which corresponds to the hydrocarbon core permeability

*P*being much greater than the interfacial permeability

_{C}*P*. The open squares in Fig. 2 show that this first term that involves the area

_{H}*A*already gives fairly good theoretical values; this reflects the point made by Mathai et al. (2007) that the best correlation of permeability is with

*A*. However, when γ is set to 0, the predicted permeability for the thickest bilayer diC22:1 is too large and the predicted permeability for the thinnest bilayer DLPC is too small. This discrepancy can clearly be alleviated by inclusion of the second term. The red circles in Fig. 2 show the best fit of the theory using Eq. 3. Inclusion of the second term does indeed alleviate the aforementioned discrepancy. The legend to Fig. 2 also shows that the values of the parameter

*A*that are given by the best fits are consistent with negligible permeability of the gel phase which has an area 48 Å

_{0}^{2}for PC bilayers. The somewhat larger values of

*A*

_{0}in the legend in Fig. 2 can be justified as accounting for the steric area of a water molecule. Another way that

*A*

_{0}could be increased for water transport is that “ethanol may block water diffusion pathways by occupying points of water entry into bilayers at the interface” (Huster et al., 1997).

Motivated by simulations (Marrink and Berendsen, 1994, 1996) and also by an electron spin resonance (ESR) result (Subczynski et al., 1994) that the hydrophobicity barrier is narrower than *d _{C}*, we have also investigated a variation of Eq. 3 that models an effective hydrocarbon thickness for permeability by replacing the factor

*d*in the second term by a factor (

_{C}*d*− δ). The green triangles in Fig. 2 show that the fit is slightly improved when δ = 15 Å and the fit continues to improve as δ increases to 76 Å. However, the physical absurdity of this last result, namely, that the effective hydrocarbon thickness (

_{C}*d*− δ) becomes strongly negative, suggests that adding the fourth fitting parameter δ is not warranted by the data. Indeed, artificially reducing

_{C}*P*just for DOPC by 10%, which is close to estimated experimental uncertainties, yields a value of δ close to zero.

Fig. 3 compares the partial permeabilities *P _{H}*/

*2*(which includes both interfaces) and

*P*for the hydrocarbon core for the last two combinations of the parameters shown in Fig. 2. Of course, the ratio

_{C}*P*/

_{C}*P*varies with different lipids due to their different structural properties. The ratio

_{H}*P*/

_{C}*P*also depends upon the choice of effective thickness

_{H}*d*− δ. For both values of δ shown in Fig. 3,

_{C}*P*/

_{H}*2*is smaller than

*P*, so passage through the headgroup regions is predicted to be the slower process. Nevertheless,

_{C}*2P*/

_{C}*P*is less than 10 for the thinnest DLPC bilayer and is less than 3 for the thickest diC22:1 bilayer, so the hydrocarbon core permeability plays a role, even though it is secondary to the role played by the headgroup regions.

_{H}### II. Two Detailed Models

The preceding section did not address the very important question regarding the role played by the partition coefficient *K* that is crucial in order for a theory to obey Overton's rule. This section analyzes two mathematical models that answer this question. The two models also predict values for the two linear parameters α and γ in Eq. 3 and this could, in principle, reduce the number of free parameters for fitting data. However, it is important to consider both models because the predicted formula for γ is different. The difference shows that this result of mathematical modeling is not robust, so this comparison prevents the drawing of unwarranted numerical conclusions.

For both mathematical models we will refer to Fig. 4 for the free energy landscape which is the local (noncratic) part of the chemical potential. The free energy of water is assumed to be high and constant in the hydrocarbon core and low in the water. These two regions are separated by the interfacial headgroup regions, which are generally quite complicated. For simplicity, linear forms for the free energy will be assumed. It may be noted that this free energy landscape is qualitatively similar to the hydrophobicity landscapes obtained from spin labeling experiments (Subczynski et al., 1994) and from simulations (Fig.6 of Marrink and Berendsen,1994).

#### A: Chemical Kinetics Model.

As advocated long ago by Zwolinski et al. (1949), one may consider a chemical kinetics description of transport and permeability. The simplest mathematical way to describe the physical model shown in Fig. 1 employs four states as shown in Fig. 4. The bulk water phases are represented by states 1 and 4 with concentrations c_{1} and c_{4}, respectively. The hydrocarbon core is represented by states 2 and 3 with concentrations c_{2} and c_{3}, respectively. The physical locations of states 2 and 3 are just within the ends of the hydrocarbon core closest to the bulk water states 1 and 4, respectively. The distance between states 2 and 3 is the thickness d_{C} of the hydrocarbon region. The distance between states 1 and 2 (and between states 3 and 4) is the thickness d_{H} of the interfacial headgroup region. The kinetics of water or other solute flow through the membrane are given by the first order kinetics scheme(4)where the forward rate constants between the states can be written *k*_{12}, *k*_{23}, and *k*_{34}, and the backward rate constants are *k*_{21}, *k*_{32}, and *k*_{43} as shown in Fig. 4. The ratios of backward and forward rate constants are given by equilibrium free energy considerations. For symmetric lipid bilayers(5)where *K* = exp(−βΔ*F*) is the partition coefficient for water in the hydrocarbon core. It will be convenient to use the simplified notation,(6)

In steady state, all concentrations *c _{i}* are constant in time. The net forward currents between pairs of contiguous states are given by(7a)(7b)(7c)

In steady state, *J*_{12} = *J*_{23} = *J*_{34} = *J*. Addition of *J*_{12} and *J*_{34} followed by elimination of *c _{2}* −

*c*using Eq. 7b then gives(8)with(9)

_{3}Correspondence with Eq. 2 in the text follows by identifying the hydrocarbon core permeability *P _{C}* =

*d*and the headgroup permeability

_{C}k_{C}K*P*=

_{H}*d*. Of course,

_{H}k_{H}K*P*is usually written as

_{C}*KD*/

_{C}*d*and this identifies the coefficient of diffusion in the hydrocarbon regions as

_{C}*D*=

_{C}*d*

_{C}^{2}

*k*, which is the usual formula from random walk theory that gives the coefficient of diffusion as the hopping distance squared divided by hopping time. We next recognize that

_{C}*k*should contain the obstruction factor (

_{H}*A*−

*A*

_{0})/

*A*, which we wish to display explicitly. The local coefficient of diffusion

*D*within the unobstructed part of the headgroup region, that should be comparable numerically to

_{H}*D*, should not contain an area- dependent factor. It is then given as

_{C}*D*=

_{H}*d*

_{H}^{2}

*k*/(

_{H}A*A*−

*A*) because

_{0}*k*contains the factor (

_{H}*A*−

*A*

_{0})/

*A.*We therefore have(10)

There are two differences between the preceding kinetic modeling and that of Zwolinski et al. (1949). The first is unimportant. They included *m* − 1 intermediate states in the hydrocarbon region between our states 2 and 3, but as they showed in their Eq. 33, if the additional rate constants are all equal, corresponding to a homogeneous region, and each is scaled by the appropriate multiple of our *k _{C}*, there is no difference in the final equation for the permeability. The second difference is quite important. Zwolinski et al. (1949) supposed a large free energy barrier to entry of water into the hydrocarbon region in addition to the increase in free energy |ΔF| shown in Fig. 4. In their presentation they did not display the factor of

*K*that must be present even when there is the extra barrier they assumed. In our presentation we have not included any extra free energy barrier. This means that our

*k*=

_{H}*k*

_{21}=

*k*

_{34}models transition over negligible barriers into states with considerably lower free energies, so

*k*should not depend upon

_{H}*K*. In the Eyring absolute rate theory (Glasstone et al., 1941) when there is no barrier,

*k*=

_{H}*kT/h*, where

*h*is Planck's constant,

*kT*is thermal energy, and the entire

*K*dependence resides in the rate constants

*k*

_{12}=

*k*

_{43}=

*Kk*. Up to this point, our free energy profile across the bilayer has the shape of a mesa with a high flat plateau in the hydrocarbon region with steeply sloping sides into the low plains for the bulk water (Fig. 4). Our innovation in Section I is that, rather than imposing an extra free energy barrier, we impose a geometric obstruction factor on

_{H}*P*, given in Eqs. 3 and 10, that impedes diffusion through a fraction of the bilayer area. This factor may be thought of as a high picket fence on the mesa slope where the pickets represent the headgroup obstructions schematically shown in Fig. 1.

_{H}The most serious objection to the model as developed by Zwolinski et al. (1949) comes from reconciling it to Overton's rule. To effect such a reconciliation following their discussion of their Eq. 34, one would have to conclude that *P _{C}* had to be the rate-limiting step for permeability, as they did on their page 1444. In contrast, our presentation has a factor of K in both

*P*and

_{H}*P*, and therefore in

_{C}*P*, so it satisfies Overton's rule without forcing

*P*to be rate limiting. Furthermore, it allows a strong area dependence by making solute entry into the hydrocarbon core (up a mesa slope) slower than diffusion through the hydrocarbon core (across a flat mesa).

_{C}It may also be noted that Dix et al. (1978) discussed a three layer model in the mathematical framework of Zwolinski et al. (1949). However, they ended their paper with the opposite conclusion, namely, that the rate limiting step was the interfacial resistance and that 2/*P _{H}* was higher by several orders of magnitude than diffusional resistance 1/

*P*within the hydrocarbon core. While closer to our conclusion, our Fig. 3 has the ratio within a factor of 10 for fluid phase lipid bilayers. The conclusion of Dix et al. (1978) was based on residency times of water of 100 ns in the membrane. However, it is well known that PC lipid headgroups bind at least one or two water molecules so tightly that they are difficult to remove even by extensive drying (Jendrasiak and Hasty, 1974). We suggest that these strongly bound waters may account for the long residency. To include this in a kinetic model, a state 2b would be added to the left headgroup region that would have a maximum capacity of a few water molecules per lipid and would have very low free energies. State 2b would not be on the linear pathway shown in Fig. 4. Rather, it could be on an alternative branched pathway between states 1 and 2 or it could just be a dead end side path connected only to state 1 or to state 2. As such, it, and its symmetrically equivalent 5b state, would hardly perturb the previous analysis while providing an explanation for long residency times for water molecules in a nonbulk water environment.

_{C}This chemical kinetics model makes specific predictions about the two linear parameters in the general theory in Section II. Comparing Eq. 10 with Eq. 3 and the two components *P _{H}* and

*P*defined in Eq. 2 gives(11a)(11b)Assuming that

_{C}*D*= 2 × 10

_{C}^{−5}cm

^{2}/s, Eq. 11b gives

*K*= 5.4 × 10

^{−4}from the value of γ for the δ = 0 case in Fig. 2 and

*K*= 4.3 × 10

^{−4}for the δ = 15 Å case. For comparison, the partition coefficient for water in hexadecane is 4.2 × 10

^{−5}(Walter and Gutknecht, 1986). Then, if we also assume that

*D*= 2 × 10

_{H}^{−5}cm

^{2}/s, Eq. 11a gives the thickness of the headgroup region to be

*d*= 6.1 Å for either δ = 0 or δ = 15 Å. These are quite reasonable values of

_{H}*K*and

*d*that could be further tuned by modest changes in

_{H}*D*and

_{H}*D*. For example, if we arbitrarily set

_{C}*D*= 10

_{C}^{−5}cm

^{2}/s and

*D*= 1.5 × 10

_{H}^{−5}cm

^{2}/s, then

*K*∼ 0.001 and

*d*= 9.3 Å, which is close to the thickness of the interfacial headgroup region (Nagle and Tristram-Nagle, 2000).

_{H}#### B: Continuum Model.

As a mathematical model, the chemical kinetics model in the previous subsection is rather primitive because the interfacial headgroup region is represented only by one reaction pathway involving only two states, one at each edge of the region. One can ask what the effect would be to have additional states on a linear kinetics pathway within the headgroup region, and the answer is that the final equations change. Rather than adding a few more states, it is more efficient to proceed to the opposite extreme that consists of an infinite number of states; this is the continuum model.

The continuum model is treated by the Nernst-Planck extension of Fick's law for diffusion. Let *x* be the position perpendicular to the bilayer and let *x _{i}* be the particular values for the positions labeled

*i*= 1,2,3,4 in Fig. 4, so the headgroup thickness

*d*=

_{H}*x*−

_{2}*x*=

_{1}*x*−

_{4}*x*and the hydrocarbon core thickness

_{3}*d*=

_{C}*x*−

_{3}*x*. Let the free energy difference

_{2}*F*(

*x*) −

_{2}*F*(

*x*) be Δ

_{1}*F*and the magnitude of the corresponding force be

*f*= Δ

*F*/

*d*, noting that

_{H}*f*is negative when

*x*

_{1}<

*x*<

*x*

_{2}. Let

*c*(

*x*) be the concentration of solute and β be the inverse thermal energy 1/

*kT*. Then, for steady state the solute current

*J*is constant as a function of

*x*and is given by the Nernst-Planck equation(12)

where *D* is the coefficient of diffusion. It has been emphasized that *D* should be a nonconstant function of *x* (Diamond and Katz, 1974; Marrink and Berendsen, 1994), but to keep the model reasonably simple and calculable, we will assume a constant *D _{C}* in the hydrocarbon chain region

*x*

_{2}<

*x*<

*x*

_{3}where

*f*= 0. In the headgroup regions,

*x*

_{1}<

*x*<

*x*

_{2}and

*x*

_{3}<

*x*<

*x*

_{4}, it is convenient to factor

*D*into the headgroup obstruction factor

*(A*−

*A*

_{0})/

*A*and a coefficient of diffusion

*D*in the unobstructed part of the region, with a value of

_{H}*D*that is comparable to

_{H}*D*.

_{C}To obtain the permeability, *c*(*x*) is first noted to have the following forms in the three separate regions(13a)(13b)(13c)

where *b* = *(Jd _{H}*/βΔ

*FD*)(

_{H}*A*/(

*A*−

*A*

_{0})). The parameters

*a*and

_{1}*a*are related to the known concentration differences in the bulk phase by(14a)and also to(14b)Elimination of (

_{4}*a*−

_{1}*a*) then gives(15)where

_{4}*K*= exp(−βΔ

*F*) is the partition coefficient. The factor in square brackets is just

*1/P*by the definition of permeability and the inverses of the two individual terms therein can be identified asand(16)

The result in Eq. 16 is identical to Eq. 10 for the chemical kinetics model except for the final factor (−*ln*(*K*)/(1 − *K*)) in *P _{H}*. This factor depends only weakly on

*K*, varying by only about one order of magnitude as

*K*varies by five orders of magnitude for hydrophilic solutes with

*K*< 0.1, so the basic Overton rule dependence of

*P*on

*K*continues to hold.

We next follow the discussion in the last paragraph of the previous subsection. Again, assuming that *D _{C}* = 2 × 10

^{−5}cm

^{2}/s,

*P*in Eq. 16 gives

_{C}*K*= 5.4 × 10

^{−4}from the value of γ for the δ = 0 case in Fig. 2. But if we also assume that

*D*= 2 × 10

_{H}^{−5}cm

^{2}/s, then Eq. 16 gives the thickness of the headgroup region to be

*d*= 46 Å, which is clearly an unphysically large value. However, setting

_{H}*D*= 0.4 × 10

_{H}^{−5}cm

^{2}/s obtains a structurally acceptable value of

*d*= 9 Å. It may be noted that the simulation of Marrink and Berendsen (1994) gives a smaller coefficient of diffusion in the headgroup region than in the hydrocarbon core region.

_{H}## Discussion

The general phenomenological theory presented in Section I was motivated by the correlation of recently measured water permeability (Mathai et al., 2007) with the structure of lipid bilayers. At the core of this theory is a free area factor (*A* − *A*_{0})/*A*, introduced in Eqs. 3, 10, and 16, that is open for permeation. Free area and free volume concepts have been criticized when the free quantities are much smaller than molecular sizes (Edholm and Nagle, 2005). However, the free area concept gains traction when the quantized open area is larger than the area of a water molecule, as it is for typical water pores. This is also the case for the quantity *A* − *A*_{0}, which is the open space locally available in our theory and which is not much smaller than water molecules.

While quite general, it is important that this essentially postulated theory be consistent with more specific, microscopic models and calculations. Section II shows that there are at least two different microscopic models from which the phenomenological theory is derivable. The phenomenological theory in Eq. 3 also did not explicitly include any role for the partition coefficient *K* or coefficients of diffusion, but this is provided by the detailed models. Both the chemical kinetics model (Eq. 10) and the continuum Nernst-Planck model (Eq. 16) have a linear *K* factor in both the headgroup permeability *P _{H}* and in the chain permeability

*P*In contrast to the coefficients of diffusion, which can be different in the core and headgroup regions, there is only one partition coefficient given in Eq. 5 by the Boltzmann factor

_{C}.*K*= exp(−βΔ

*F*) for the free energy difference Δ

*F*of the solute in the hydrocarbon core versus water. This is an important result because it shows that a three layer theory is consistent with Overton's rule. The three layer theory also removes the discrepancy that the hydrocarbon core thickness is too large in the single layer solubility-diffusion theory. Fig. 3 shows that

*P*can be quite large as is required in order to have realistic values of

_{C}*d*because the experimental permeability is primarily determined in Eq. 2 by the smaller

_{C}*P*, which provides the greater resistance.

_{H}The theoretical result for the continuum model (Eq. 16) is different from the chemical kinetics model (Eq. 10) by having a weakly varying logarithmic K factor in the headgroup permeability *P _{H}*. The last paragraphs of the two subsections in Section II show that either model leads to reasonable results for the thickness of the headgroup region

*d*provided that the unknown coefficients of diffusion

_{H}*D*and

_{H}*D*are chosen appropriately. However, because of the ln

_{C}*K*factor in the continuum model, the ratio

*D*is different for the two models. The smaller value of

_{H}/D_{C}*D*required for the continuum model is consistent with the presence of local free energy minima within the heterogeneous headgroup region that could trap the solute for periods of time long compared with free diffusion in the more homogeneous hydrocarbon chain environment as suggested by Marrink and Berendsen (1994). While quite plausible, our results may not warrant such a firm conclusion. We assumed in the continuum model that the free energy profile is linear in the headgroup region (Eq. 13), but this leads to an exponential water concentration profile, whereas computer simulations suggest a more nearly linear water profile, as indicated in Fig. 1. Any continuum model requires detailed assumptions about the free energy profile that can be quite complicated and uncertain and obscure the main ideas, so we have chosen not to pursue variations of the continuum model. The chemical kinetics model avoids such complications by incorporating all the details of the headgroup region into a single rate constant, which has the merit of simplicity. Chemical kinetics models also allow for easy variations in the free energy landscape to treat detailed aspects of other solutes, as shown in the online supplemental materials (available at http://www.jgp.org/cgi/content/full/jgp.200709849/DC1).

_{H}/D_{C}All the fitting to water permeability data in this paper assumed that the partition coefficient *K* is the same in the five lipid bilayers. One might suppose that *K* for water would be larger for lipids with more polarizable unsaturated double bonds, as appears to be the case for polyunsaturated lipids (Huster et al., 1997; Olbrich et al., 2000). This would account for the theoretical permeability being too low for DOPC in Fig. 2 but it would make the fit worse for diC22:1PC. Also, electron spin resonance (ESR) measurements suggest that DOPC is more, rather than less, hydrophobic than lipids with saturated chains (Subczynski et al., 1994), so we have chosen not to allow variations in K, which is consistent with all the lipids having the same density of packing, i.e., the same *V _{CH2}*.

The theory as presented uses average structural quantities, such as the average area *A* of the headgroups. Of course, there are fluctuations in the local *A* in the fluid phase of bilayers, and the permeability will be transiently enhanced locally when *A* fluctuates to a larger value. Indeed, it has been suggested that the anomalously large permeability of bilayers to Na^{+} ions near the main chain melting phase transition is due to the nonlinear effect of greater fluctuations in the local area that must occur when the lateral area modulus *K _{A}* becomes small near a higher order phase transition (Nagle and Scott, 1978). However, none of the bilayers discussed here were in critical regions near the chain-melting transition temperature and all had values of

*K*that were substantially the same (Rawicz et al., 2000). The lack of empirical correlation of

_{A}*P*with

*K*(Mathai et al., 2007) suggests that average structural quantities suffice.

_{A}The bilayers used in Figs. 2 and 3 all had the same headgroup. Water permeability data for DLPE and DOPS are also presented by Mathai et al. (2007) and compared with structural data. Of course, different head groups should require different values of *A*_{0} and possibly different values of the coefficient of diffusion *D _{H}* in the headgroup region, so data from at least two different lipids with the same headgroup are required to obtain both parameters to enable a comparison to the PC lipids. Since we do not have those data, let us assume that

*D*is the same as for PC lipids. Then, the values of

_{H}*A*

_{0}required to match theory, using Eq. 3, to experiment are

*A*

_{0}= 51.2 Å

^{2}for DOPS and

*A*

_{0}= 50.1 Å

^{2}for DLPE. As would be expected, these values are smaller than the

*A*

_{0}= 53.6 Å

^{2}given in Fig. 2 for the PC headgroups, but they are not as much smaller as would be expected by the gel phase areas that are 41.0 Å

^{2}for DLPE (McIntosh and Simon, 1986) and 40.8 Å

^{2}for DMPS (Petrache et al., 2004), ∼7 Å

^{2}less than the 47–48 Å

^{2}for PC headgroups. However, compared with PC headgroups, PE and PS headgroups have additional hydrogen bonding opportunities that could be modeled either as blocking some of the area available for water permeation (

*i.e*., increasing

*A*

_{0}) or as providing local traps that would reduce

*D*(Marrink and Berendsen, 1994). Water permeability and structural data for DOPC with 10, 20, and 40% cholesterol were also presented by Mathai et al. (2007). Incorporation of cholesterol into our theory requires additional choices. Cholesterol might additionally obstruct entry of the water into the hydrocarbon region, or it might not, according to the theory of Huang and Feigenson (1999) that the headgroups shield the cholesterol from water. Also, the rigid ring structure of cholesterol might obstruct the diffusion within the hydrocarbon region. With enough cholesterol, the hydrocarbon chains become more ordered, like a gel phase, and less mobile, so

_{H}*D*might become smaller. Furthermore, it has been suggested that

_{C}*K*should be reduced by cholesterol (DeYoung and Dill, 1990; Xiang and Anderson, 1997), as seems plausible as the phase becomes liquid ordered instead of fully fluid. These are issues that are difficult to model, and we have chosen not to include cholesterol data in the fits in this paper. However, if we assume that α and γ in Eq. 3 are the same as for fully fluid phase lipids, then the values of

*A*

_{0}required to match theory and experiment in Fig. 2 are

*A*

_{0}= 53.1 Å

^{2}for DOPC with 10% cholesterol,

*A*

_{0}= 55.2 Å

^{2}for DOPC with 20% cholesterol, and

*A*

_{0}= 58.0 Å

^{2}for DOPC with 40% cholesterol.

While this theory has been motivated by water permeability measurements and while the tests presented use only these data, we suggest that the general theory may apply more generally to other solutes. Two classes of solute are considered in detail in the online supplemental material (http://www.jgp.org/cgi/content/full/jgp.200709849/DC1). The first is solutes, like acetic acid, that have been suggested to have strong binding to the interfacial region of bilayers (Xiang and Anderson, 1995). The second class is hydrophobic solutes whose partition coefficients into oil are greater than unity. We suggest that studies with different solutes concentrate primarily on bilayers with lipids that share the same headgroup and whose structures have been determined. Even with this constraint, one should expect some of the parameters and even the underlying free energy landscapes to be different from Fig. 4, as discussed in the online supplemental material.

Even homogeneous lipid bilayers have more complexity than can readily be included in a simple theory for passive permeability, so perfect agreement with experiment is not a realistic goal. As was emphasized by Diamond and Katz (1974) and mentioned many times since, the most realistic models would include partition coefficients and coefficients of diffusion that would vary continuously through the bilayer. However, an appropriate goal should still be a simple theory that can provide insight while accommodating the most significant permeability data with a reasonably small number of measurable parameters. With more precise structural data on lipid bilayers now available (Mathai et al., 2007), we believe that it is warranted to return to the approach of Zwolinski et al. (1949) and Diamond and Katz (1974) and try to improve the theory beyond the single layer solubility-diffusion model while stopping short of the continuous description with infinitely many parameters. We offer the present three layer theory, which should be tested further experimentally with other solutes and with other lipid bilayers when their structures are determined.

## Acknowledgments

This research was supported by the National Institutes of Health, grants GM 44976 (J.F. Nagle) and DK43955 (M.L. Zeidel, J.C. Mathai, and S. Tristram-Nagle).

Olaf S. Andersen served as editor.

## Footnotes

Abbreviations used in this paper: PC, phosphatidylcholine; SD, solubility-diffusion.

- Submitted: 26 June 2007
- Accepted: 7 December 2007