# Theory of Lipid Polymorphism: Application to Phosphatidylethanolamine and Phosphatidylserine

## I Introduction

Biological lipids in solution display several different lyotropic phases, and the implications this may have for biological function has been a subject of speculation for many years (Cullis et al., 1985; de Kruijff 1997). Lipid phase behavior depends upon several factors, some of which are intrinsic to the lipid architecture itself. For example, an increase in the length of the hydrocarbon tails brings about transitions from lamellar, $`L_\alpha `$, to inverted hexagonal, $`H_{II}`$, phases (Seddon, 1990), while an increase in the volume of the headgroup brings about the reverse (Gruner, 1989). Other factors regulating phase behavior are externally controlled, such as temperature, solvent concentration, and solvent pH (Hope and Cullis, 1980, Seddon et al., 1983, Bezrukov et al., 1999). It is these factors which are the focus of this paper.

Lipid phase behavior has been addressed extensively by the construction of phenomenological free energy functions which contain terms describing, inter alia, bending, hydration, and interstitial energies (Helfrich, 1973; Kirk et al., 1984; Rand and Parsegian, 1989; Kozlov et al., 1994). Such approaches, which obtain their several parameters from experimental measurement of various quantities, are quite useful, particularly in correlating phase behavior with other thermodynamic properties. Nonetheless, it would clearly be desirable to derive all thermodynamic quantities, including the phase behavior, by applying statistical mechanics to a microscopic model of the system. In addition to simplifying the description considerably, such approaches would correlate phase behavior with the the architectural properties of the lipid itself and its solvent.

Analytic, mean-field, approaches of statistical mechanics have been applied to anhydrous lipids to investigate behavior of increasing complexity. Such methods have been combined with realistic models of lipid tails to determine how the hydrocarbon chains pack in aggregates and in bilayers (Marcelja, 1974; Gruen, 1981 and 1985; Ben-Shaul et al. 1985; Fattal and Ben-Shaul, 1994). Results for the bilayer are in good agreement with molecular dynamic simulation (Heller et al., 1993). These methods have shown that in a neutral, anhydrous system, the entropy of the lipid tails always favors the $`H_{II}`$ over the $`L_\alpha `$ phase, and that a change in area per headgroup could bring about a transition between them (Steenhuizen et al., 1991).

Aggregates, such as the lipid bilayer, in the presence of solvent have also been considered within the mean-field approach applied to lattice models (Leermakers and Scheutjens, 1988). In addition to the tails, one must now model the solvent and the headgroups, and phosphatidylcholine and phosphatidylserine headgroups are among those which have been described (Meijer et al., 1994). The method is flexible and has been applied to many different systems, including bilayers with trans-membrane guest molecules (Leermakers et al., 1990). Results are quite good, with the exception that the local volume fraction of solvent inside the bilayer is rather large, several orders of magnitude greater than that observed in experiment (Jacobs and White, 1989). Lattice models, however, are not well-suited to the description of transitions between phases of different symmetry.

Recently similar methods were applied to a system of solvent and monoacyl lipid embedded in a continuous space, and the phase diagram was obtained by solving the mean-field theory exactly (Müller and Schick, 1998). It displayed both $`L_\alpha `$ and $`H_{II}`$ phases, so that the transition between them could be studied as a function of lipid architecture. The dependence of the transition on the architectural parameters, length of tail and volume of headgroup, was that observed in experiment. However, the fraction of solvent within the bilayers was again too large.

In this paper, we introduce a computationally more tractable model of a lipid than that employed by Müller and Schick, one whose hydrocarbon tails are modeled as flexible chains rather than within the rotational isomeric states framework (Flory, 1969; Mattice and Suter, 1994) employed earlier. We first study the model with an uncharged headgroup. Its phase behavior, both with respect to variations in architecture and variations in solvent concentration, is as expected, and in agreement with experiment. In particular, choosing model parameters appropriate to dioleoylphosphatidylethanolamine (DOPE), we obtain a phase diagram similar to that observed (Gawrisch et al. 1992, Kozlov et al. 1994). We extract the variation with temperature and solvent concentration of the lattice parameter of the inverted hexagonal phase, and compare it to experiment (Tate and Gruner, 1989, Rand and Fuller, 1994). The agreement is excellent. We also find that the concentration of solvent within the bilayer is vanishingly small. We then allow the headgroup to be negatively charged, introduce counter ions into the system, and include the Coulomb interaction between all charges. We show that the Coulomb repulsion between the headgroups does not tend to stabilize the $`L_\alpha `$ phase at the expense of the $`H_{II}`$, but has the opposite effect. Lastly we turn on a short-ranged interaction between charges and neutral solvent, an interaction which models the thermally averaged interaction between charges and the dipole of water. We then find that as the charge on the headgroup is turned on, the $`L_\alpha `$ phase is stabilized with respect to the $`H_{II}`$. In effect, as the charge on the headgroup increases, so too do the waters of hydration. In addition, the counter ions which are attracted to the headgroup are also enlarged by their own waters of hydration. It is the totality of these waters which effectively increases the headgroup volume and therefore stabilizes the lamellar phase.

The paper is organized as follows. In the next section, we introduce the model for the charged lipid, the solvent, and counter ions, specify all the interactions between them and set up the partition function of the system. In section III, we derive the self-consistent field theory for it. At the heart of the theory are four self-consistent equations for the electrostatic potential of the system and the three effective fields which determine the headgroup, tail, and solvent densities. One of these self-consistent conditions is simply the non-linear Poisson-Boltzmann equation. In section IV, we expand all functions of position into a complete set of functions having a specified space-group symmetry, and rewrite the self-consistent equations in terms of the coefficients of these expansions. These equations are solved numerically, and the free energies of the various phases computed. A comparison of the free energies yields the phase diagram. In Section V, we present the phase diagram for the neutral lipid as a function of temperature and one architectural parameter. We include here only the classical phases, lamellar, inverted and normal hexagonal, and inverted and normal body-centered-cubic, as well as the disordered phase. For the remainder of the section, we choose an architecture such that the anhydrous, neutral lipid orders into the $`H_{II}`$ phase. Results for the system in the presence of a neutral solvent, along with comparisons to experiment, are presented next.

In Section VI, we consider the charged lipid. We choose a water concentration such that the neutral lipid remains in the $`H_{II}`$ phase. By varying the counter ion concentration, we turn on the charge on the headgroup, and thus all Coulomb interactions, but we keep the short-ranged interaction between charges and neutral solvent set to zero. We show that the $`H_{II}`$ phase is stabilized with respect to the $`L_\alpha `$ rather than destabilized, a result completely analogous to that of polyelectrolytes (Nyrkova et al., 1994). Lastly we present our results for the case in which, in addition to the Coulomb interaction, a short-ranged interaction between charges and solvent is turned on. Here we find that the $`L_\alpha `$ is indeed stabilized with respect to the $`H_{II}`$ phase, in agreement with experiment (Hope and Cullis, 1980, Bezrukov et al., 1999).

## II The Model

We consider a system composed of charged lipids, neutral solvent, and counter ions in a volume $`V`$. There are $`n_L`$ lipids, each of which consists of a head, with volume $`v_h`$, and two equal length, completely flexible, tails each consisting of $`N`$ segments of volume $`v_t`$. Each lipid tail is characterized by a radius of gyration $`R_g=(Na^2/6)^{1/2}`$, with $`a`$ the statistical segment length. The heads carry a negative charge $`eQ_h`$. The solvent consists of $`n_s`$ neutral particles of volume $`v_s`$, while the $`n_c`$ counter ions have charge $`+e`$ and negligible volume, $`v_c=0`$. There are five dimensionless densities which totally specify the state of the system; the number density of the headgroups, $`\widehat{\mathrm{\Phi }}_h`$, of the tail segments, $`\widehat{\mathrm{\Phi }}_t`$, and of the solvent, $`\widehat{\mathrm{\Phi }}_s`$, and the charge density of the headgroups, $`e\widehat{P}_h`$, and of the counter ions, $`e\widehat{P}_c`$. They can be written as

$`\widehat{\mathrm{\Phi }}_h(𝐫)`$ $`=`$ $`v_h{\displaystyle \underset{l=1}{\overset{n_l}{}}}\delta (𝐫𝐫_l(1/2)),`$ (1)   
$`\widehat{\mathrm{\Phi }}_t(𝐫)`$ $`=`$ $`v_h{\displaystyle \underset{l=1}{\overset{n_l}{}}}{\displaystyle _0^1}\delta (𝐫𝐫_l(s))𝑑s,`$ (2)   
$`\widehat{\mathrm{\Phi }}_s(𝐫)`$ $`=`$ $`v_h{\displaystyle \underset{j=1}{\overset{n_s}{}}}\delta (𝐫𝐑_{s,j}),`$ (3)   
$`\widehat{P}_h(𝐫)`$ $`=`$ $`v_h{\displaystyle \underset{l=1}{\overset{n_l}{}}}Q_{h,l}\delta (𝐫𝐫_l(1/2)),`$ (4)   
$`\widehat{P}_c(𝐫)`$ $`=`$ $`v_h{\displaystyle \underset{i=1}{\overset{n_c}{}}}\delta (𝐫𝐑_{c,i}).`$ (5)   
We have chosen $`v_h`$ as a convenient volume to make all densities dimensionless. In the above, $`𝐑_{s,j}`$ is the position of the $`j`$’th solvent particle, and $`𝐑_{c,i}`$ the position of the $`i`$’th counter ion. The configuration of the $`l`$’th lipid is described by a space curve $`𝐫_l(s)`$, where $`s`$ ranges from 0, at the end of one tail, through $`s=1/2`$ at which the head is located, to $`s=1`$, the end of the other tail. The nominal probability that the charge on the headgroup of the $`l`$’th lipid, $`eQ_{h,l}`$, is equal to $`e`$ or 0 is $`p`$ or $`1p`$ respectively. As we model the case in which charges can associate or disassociate from the headgroup, it will be necessary to average the partition function of the system with respect to the charge distribution. This corresponds to an annealed distribution in the nomenclature of Borukhov (Borukhov et al., 1998). The concentrations of lipid, solvent, and free counter ions are controlled by chemical potentials. In particular, increasing the number of free, positive, counter ions implies, by charge neutrality, an increase in the negative charge on the headgroups, and thus corresponds to an increase in the pH of the system.

The interactions among these elements are as follows. First there is a repulsive, contact interaction between headgroup and tail segments, and also between solvent and tail segments. The strength of the interaction is $`kTv_h\chi `$, where $`k`$ is Boltzmann’s constant and $`T`$ the absolute temperature. Second there is the Coulomb interaction between all charges. The dielectric constant of the solvent is denoted $`ϵ`$. Finally there is a contact interaction between all charges and the neutral solvent whose strength is $`kTv_h\lambda `$. This is to model the short-ranged, thermally averaged, interaction between charges and the dipole of water, an attractive interaction which decreases like $`r^4`$ and is of strength $`e^2u^2/6ϵ^2kT`$, where $`u`$ is the dipole moment of water (Israelachvili, 1985). Thus the energy per unit volume of the system, $`E/V`$, can be written

$`{\displaystyle \frac{v_h}{kT}}{\displaystyle \frac{E}{V}}[\widehat{\mathrm{\Phi }}_h,\widehat{\mathrm{\Phi }}_t,\widehat{\mathrm{\Phi }}_s,\widehat{P}_h,\widehat{P}_c]`$ $`=`$ $`2\chi N{\displaystyle \frac{d𝐫}{V}[\widehat{\mathrm{\Phi }}_h(𝐫)+\widehat{\mathrm{\Phi }}_s(𝐫)]\widehat{\mathrm{\Phi }}_t(𝐫)}`$ (8)   
$`+{\displaystyle \frac{\beta ^{}}{8\pi }}{\displaystyle \frac{d𝐫}{V}\frac{d𝐫^{}}{R_g^2}[\widehat{P}_h(𝐫)+\widehat{P}_c(𝐫)]\frac{1}{|𝐫𝐫^{}|}[\widehat{P}_h(𝐫^{})+\widehat{P}_c(𝐫^{})]}`$   
$`\lambda {\displaystyle \frac{d𝐫}{V}\widehat{\mathrm{\Phi }}_s(𝐫)[\widehat{P}_c(𝐫)\widehat{P}_h(𝐫)]},`$   
where

$$\beta ^{}\frac{4\pi e^2R_g^2}{v_hϵkT}$$
(9)   
is a dimensionless measure of the strength of the Coulomb interaction. The grand partition function (Matsen, 1995) of the system is

$`𝒵`$ $`=`$ $`{\displaystyle \underset{n_l,n_c,n_s}{}}{\displaystyle \frac{z_l^{n_l}z_c^{n_c}z_s^{n_s}}{n_l!n_c!n_s!}}{\displaystyle \underset{l=1}{\overset{n_l}{}}\stackrel{~}{𝒟}𝐫_l\stackrel{~}{𝒟}Q_{h,l}\underset{i=1}{\overset{n_c}{}}d𝐑_{c,i}\underset{j=1}{\overset{n_s}{}}d𝐑_{s,j}}`$ (10)   
$`\times `$ $`\mathrm{exp}\left\{E[\widehat{\mathrm{\Phi }}_h,\widehat{\mathrm{\Phi }}_t,\widehat{\mathrm{\Phi }}_s,\widehat{P}_h,\widehat{P}_c]/kT\right\}\delta (1\widehat{\mathrm{\Phi }}_h\gamma _s\widehat{\mathrm{\Phi }}_s\gamma _t\widehat{\mathrm{\Phi }}_t).`$ (11)   
Here $`\stackrel{~}{𝒟}𝐫_l`$ denotes a functional integral over the possible configurations of the $`l`$’th lipid and in which, in addition to the Boltzmann weight, the path is weighted by the factor $`𝒫[𝐫_{t,l}(s);0,1],`$ with

$$𝒫[𝐫,s_1,s_2]=𝒩\mathrm{exp}\left[\frac{1}{8R_g^2}_{s_1}^{s_2}𝑑s|\frac{d𝐫(s)}{ds}|^2\right],$$
(12)   
with $`𝒩`$ an unimportant normalization constant. The notation $`\stackrel{~}{𝒟}Q_{h,l}`$ denotes an integral over the probability distribution of the charge on the headgroup of the $`l`$’th lipid. We have enforced an incompressibility constraint on the system with the aid of the delta function $`\delta (1\widehat{\mathrm{\Phi }}_h\gamma _s\widehat{\mathrm{\Phi }}_s\gamma _t\widehat{\mathrm{\Phi }}_t)`$, where $`\gamma _s=v_s/v_h`$, and $`\gamma _t=2Nv_t/v_h`$. The latter parameter is the lipid architectural parameter. The relative volume of the headgroup with respect to that of the entire molecule is $`1/(1+\gamma _t)`$.

The model is now completely defined. The solvent is specified by $`\gamma _s`$, its volume per particle relative to that of the headgroup, and the architecture of the lipid is characterized by $`\gamma _t`$. There are three interactions, hydrophobic-hydrophilic, charge-charge, and charge-solvent, whose strengths are given by $`\chi `$, $`\beta ^{}`$, and $`\lambda `$ respectively. The external parameters are the temperature, conveniently specified in terms of a dimensionless temperature $`T^{}(2\chi N)^1`$, the fugacity of the solvent, $`z_s`$, and the fugacity of the free counter ions, $`z_c`$, which, by charge neutrality, controls the charge on the lipid headgroups. The characteristic length in the system is the radius of gyration, $`R_g`$. In the next two sections we derive the self consistent field theory for the model, first in real space, and then in the following section, in Fourier space.

## III Theory: real space

Evaluation of the partition function of Eq. 10 is difficult because the interactions are products of densities, each of which depends on the specific coordinates of one of the elements of the system. This dependence is eliminated in a standard way. We illustrate it on $`\widehat{\mathrm{\Phi }}_h(𝐫)`$ which, from its definition in Eq. 1, depends on the coordinates of the headgroup, $`𝐫_l(1/2)`$. One introduces into the partition function the identity

$`1`$ $`=`$ $`{\displaystyle 𝒟\mathrm{\Phi }_h\delta (\mathrm{\Phi }_h\widehat{\mathrm{\Phi }}_h)},`$ (13)   
$`=`$ $`{\displaystyle 𝒟\mathrm{\Phi }_h𝒟W_h\mathrm{exp}\left\{\frac{1}{v_h}W_h(𝐫)[\mathrm{\Phi }_h(𝐫)\widehat{\mathrm{\Phi }}_h(𝐫)]𝑑𝐫\right\}},`$ (14)   
in which $`\mathrm{\Phi }_h(𝐫)`$ does not depend on any specific coordinates of one of the elements of the system, but is simply a function of $`𝐫`$. The integration on $`W_h`$ extends up the imaginary axis. Inserting such identities for the five densities $`\widehat{\mathrm{\Phi }}_h`$, $`\widehat{\mathrm{\Phi }}_t`$, $`\widehat{\mathrm{\Phi }}_s`$, $`\widehat{P}_h`$, and $`\widehat{P}_c`$, and a similar identity for the delta function expressing the incompressibility condition, one rewrites the partition function, Eq. 10, as

$`𝒵`$ $`=`$ $`{\displaystyle 𝒟\mathrm{\Phi }_h𝒟W_h𝒟\mathrm{\Phi }_t𝒟W_t𝒟\mathrm{\Phi }_s𝒟W_s𝒟P_h𝒟U_h𝒟P_c𝒟U_c𝒟\mathrm{\Xi }}`$ (15)   
$`\times `$ $`\mathrm{exp}\{z_l𝒬_l[W_h,W_t,U_h]+z_c𝒬_c[U_c]+z_s𝒬_s[W_s]E[\mathrm{\Phi }_h,\mathrm{\Phi }_t,\mathrm{\Phi }_s,P_h,P_c]/kT\}`$ (16)   
$`\times `$ $`\mathrm{exp}\left\{{\displaystyle \frac{1}{v_h}}{\displaystyle \left[W_h\mathrm{\Phi }_h+W_t\mathrm{\Phi }_t+W_s\mathrm{\Phi }_s+U_hP_h+U_cP_c+\mathrm{\Xi }(1\mathrm{\Phi }_h\gamma _s\mathrm{\Phi }_s\gamma _t\mathrm{\Phi }_t)\right]𝑑𝐫}\right\},`$ (17)   
where

$$𝒬_l[W_h,W_t,U_h]=\stackrel{~}{𝒟}𝐫_l\stackrel{~}{𝒟}Q_h\mathrm{exp}\left\{W_h(𝐫_l(1/2))+Q_hU_h(𝐫_l(1/2))_0^1𝑑sW_t(𝐫_l(s))\right\},$$
(18)   
is the partition function of a single lipid in external fields $`W_h`$, $`W_t`$, and $`U_h`$,

$$𝒬_c[U_c]=𝑑𝐑_c\mathrm{exp}[U_c(𝐑_c)],$$
(19)   
is the partition function of a single counter ion of unit positive charge in an external potential $`U_c`$, and

$$𝒬_s[W_s]=𝑑𝐑_s\mathrm{exp}[W_s(𝐑_s)],$$
(20)   
is the partition function of a single solvent particle in the external field $`W_s`$. It is convenient to shift the zero of all chemical potentials so that $`z_l1/v_h`$, $`z_cz_c/v_h`$, and $`z_sz_s/v_h`$. The partition function, Eq. 15, can then be written in the form

$$𝒵=𝒟\mathrm{\Phi }_h𝒟W_h𝒟\mathrm{\Phi }_t𝒟W_t𝒟\mathrm{\Phi }_s𝒟W_s𝒟P_h𝒟U_h𝒟P_c𝒟U_c𝒟\mathrm{\Xi }\mathrm{exp}[\mathrm{\Omega }/kT],$$
(21)   
with

$`{\displaystyle \frac{v_h}{kTV}}\mathrm{\Omega }`$ $`={\displaystyle \frac{𝒬_l[W_h,W_t,U_h]}{V}}z_c{\displaystyle \frac{𝒬_c[U_c]}{V}}z_s{\displaystyle \frac{𝒬_s[W_s]}{V}}+{\displaystyle \frac{v_h}{kTV}}E[\mathrm{\Phi }_h,\mathrm{\Phi }_t,\mathrm{\Phi }_s,P_h,P_c]`$ (22)   
$``$ $`{\displaystyle \frac{d𝐫}{V}[W_h\mathrm{\Phi }_h+W_t\mathrm{\Phi }_t+W_s\mathrm{\Phi }_s+U_hP_h+U_cP_c+\mathrm{\Xi }(1\mathrm{\Phi }_h\gamma _s\mathrm{\Phi }_s\gamma _t\mathrm{\Phi }_t)]}.`$ (23)   
No approximations have been made to this point. What has been accomplished is a rewriting of the partition function from a form, Eq 10, in which all entities interact directly with one another, to a form, Eqs. 21 and 22, in which they interact indirectly with one another via fluctuating fields. Although the integrals in Eq. 21 over $`\mathrm{\Phi }_h,\mathrm{\Phi }_t,\mathrm{\Phi }_s,P_h,P_c`$ and $`\mathrm{\Xi }`$ could all be carried out, as they are no worse than Gaussian, the integrals over the fields $`W_h,W_t,W_s,U_h,`$ and $`U_c`$ cannot. Therefore we employ the self-consistent field theory in which we replace the integral in Eq. 21 by its integrand evaluated at its extremum. The values of $`W_h`$, $`\mathrm{\Phi }_h`$, etc. which satisfy the extremum conditions will be denoted by the corresponding lower case letters $`w_h`$, and $`\varphi _h`$, etc. The equations which determine them are six self-consistent equations for the six fields $`w_h`$, $`w_t`$, $`w_s`$, $`u_h`$, $`u_c`$, and $`\xi `$. They are

$`w_h(𝐫)`$ $`=`$ $`2\chi N\varphi _t(𝐫)+\xi (𝐫)`$ (24)   
$`w_t(𝐫)`$ $`=`$ $`2\chi N(\varphi _h(𝐫)+\varphi _s(𝐫))+\gamma _t\xi (𝐫)`$ (25)   
$`w_s(𝐫)`$ $`=`$ $`2\chi N\varphi _t(𝐫)\lambda (\rho _c(𝐫)\rho _h(𝐫))+\gamma _s\xi (𝐫)`$ (26)   
$`u(𝐫)`$ $``$ $`{\displaystyle \frac{u_h(𝐫)+u_c(𝐫)}{2}}`$ (27)   
$`=`$ $`{\displaystyle \frac{\beta ^{}}{4\pi }}{\displaystyle \frac{d𝐫^{}}{R_g^2}\frac{\rho _h(𝐫^{})+\rho _c(𝐫^{})}{|𝐫𝐫^{}|}}`$ (28)   
$`u_s(𝐫)`$ $``$ $`{\displaystyle \frac{u_h(𝐫)u_c(𝐫)}{2}}`$ (29)   
$`=`$ $`\lambda \varphi _s(𝐫)`$ (30)   
$`1`$ $`=`$ $`\varphi _h(𝐫)+\gamma _t\varphi _t(𝐫)+\gamma _s\varphi _s(𝐫).`$ (31)   
As the field $`\xi `$ is easily eliminated, the six equations readily reduce to five. The simplicity of Eq. 29 reduces this, in practice, to a set of four equations. The five densities $`\varphi _h`$, $`\varphi _t`$, $`\varphi _s`$, $`\rho _h`$, and $`\rho _c`$ are functionals of all of the above fields except $`\xi `$, and therefore close the cycle of self-consistent equations:

$`\varphi _h(𝐫)[w_h,w_t,u_h]`$ $`=`$ $`{\displaystyle \frac{\delta 𝒬_l[w_h,w_t,u_h]}{\delta w_h(𝐫)}}`$ (32)   
$`\varphi _t(𝐫)[w_h,w_t,u_h]`$ $`=`$ $`{\displaystyle \frac{\delta 𝒬_l[w_h,w_t,u_h]}{\delta w_t(𝐫)}}`$ (33)   
$`\varphi _s(𝐫)[w_s]`$ $`=`$ $`z_s{\displaystyle \frac{\delta 𝒬_s[w_s]}{\delta w_s(𝐫)}}`$ (34)   
$`=`$ $`z_s\mathrm{exp}[w_s(𝐫)]`$ (35)   
$`\rho _h(𝐫)[w_h,w_t,u_h]`$ $``$ $`{\displaystyle \frac{\delta 𝒬_l[w_h,w_t,u_h]}{\delta u_h(𝐫)}}`$ (36)   
$`\rho _c(𝐫)[u_c]`$ $`=`$ $`{\displaystyle \frac{\delta 𝒬_c[u_c]}{\delta u_c(𝐫)}}`$ (37)   
$`=`$ $`z_c\mathrm{exp}[u_c(𝐫)]`$ (38)   
The density $`\varphi _h(𝐫)`$ is simply the expectation value of $`\widehat{\mathrm{\Phi }}_h(𝐫)`$ in the single lipid ensemble. Similar interpretations follow for the other densities. Note that one of the self-consistent equations, Eq. 27, is simply the non-linear Poisson-Boltzmann equation, and $`u(𝐫)`$ the electric potential.

With the aid of the above equations, the mean-field free energy, $`\mathrm{\Omega }_{mf}`$, which is the free energy function of Eq. 22 evaluated at the mean-field values of the densities and fields, can be put in the form

$`\mathrm{\Omega }_{mf}`$ $`=`$ $`{\displaystyle \frac{kT}{v_h}}\left(𝒬_l[w_h,w_t,u_h]+z_c𝒬_c[u_c]+z_s𝒬_s[w_s]\right)+E[\varphi _h,\varphi _t,\varphi _s,\rho _h,\rho _c],`$ (39)   
$`=`$ $`kT(n_l+n_c+n_s)+E[\varphi _h,\varphi _t,\varphi _s,\rho _h,\rho _c],`$ (40)   
with $`E`$ given by Eq. 8. The thermodynamic potential, $`\mathrm{\Omega }`$, is that appropriate to an incompressible system calculated in the grand ensemble; the negative of the osmotic pressure multiplied by the volume. Thus the above equation states that the osmotic pressure is the sum of the ideal, partial osmotic pressures plus a correction due to the interactions. Within mean field theory, this correction is simply the energy per unit volume of the system.

We now specify that the charges in the system can associate with or disassociate from the headgroup in response to the local electrostatic potential. This implies that the partition function of a single lipid, $`𝒬_l`$ is to be averaged over the nominal charge distribution that $`Q_h=1`$ with probability $`p`$, and $`Q_h=0`$ with probability $`1p`$ (Borukhov et al., 1998). The consequence of this averaging is that $`𝒬_l[w_h,w_t,u_h]`$ of Eq. 18 becomes

$$𝒬_l[w_{h,eff},w_t]=\stackrel{~}{𝒟}𝐫_l\mathrm{exp}\left\{w_{h,eff}(𝐫_l(1/2))_0^1𝑑sw_t(𝐫_l(s))\right\},$$
(41)   
where

$`w_{h,eff}(𝐫)`$ $``$ $`w_h(𝐫)\mathrm{ln}{\displaystyle \stackrel{~}{𝒟}Q_h\mathrm{exp}[Q_hu_h(𝐫)]}`$ (42)   
$`=`$ $`w_h(𝐫)\mathrm{ln}[1+p(\mathrm{exp}[u_h(𝐫)]1)].`$ (43)   
Although this appears to introduce an unknown parameter $`p`$ into the problem, the condition of charge neutrality,

$$𝑑𝐫[\rho _h(𝐫)+\rho _c(𝐫)]=0,$$
(44)   
relates this parameter to the fugacity of the counter ions, $`z_c`$. In practice, we use this fugacity to control the pH and the amount of charge on the lipids.

There remains only to specify how the single-lipid partition function is obtained. One defines the end-segment distribution function

$$q(𝐫,s)=𝒟𝐫_l(s)\delta (𝐫𝐫_l(s))\mathrm{exp}\left\{_0^s𝑑t\left(\left[\frac{1}{8R_g^2}|\frac{d𝐑(t)}{dt}|^2\right]+w_{h,eff}(𝐫_l(t))\delta (t1/2)+w_t(𝐫_l(t))\right)\right\},$$
(45)   
which satisfies the equation

$$\frac{q(𝐫,s)}{s}=2R_g^2^2q(𝐫,s)[w_{h,eff}(𝐫)\delta (s1/2)+w_t(𝐫)]q(𝐫,s),$$
(46)   
with initial condition

$$q(𝐫,0)=1.$$
(47)   
The partition function of the lipid is then

$$𝒬_l=𝑑𝐫q(𝐫,1).$$
(48)   
From this expression for the single-lipid partition function and Eqs. 32, 33, and 36, one obtains expressions for the local density of the lipid heads

$$\varphi _h(𝐫)=\mathrm{exp}[w_{h,eff}(𝐫)]q(𝐫,\frac{1}{2})q(𝐫,\frac{1}{2}),$$
(49)   
of the lipid tails

$$\varphi _t(𝐫)=_0^1𝑑sq(𝐫,s)q(𝐫,1s)$$
(50)   
and of the charge density on the lipid heads

$$\rho _h(𝐫)=\frac{p\mathrm{exp}[u_h(𝐫)]}{1+p(\mathrm{exp}[u_h(𝐫)]1)}\varphi _h(𝐫).$$
(51)   
To summarize: there are four self-consistent equations to be solved for the fields $`w_h`$, $`w_t`$, $`w_s`$, and electrostatic potential $`u`$. These equations, obtained from simple algebraic manipulation of Eqs. 24 to 31, can be taken to be

$`\gamma _tw_h(𝐫)w_t(𝐫)`$ $`=`$ $`2\chi N[\gamma _t\varphi _t(𝐫)\varphi _h(𝐫)\varphi _s(𝐫)],`$ (52)   
$`\gamma _sw_h(𝐫)w_s(𝐫)`$ $`=`$ $`2\chi N(\gamma _s1)\varphi _t(𝐫)+\lambda (\rho _c(𝐫)\rho _h(𝐫)),`$ (53)   
$`1`$ $`=`$ $`\varphi _h(𝐫)+\gamma _t\varphi _t(𝐫)+\gamma _s\varphi _s(𝐫),`$ (54)   
$`R_g^2^2u(𝐫)`$ $`=`$ $`\beta ^{}(\rho _h(𝐫)+\rho _c(𝐫)).`$ (55)   
Note that we have chosen here to write the Poisson-Boltzmann equation, Eq. 27, in its local, rather than its integral form. When the four fields are known, the corresponding densities follow from Eqs. 35, 38, 49, 50, and 51.

Rather than attempt to solve these equations in real space, a difficult task for the periodic phases in which we are interested such as $`H_{II}`$, we recast the equations in a form which makes straightforward their solution for a phase of arbitrary space-group symmetry (Matsen and Schick, 1994).

## IV Theory: Fourier Space

We note that the fields, densities, and the end point distribution function depend only on one coordinate $`𝐫`$. Therefore in an ordered phase, these functions reflect the space-group symmetry of that phase. To make this symmetry manifest in the solution, we expand all functions of position in a complete, orthonormal, set of functions, $`f_i(𝐫),i=1,2,3\mathrm{}`$, each of which have the desired space group symmetry; e.g.

$`\varphi _h(𝐫)`$ $`=`$ $`{\displaystyle \underset{i}{}}\varphi _{h,i}f_i(𝐫),`$ (56)   
$`\delta _{i,j}`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle 𝑑𝐫f_i(𝐫)f_j(𝐫)}.`$ (57)   
Furthermore we choose the $`f_i(𝐫)`$ to be eigenfunctions of the Laplacian

$$^2f_i(𝐫)=\frac{\lambda _i}{D^2}f_i(𝐫),$$
(58)   
where $`D`$ is a length scale for the phase. The functions for the lamellar phase are clear. They can be taken to be

$`f_1(𝐫)`$ $`=`$ $`1,`$ (59)   
$`f_i(𝐫)`$ $`=`$ $`\sqrt{2}\mathrm{cos}[2\pi (i1)x/D],i2,`$ (60)   
Expressions for the unnormalized basis functions for other space-group symmetries can be found in X-ray tables (Henry and Lonsdale, 1969) as they are intimately related to the Bragg peaks. In the tables cited, those for the hexagonal phase, space group (p6m) can be found on page 372, and that of the bcc phase, space group (Im3m) on page 524.

The four self-consistent equations become

$`\gamma _tw_{h;i}w_{t;i}`$ $`=`$ $`2\chi N[\gamma _t\varphi _{t;i}\varphi _{h;i}\varphi _{s;i}],`$ (61)   
$`\gamma _sw_{h;i}w_{s;i}`$ $`=`$ $`2\chi N(\gamma _s1)\varphi _{t;i}+\lambda (\rho _{c;i}\rho _{h;i}),`$ (62)   
$`\delta _{1,i}`$ $`=`$ $`\varphi _{h;i}+\gamma _t\varphi _{t;i}+\gamma _s\varphi _{s;i},`$ (63)   
$`{\displaystyle \frac{\lambda _iR_g^2}{D^2}}u_i`$ $`=`$ $`\beta ^{}(\rho _{h;i}+\rho _{c;i})`$ (64)   
To obtain the partition functions and densities, we proceed as follows. For any function $`G(𝐫)`$, we can define a symmetric matrix

$$(G)_{ij}\frac{1}{V}f_i(𝐫)G(𝐫)f_j(𝐫)𝑑𝐫$$
(65)   
Note that $`(G)_{1i}=(G)_{i1}=G_i`$, the coefficient of $`f_i(𝐫)`$ in the expansion of $`G(𝐫)`$. Matrices corresponding to functions of $`G(𝐫)`$, such as

$$\left(e^G\right)_{ij}\frac{1}{V}f_i(𝐫)e^{G(𝐫)}f_j(𝐫)d(𝐫),$$
(66)   
are evaluated by making an orthogonal transformation which diagonalizes $`(G)_{ij}`$. With this definition, Eqs. 35 and 38 yield the solvent density and counter ion charge density

$`\varphi _{s;i}`$ $`=`$ $`z_s\left(e^{w_s}\right)_{i,1},`$ (67)   
$`\rho _{c;i}`$ $`=`$ $`z_c\left(e^{u_c}\right)_{i,1},`$ (68)   
$`=`$ $`z_c\left(e^{(uu_s)}\right)_{i,1}.`$ (69)   
To obtain the remaining densities, we need the end-point distribution function. From Eq. 46 we obtain

$`{\displaystyle \frac{dq_i(s)}{ds}}`$ $`=`$ $`{\displaystyle \underset{j}{}}[A_{ij}+(w_{h,eff})_{ij}\delta (s1/2)]q_j(s),`$ (71)   
$`A_{ij}`$ $`=`$ $`{\displaystyle \frac{2R_g^2}{D^2}}\lambda _i\delta _{ij}+(w_t)_{ij},`$ (72)   
with initial condition $`q_i(0)=\delta _{i,1}`$. The solution of this equation is

$`q_i(s)`$ $`=`$ $`\left(e^{As}\right)_{i,1},\mathrm{if}s<1/2`$ (73)   
$`=`$ $`{\displaystyle \underset{j}{}}\left(e^{w_{h,eff}}\right)_{ij}\left(e^{A/2}\right)_{j,1},s=1/2`$ (74)   
$`=`$ $`{\displaystyle \underset{j,k}{}}\left(e^{A(s1/2)}\right)_{i,j}\left(e^{w_{h,eff}}\right)_{jk}\left(e^{A/2}\right)_{k,1},s>1/2.`$ (75)   
From this, the remaining densities follow from Eqs. 49, 50, and 51:

$`\varphi _{h;i}`$ $`=`$ $`{\displaystyle \underset{jkl}{}}\left(e^{w_{h,eff}}\right)_{ij}\mathrm{\Gamma }_{jkl}q_k({\displaystyle \frac{1}{2}})q_l({\displaystyle \frac{1}{2}}),`$ (76)   
$`\varphi _{t;i}`$ $`=`$ $`{\displaystyle _0^1}𝑑s{\displaystyle \underset{jk}{}}\mathrm{\Gamma }_{ijk}q_j(s)q_k(1s),`$ (77)   
$`\rho _{h;i}`$ $`=`$ $`{\displaystyle \underset{j}{}}\left({\displaystyle \frac{pe^{u_h}}{1+p(e^{u_h}1)}}\right)_{i,j}\varphi _{h;j},`$ (78)   
with

$$\mathrm{\Gamma }_{ijk}\frac{1}{V}f_i(𝐫)f_j(𝐫)f_k(𝐫).$$
(79)   
The mean-field free energy, Eq. 40, takes the form

$$\mathrm{\Omega }_{mf}=\frac{kTV}{v_h}(\varphi _{t;1}+\varphi _{s;1}+\rho _{c;1})+E,$$
(80)   
with the mean-field energy being given by

$$E=\frac{kTV}{v_h}\underset{i}{}[2\chi N(\varphi _{h;i}+\varphi _{s;i})\varphi _{t;i}+\frac{1}{2}(\rho _{h;i}+\rho _{c;i})u_i\lambda (\rho _{c;i}\rho _{h;i})\varphi _{s;i}].$$
(81)   
We have expressed the Coulomb energy as a product of the charge densities and electrostatic potential. Note that this free energy still depends parametrically on $`D`$, the length scale of the phase, so that the value of $`D`$ which minimizes it must be determined. Once this is done, we compare the free energies obtained for phases of different space-group symmetry, and thereby determine the phase diagram of our model lipid system.

The infinite set of self-consistent equations, Eqs. 61 to 64 must be truncated in order to be solved numerically. We have employed up to 50 basis functions. This truncation is sufficient to ensure, for $`T^{}>0.03`$ and $`1/(1+\gamma _t)<0.66`$, an accuracy of $`10^4`$ in the free energy $`v_h\mathrm{\Omega }_{mf}/kTV`$. As noted, one must also determine the length scale which minimizes the free energy. This is usually straightforward as there is a single well-defined minimum for a phase of given symmetry at given thermodynamic parameters; temperature, and chemical potentials. Were there more than one minimum, this would reflect a tendency for the system to phase separate into two phases with the same non-trivial space-group symmetry, an extremely unusual occurrence. The single minimum that one finds normally is sharp, that is, the free energy varies rather rapidly with $`D`$. Only in cases in which phases are greatly swollen is the minimum extremely shallow and difficult to locate.

## V Results: The Neutral Lipid

We first apply our method to a neutral lipid. We show here the phase behavior of the neutral lipid, in the absence and in the presence of solvent. Figure 1 shows the phase diagram of the pure lipid as a function of the dimensionless temperature $`T^{}`$, and the architecture of the lipid. The latter is characterized by the single parameter $`1/(1+\gamma _t)`$ which is the relative volume of the headgroup to that of the entire lipid. It is analogous to, but not the same as, the single parameter used by Israelachvili to characterize the geometry of lipids (Israelachvili, 1985). Shown are the lamellar phase, $`L_\alpha `$, the normal and inverted hexagonal phases, $`H_I`$ and $`H_{II}`$, and the normal and inverted body-centered cubic phases $`bcc_I`$ and $`bcc_{II}`$. The occurrence of bicontinuous phases will be discussed in a later paper.

One sees that the phase behavior is reasonable and in accord with packing considerations; as the headgroup increases in volume, the system passes through a series of phases from the inverted ones with most curvature, through the lamellar phase, to the normal ones of most curvature. We also note that as the temperature increases, the lamellar phase becomes unstable to one or the other of the hexagonal phases. This is very clear for the $`H_I`$ phase, less so for the $`H_{II}`$ phase, but is true for the latter over much of the temperature range as the slope of the $`L_\alpha /H_{II}`$ phase boundary is slightly positive below $`T^{}`$ of about 0.045. This transition with temperature is well known in the analogous case of diblock copolymers (Leibler, 1980). In order to model a lipid which, like phosphatidylserine, adopts the $`H_{II}`$ configuration when essentially neutral (Cullis et al. 1985, Bezrukov et al., 1999), we have chosen $`1/(1+\gamma _t)=0.24`$ in our subsequent studies. For comparison, the value appropriate for DOPE, calculated from the molecular volumes in the literature (Rand and Fuller, 1994) is $`0.254`$. For the convenience of the reader interested in carrying out similar calculations, we present in Table I the values of the first three non-trivial Fourier components of the headgroup density $`\varphi _h(𝐫)`$, the lattice parameter $`D/R_g`$, and the free energy $`\mathrm{\Omega }_{mf}v_h/kTV`$ for the $`L_\alpha `$, $`H_{II}`$ and $`bcc_{II}`$ phases at $`T^{}=0.04`$. We note, in passing, that the relative intensities of X-ray Bragg peaks can be determined directly from the Fourier components of the various densities with which they are associated.

The effect on this neutral lipid of the addition of a solvent of small volume, characterized by $`\gamma _s=v_s/v_h=0.1`$, close to the value of $`0.096`$ (Rand and Fuller, 1994; Kozlov et al., 1994) appropriate to water and a phosphatidylethanolamine headgroup, is shown in Fig. 2. There is a lamellar phase at small solvent volume fractions and low temperatures. This phase becomes unstable with respect to the $`H_{II}`$ phase which envelops it at higher temperatures. There is a large region of two-phase coexistence between the ordered lipid-rich phases and an almost pure solvent phase. These features are reasonable, and are observed in the systems of aqueous dialkyl didodecylphosphatidylethanolamine and of diacyl diarachinoylphosphatidylethanolamine (Seddon et al, 1984). Of particular interest is that we find a small temperature region of re-entrant hexagonal-lamellar-hexagonal transitions, an unusual feature which has been observed in DOPE (Gawrisch et al. 1992, Kozlov et al. 1994). As a consequence, there is an azeotrope at which the transition between lamellar and hexagonal phases occurs without a change in the concentration of water. We have used the coordinates of this point, $`T^{}=0.06`$ and $`\varphi _s=1.42`$, denoted $`T_0`$ and $`(\varphi _s)_0`$, to normalize the temperature and solvent-density axes. There is a small region of $`bcc_{II}`$ in our phase diagram. Again, the possible occurrence of bicontinuous phases will be examined in a later publication. The uncertainty in the temperature of the phase boundaries, $`\delta T/T`$ introduced by the truncation of the number of basis functions, is approximately $`2\times 10^3`$.

As the volume fraction of water is increased, we find that the period of all structures increases, as is expected. In Fig. 3, we compare experimental results on DOPE taken in the inverted hexagonal phase at a temperature $`T=22^{}C`$, just above that of the azeotrope (Rand and Fuller, 1994), to our values calculated just above the azeotrope. Knowing the molecular weight of DOPE, we convert the volume fraction of solvent, $`\gamma _s\varphi _s`$, which occurs in the calculation, to the experimental variable of weight fraction of water. The lattice parameter of the hexagonal phase in the calculation, however, is measured in units of the radius of gyration of either lipid tail. What value should be taken to model DOPE is unknown. Hence we have used in the comparison the lattice parameter $`D`$, in units of $`D_0`$, the lattice parameter at the azeotrope. The agreement is excellent.

An effective value of the radius of gyration can be defined as that value which brings agreement between the calculated and measured lattice parameters. As the former, at the azeotrope, is $`D(T_0)D_0=4.79R_g^0`$, and the latter is 58.9Å (Rand and Fuller, 1994), the equivalent radius of gyration at the temperature of the azeotrope, $`R_g^0`$, is 12.3Å for a single tail, not unreasonable when compared to the extended length of a single chain of DOPE which is approximately 26Å.

As the temperature of the system is lowered, the period of all structures increases, which is due to the lengthening of the tails as their entropy decreases. We would again like to compare our results with those of the DOPE system. In order to do so, we must address the temperature dependence of the radius of gyration, $`R_g`$, which appears in the theory, and which is not given a priori. Because the model chain is flexible, the radius of gyration is related to the mean square end-to-end distance, $`\overline{R}`$, by a numerical constant, $`R_g=\overline{R}/\sqrt{6}`$, so their dependence on temperature is the same. To compare our results to DOPE, we shall assume that the temperature dependence of the radius of gyration which appears in the calculation is the same as that given for lipid chains by the Rotational Isomeric States Model, a model which describes the properties of such chains very well (Flory, 1969; Mattice and Suter, 1994). Thus we assume

$$R_g(T)=c\left(\frac{1<\mathrm{cos}\varphi >}{1+<\mathrm{cos}\varphi >}\right)^{1/2},$$
(82)   
where the angle $`\varphi `$ takes the values $`180^{}`$ and $`\pm 70^{}`$ corresponding to $`trans`$, $`gauche^+`$ and $`gauche^{}`$ configurations, and $`c`$ is a constant. The statistical average of $`\mathrm{cos}\varphi `$ is

$$<\mathrm{cos}\varphi >=\frac{\mathrm{cos}(180^{})+\sigma \mathrm{cos}(70^{})+\sigma \mathrm{cos}(70^{})}{1+2\sigma }$$
(83)   
with $`\sigma =\mathrm{exp}(T_{rism}/T)`$ and $`T_{rism}=280.25^{}`$K. From the behavior with temperature of $`R_g(T)`$, the lattice parameter, $`D(T)`$, at any temperature can be obtained from

$$\frac{D(T)}{D(T_0)}=\frac{[D(T)/R_g(T)]}{[D(T_0)/R_g(T_0)]}\frac{R_g(T)}{R_g(T_0)}.$$
(84)   
Again, it is the factor $`D(T)/R_g(T)`$ which occurs naturally in the calculation.

A comparison of the experimentally measured (Tate and Gruner, 1989) and theoretically calculated lattice parameters vs. temperature is shown in Fig. 4. In part (a), the variation of the parameter of the $`H_{II}`$ is shown at two different lipid weight fractions. The agreement is very good. In part (b), the comparison is made of the $`H_{II}`$ and $`L_\alpha `$ parameters along the coexistence with excess water. Note that this comparison is a much more stringent test, for it requires not only that the dependence of the lattice parameters on solvent concentration and on temperature be reproduced well by the calculation, but also that the phase boundaries be given well. Considering these requirements, the agreement is rather good. It should be noted that the agreement in 4(b)does not depend on the exact temperature of the triple point, which is difficult to locate precisely, but only on the existence of stable $`H_{II}`$ and $`L_\alpha `$ phases which coexist with excess water.

As seen in Fig. 4(b), the lattice parameter of the $`H_{II}`$ phase is much larger than that of the $`L_\alpha `$ at the triple point. This is due to the coexistence with excess solvent, which swells the hexagonal cores, but which is only weakly present between the lamellae. In contrast, when the lamellar and hexagonal phases are only in two-phase coexistence with one another and there is no excess solvent, the hexagonal phase in general has a smaller lattice parameter than the coexisting lamellar phase, as shown in Fig. 5. This is due to the fact that over almost all of their coexistence region, the $`H_{II}`$ phase has a smaller volume fraction of solvent than does the $`L_\alpha `$ phase, as can be seen from the phase diagram of Fig. 2. Only over the re-entrant region, which occurs as the triple point is approached, does this balance shift. This shift in the relative size of the latter parameters through the re-entrant region is in agreement with experiment (Kozlov et al., 1994).

It is of interest to determine if any one effect can be said to drive the transition from the $`H_{II}`$ to $`L_\alpha `$ phase in the neutral system. To this end we examine the individual terms in the thermodynamic potentials per unit volume $`\mathrm{\Omega }_{mf}v_h/kTV=Ev_h/kTVS_lv_h/kV(S_sv_h/kV+\varphi _s\mathrm{ln}z_s)`$ of the $`L_\alpha `$ and $`H_{II}`$ phases as the transition is crossed by increasing the solvent fugacity, $`z_s`$, at constant temperature $`T/T_0=0.67`$. In Table II we show the contributions to the free energy per unit volume of the $`L_\alpha `$ phase and that of the $`H_{II}`$ phase coming from the interaction energy, $`Ev_h/kTV`$, the lipid tails, $`S_lv_h/kV`$, and the solvent $`S_sv_h/kV\varphi _s\mathrm{ln}z_s.`$ All contributions are evaluated at the transition itself, which occurs at $`z_s3.14`$, and are measured with respect to the free energy per unit volume of the disordered phase. We also show the difference in the contribution of each term to the free energies of each phase, and the derivative of this difference with respect to the solvent fugacity. The difference in the contribution of the entropy of the lipid tails is positive because, with the lipid architecture we have chosen, the large tail volume relative to that of the head favors the hexagonal phase. The interaction energy favors the lamellar phase, as does the solvent, presumably because the interstices of that phase are two-dimensional, while those of the inverted hexagonal phase are one-dimensional. The difference between the lipid entropy contributions decreases with increasing solvent concentration because the packing constraints in the $`H_{II}`$ phase become more severe as the size of the cores increases (Gruner, 1989). However it is apparent that neither this term, nor either of the others, changes so rapidly with solvent concentration compared to the others that any particular effect can be said to “drive” this transition.

## VI Results: The Charged Lipid

### A Coulomb Interactions Only

We now turn on the negative charge of the headgroups by varying the chemical potential of the free counter ions while enforcing charge neutrality. Increasing the density of free, positive, counter ions in our closed system is equivalent to increasing the magnitude of the negative charge density on the headgroups. It therefore corresponds to an increase in the pH of an experimental system. The charge on the headgroups is annealed, meaning that it is determined by the local value of the electrostatic potential, and therefore the headgroup charge varies with the location of that group. The parameter $`\beta ^{}`$, defined in Eq. 9, measures the strength of the Coulomb interaction. It can be written as the ratio of two lengths, $`\beta ^{}=\xi /L_l`$, where $`\xi e^2/ϵkT`$ is the Bjerrum length, and $`L_lv_h/4\pi R_G^2`$ is a length characterizing the architecture of the lipid. It is reasonable that $`\beta ^{}`$ be larger than or of order unity, and we have arbitrarily taken $`\beta ^{}=1.`$

We find for this system that the effect of increasing headgroup charge is to decrease the temperature of all transitions. As an example, we show in Fig. 6 the way the temperature, $`T^{}`$, of the transition between $`L_\alpha `$ and $`H_{II}`$ phases varies with the magnitude of the average charge density of the headgroups $`\overline{\rho }_h\rho _{h;1}`$. The range of charge density corresponds to the headgroups varying from being neutral to fully charged. The region beyond the almost vertical line at $`|\overline{\rho }_h|0.24`$ would correspond to the headgroups being charged with a probability greater than unity, and is therefore unphysical. The maximum value of $`|\overline{\rho }_h|`$ is slightly different in the two phases. Over the physical range, one sees a 10 % decrease in the transition temperature, and therefore, a stabilization of the $`H_{II}`$ phase with respect to the $`L_\alpha `$ phase. Although at first surprising, this is a reasonable result, and can be understood as follows. In the neutral system, the single interaction parameter in the system, $`\chi `$, is a measure of the tendency for the system to order, and is proportional to the inverse temperature. When the charges are turned on, the Coulomb repulsion between headgroups opposes ordering, and therefore has an effect similar to a decrease in $`\chi `$ or an increase in temperature. As we have noted earlier, increasing the temperature causes the lamellar phase to become unstable to the hexagonal phase. A complimentary view is to note that increasing the Coulomb repulsion between headgroups increases the area per headgroup, and thus the area between headgroups, an effect which, again, is equivalent to increasing the temperature (Seddon and Templer. 1995). The increase in area permits the tails to splay further out, favoring the formation of the inverted hexagonal phase. The very same result of increasing the density of charges in a system is seen in the lyotropic phases of polyelectrolytes. There one finds that the phase diagram depends only on the ratio of $`\chi `$ to the charge on the polymer, at least in the absence of fluctuations (Nyrkova et al., 1994). Hence an increase in the charge is equivalent to a decrease in $`\chi `$ or an increase in temperature, leading to an instability of the lamellar phase.

We are confident, therefore, that the transition from $`L_\alpha `$ to $`H_{II}`$ phases with increasing charge on the headgroups is the correct result of the model. However it is completely opposite to the experimental results on phosphatidylserine in water (Hope and Cullis, 1980, Bezrukov et al., 1999). Hence there must be a crucial physical mechanism which is not included in the model which includes Coulomb interactions only between objects with a net charge. One such mechanism is the short-ranged attraction between the charges of the headgroup and the dipoles of the neutral water molecules which will cause the latter to aggregate around the former. Furthermore the counterions attracted to the head group will now also be associated with waters of hydration endowing them with a non-negligible volume. The additional volume of all these waters will increase the effective volume of the headgroup, and therefore tend to stabilize the lamellar phase. This additional volume will increase with the charge of the headgroup, and could be sufficient to counteract the Coulomb repulsion tending to stabilize the hexagonal phase. This competition is investigated next.

### B  Coulomb and Short-Ranged Solvent Interactions

We now include in the system an attractive contact interaction between charges and the neutral solvent. This models the short-ranged, thermally averaged, interaction between charges and the water dipole, one which varies with separation $`r`$ as $`(kT/6)(u/e\xi )^2(\xi /r)^4\omega (r)`$, with $`u`$ the dipole moment of water and $`\xi `$ the Bjerrum length. The above expression is valid for distances such that $`r/\xi >(u/e\xi )^{1/2}.`$ For water, $`\xi 7`$ Å, and $`\omega (r)/kT4.6\times 10^4(\xi /r)^4`$. To approximate this short-ranged interaction by a contact interaction of dimensionless strength $`\lambda `$ is equivalent to employing $`\omega /kT`$ evaluated at some fixed distance. Any reasonable choice shows that $`\lambda `$ is small. We have arbitrarily chosen $`\lambda =0.1`$, which corresponds to $`\omega (r)/kT`$ evaluated at $`1.9`$Å, a distance within the regime in which the approximate expression for the charge-dipole interaction is valid.

The results for the system with the Coulomb interaction of strength $`\beta ^{}=1`$ and which includes the interaction between charges and solvent of strength $`\lambda =0.1`$ are qualitatively different from those of the charged lipids without this interaction. The temperatures of all transitions between ordered phases now increase with increasing headgroup charge, whereas previously they decreased. Shown in Fig. 7 is the temperature of the transition between $`H_{II}`$ and $`L_\alpha `$ phases as a function of the magnitude of the average charge density on the headgroups. One sees that with increasing headgroup charge, the inverted hexagonal phase becomes unstable with respect to the lamellar phase, just as in the experiments on phosphatidylserine in water (Hope and Cullis, 1980, Bezrukov et al., 1999)

In Fig. 8a we show the volume fraction profiles in the $`L_\alpha `$ phase at a value of $`z_c=0.1228`$, $`z_s=3`$, and temperature $`T^{}=0.04`$ at which the $`H_{II}L_\alpha `$ transition occurs. The position through the system is divided by the lattice parameter of the phase which, in units of the radius of gyration of the lipid, is $`D_L/R_g4.14.`$ All looks reasonable. In particular, the volume fraction of solvent within the bilayer is negligible. This is true also when the lipid is neutral. In Fig. 8b, the charge density profile of the same structure is shown. One sees that the charge on the headgroup mimics, but does not reproduce, the headgroup volume fraction. This is because the charge on the headgroup is not fixed, but varies with the local electrostatic potential. The counter ion density is fairly uniform because we have employed a single dielectric constant, that of water, throughout the system. Were we to require that the dielectric constant be position-dependent, and to take the much lower value in the tail region appropriate to it, the counter ion density there would be much reduced. In Figs. 9a and 9b we show the same quantities for the $`H_{II}`$ phase at the same value of $`z_c`$. The cut through the system is taken along the nearest-neighbor direction, and the distance is normalized to its lattice parameter $`D_H3.83R_g.`$ Again the wavelength of the $`H_{II}`$ phase is smaller than that of the $`L_\alpha `$ phase because, in two-phase coexistence, i.e. in the absence of a reservoir of excess water, the cores of the cylinders are not swollen with water and the hexagonal phase contains a smaller volume fraction of water than does the lamellar phase. One can infer from Fig. 9a that in the nearest neighbor direction, there is far more interdigitation of lipid tails than in the lamellar phase. This makes sense, as the tails must certainly stretch to fill the space between cores in the second-neighbor direction, so that interdigitation is expected in the nearest-neighbor direction.

To investigate this transition further, we show in Table III the contributions in the $`L_\alpha `$ phase and in the $`H_{II}`$ phase of the various terms in the thermodynamic potential per unit volume $`\mathrm{\Omega }_{mf}v_h/kTV=(E_1+E_2+E_3)v_h/kTVS_lv_h/kV(S_sv_h/kV+\varphi _s\mathrm{ln}z_s)(S_cv_h/kV+\rho _c\mathrm{ln}z_c)`$. Here $`E_1`$ is the hydrophilic-hydrophobic interaction proportional to $`\chi N`$, $`E_2`$ is the electrostatic interaction proportional to $`\beta ^{}`$, and $`E_3`$ is the charge-solvent interaction proportional to $`\lambda `$. These contributions are evaluated at the transition itself, and are measured from the free energy per unit volume of the disordered phase. We also show the difference between these contributions to each phase, and the derivatives of each of these differences with respect to the counter ion fugacity, $`z_c`$. There are several interesting things to note. The electrostatic energy is a relatively small contribution to the free energy of each phase, and hardly differs between them. Therefore it does not have a large effect in bringing about the transition. The contribution of the counter ions to the free energy of each phase is of the same order of magnitude as the electrostatic interaction and, like it, does not change rapidly with the counter ion fugacity. The contribution of the short-range charge-solvent interaction is small, but it changes most rapidly with the counter ion fugacity, and therefore appears to be most important in actually bringing about the transition itself.

The physical mechanism in the experimental systems appears now to be clear. The lipid with an almost neutral headgroup forms the $`H_{II}`$ phase because the volume of the headgroup is relatively small compared to that of the entire lipid. As the charge on the headgroup is turned on, it attracts an increasing volume of waters of hydration via the attractive interaction between the charge and the dipoles of water. In addition, more counter ions, enlarged by their own waters of hydration, are attracted to the headgroup. Thus the headgroup becomes effectively larger, and drives the transition to the $`L_\alpha `$ phase. In doing so, it overcomes the effect of the Coulomb repulsion between headgroups which, in fact, opposes the transition to the $`L_\alpha `$ phase.

We have argued earlier that an increase in Coulomb repulsion was similar to an increase in temperature, and therefore its effect in favoring the $`H_{II}`$ phase could be understood from the phase diagram of Fig. 1. In a similar way, an increase in the strength of the charge-solvent interaction increases the volume fraction of the headgroup, so that one moves to the right in Fig. 1 and stabilizes the lamellar phase. The combined effect of increasing the headgroup charge is to increase the effective volume fraction of the headgroup and effective temperature. The phase which is ultimately stabilized results from the competition of these two effects. By changing the strength of the Coulomb interaction, we have verified that one can alter the effect of this competition, and not only stabilize the inverted hexagonal phase, but also bring about the inverted b.c.c. phase or even the disordered phase, scenarios which are again understandable from our discussion and the phase diagram of Fig. 1. But that the result of this competition can be the stabilization of the lamellar phase at the expense of the inverted hexagonal phase has been demonstrated by our model calculation, as well as by experiment.

We gratefully acknowledge useful communications with Drs. Pieter Cullis, John Seddon, and Sol Gruner. This work was supported in part by the National Science Foundation under grant numbers DMR9531161 and DMR9876864. One of us, (MS), would like to thank the CEA, Saclay, for their gracious hospitality while this paper was being written.

REFERENCES

Ben-Shaul, A., I. Szleifer, and W.M. Gelbart. 1985. Chain organization and thermodynamics in micelles and bilayers: I theory. J. Chem. Phys. 83:3597-3611.

Bezrukov, S.M., R.P. Rand, I. Vodyanoy, and V.A. Parsegian. 1999. Lipid packing stress and polypeptide aggregation: alamethicin channel probed by proton titration of lipid charge. Faraday Discuss. 111:000-000

Borukhov, I., D. Andelman, and H. Orland. 1998. Random polyelectrolytes and polyampholytes in solution. Eur. Phys. J. B 5:869-880.

Cullis, P.R., M.J. Hope, B. de Kruijff, A.J. Verkleij, and C.P.S. Tilcock. 1985. Structural properties and functional roles of phospholipids in biological membranes. in Phospholipid and Cellular Regulation, vol. 1., J.F. Kuo (ed.). CRC Press, Boca Raton, Florida 1-59.

Fattal, D.R. and A. Ben-Shaul. 1994. Mean-field calculations of chain packing and conformational statistics in lipid bilayers: comparison with experiments and molecular dynamics studies. Biophys. J. 67:983-995.

Flory, P.J. 1969 Statistical mechanics of chain molecules, Wiley Interscience, New York.

Gawrisch, K., V.A Parsegian, D.A. Hadjuk, M.W. Tate, S.M. Gruner, N.L. Fuller, and R.P. Rand. 1992. Energetics of a hexagonal-lamellar-hexagonal-phase transition sequence in dioleoylphosphatidylethanolamine membranes. Biochemistry 31:2856-2864.

Gruen, D.W.R. 1981. A statistical mechanical model of the lipid bilayer above its phase transition Biochim. Biophys. Acta 595:161-183

Gruen, D.W.R. 1985. A model for the chains in amphiphillic aggregates.1. Comparison with a molecular dynamic simulation of a bilayer J. Phys. Chem. 89:146-153.

Gruner, S.M. 1989. Stability of lyotropic phases with curved surfaces. J. Phys. Chem., 93:7652-7570.

Helfrich, W. 1973. Elastic properties of lipid bilayers: theory and possible experiments. Z. Naturforsch. 28C:693-703.

Heller, H., M. Schaefer, and K. Schulten. 1993. Molecular dynamics simulation of a bilayer of 200 lipids in the gel and in the liquid-crystal phases. J. Phys. Chem. 97:8343-8360.

Henry, N.F.M. and K. Lonsdale 1969 (eds.) International Tables for X-Ray Crystallography, Kynoch, Birmingham.

Hope, M.J., and P.R. Cullis. 1980. Effects of divalent cations and pH on phosphatidylserine model membranes: A <sup>31</sup>P NMR study. Biochemical and Biophysical Research Communications 92:846-852.

Israelachvili, J.N. 1985. Intermolecular and Surface Forces, Academic Press, San Diego.

Jacobs, R.E. and S.H. White. 1989. The nature of the hydrophobic bonding of small peptides at the bilayer interface: implications for the insertion of transbilayer helices. Biochemistry 28:3421-3437

Kirk, G.L., S.M. Gruner, and D.L. Stein. 1984 A thermodynamic model of the lamellar to inverse hexagonal phase transition of lipid membrane-water system. Biochemistry. 23:1093-1102.

Kirk, G.L. and S.M. Gruner. 1985. Lyotropic effects of alkanes and headgroup composition on the $`L_\alpha H_{II}`$ lipid liquid crystal phase transition: hydrocarbon packing versus intrinsic curvature. J. Physique 46:761-769.

Kozlov, M.M., S. Leiken, and R.P. Rand. 1994. Bending, hydration and interstitial energies quantitatively account for the hexagonal-lamellar-hexagonal re-entrant phase transition in dioleoylphosphatidylethanolamine. Biophys. J 67:1603-1611.

de Kruijff, B. 1997. Lipid polymorphism and biomembrane function. Current Opinion in Chemical Biology,1:564-569.

Leermakers, F.A.M., and J.M.H.M. Scheutjens. 1988 Statistical thermodynamics of association colloids. I. Lipid bilayer membranes. J. Chem. Phys. 89:3264-3274.

Leermakers, F.A.M., J.M.H.M. Scheutjens, and J. Lyklema. 1990. Statistical thermodynamics of association colloids. IV. Inhomogeneous membrane systems. Biochim. Biophys. Acta 1024:139-151.

Leibler, L. 1980. Theory of Microphase Separation in Block Copolymers. Macromolecules 13:1602-1617.

Marcelja, S. 1974. Chain ordering in liquid crystals II. Structure of bilayer membranes Biochim. Biophys. Acta 367:165-176

Matsen, M.W., and M. Schick. 1994. Stable and unstable phases of a diblock copolymer melt. Phys. Rev. Lett. 72:2660-2663.

Matsen, M.W. 1995. Stabilizing new morphologies by blending homopolymer with block-copolymer. Phys. Rev. Lett. 74:4225-4228.

Mattice, W.L. and U.W. Suter. 1994 Conformational theory of large molecules; the rotational isomeric state model in macromolecular systems, Wiley-Interscience, New York.

Meijer, L.A., F.A.M. Leermakers, and A. Nelson. 1994. Modelling of electrolyte ions-phospholipid layers interaction Langmuir 10:1199-1206.

Müller, M. and M. Schick, 1998. Calculation of the phase behavior of lipids. Phys. Rev. E 57:6973-6978.

Nyrkova, I.A., A.R. Khokhlov,.and M. Doi. 1994. Microdomain structures in polyelectrolyte systems: calculation of the phase diagram by direct minimization of the free energy. Macromolecules 27:4220-4230.

Rand, R.P. and N.L. Fuller. 1994. Structural dimensions and their changes in a reentrant hexagonal-lamellar transition of phospholipids. Biophys. J. 66: 2127-2138.

Rand, R.P. and V.A. Parsegian. 1989. Hydration forces between phospholipid bilayers. Biochim. Biophys. Acta. 1031:1-69.

Seddon, J.M., G. Cevc, and D. Marsh. 1983. Calorimetric studies of the gel-fluid ($`L_\alpha L_\beta `$) and lamellar-inverted hexagonal ($`L_\alpha H_{II}`$) phase transitions in dialkyl- and diacylphosphatidylethanolamines. Biochemistry, 22:1280-1289.

Seddon, J.M., G. Cevc, R.D. Kaye, and D. Marsh. 1984. X-ray diffraction study of the polymorphism of hydrated diacyl- and dialkylphosphatidylethanolamines. Biochemistry 23:2634-2644.

Seddon, J.M. 1990. Structure of the inverted hexagonal ($`H_{II}`$) phase, and non-lamellar phase transitions of lipids. Biochimica and Biophysica, 1031:1-69.

Seddon, J.M. and R. H. Templer 1995. Polymorphism of lipid-water systems in Structure and dynamics of membranes, vol. 1, R. Lipowsky and E. Sackmann (eds.). Elsevier Press, Amsterdam 97-160.

Steenhuizen,L., D. Kramer, and A. Ben-Shaul. 1991. Statistical thermodynamics of molecular organization in the inverse hexagonal phase.J. Phys. Chem. 95:7477-7483.

Tate, M.W. and S.M. Gruner 1989. Temperature dependence of the structural dimensions of the inverted hexagonal ($`H_{II}`$) phase of phosphatidylethanolamine-containing membranes. Biochemistry 28:4245-4253.

Anhydrous, neutral lipid: the lattice parameter, the free energy, and the first 3 non-trivial Fourier components of $`\varphi _h(𝐫)`$ for $`L_\alpha `$, $`H_{II}`$ and $`bcc_{II}`$ phases at $`T^{}=0.04`$ and $`1/(1+\gamma _t)=0.24`$. Recall that $`\varphi _{h,1}=\varphi _{t,1}=0.24`$ for all phases.

$`\begin{array}{cccccc}& D/R_g& v_h\mathrm{\Omega }_{mf}/kTV& \varphi _{h,2}& \varphi _{h,3}\times 10^2& \varphi _{h,4}\times 10^2\\ L_\alpha \hfill & 2.921& 0.7969& 0.2272& 6.781& 1.171\\ H_{II}\hfill & 3.167& 0.7936& 0.2368& 0.476& 2.372\\ bcc_{II}\hfill & 3.400& 0.8027& 0.2204& 1.094& 4.038\end{array}`$

Neutral lipid: contributions to the free energy per unit volume and temperature in the $`L_\alpha `$ phase and in the $`H_{II}`$ phase, the difference in these contributions, and the derivative of this difference with respect to the solvent fugacity. All contributions are evaluated at the $`L_\alpha `$, $`H_{II}`$ transition occurring on the water-poor side of the azeotrope. The solvent chemical potential at the transition is $`z_s3.15`$. All contributions are measured with respect to the free energy of the disordered phase, ($`D`$). The temperature, $`T/T_0=0.67`$. $`v_hE/VkT`$, $`𝒮_lv_hS_l/Vk`$, and $`𝒮_sv_hS_s/Vk+\varphi _s\mathrm{ln}z_s`$.

$`\begin{array}{ccccc}& \hfill L_\alpha D& \hfill H_{II}D& \hfill L_\alpha H_{II}& \hfill d(L_\alpha H_{II})/dz_s\\ \hfill & \hfill 0.6923& \hfill 0.6086& \hfill 0.0837& \hfill 0.0244\\ 𝒮_l\hfill & \hfill 1.1367& \hfill 0.9123& \hfill 0.2250& \hfill 0.0413\\ 𝒮_s\hfill & \hfill 0.6666& \hfill 0.5253& \hfill 0.1413& \hfill 0.0394\end{array}`$

Charged lipid with both Coulomb and charge-solvent interactions: contributions to the free energy per unit volume and temperature in the $`L_\alpha `$ phase and in the $`H_{II}`$ phase, the difference in these contributions, and the derivative of this difference with respect to the counter ion fugacity. All contributions are evaluated at the $`L_\alpha `$, $`H_{II}`$ transition, $`z_c=0.1228`$, and are measured with respect to the free energy of the disordered phase, ($`D`$). The temperature, $`T^{}=0.04`$, $`z_s=3`$, $`_iv_hE_i/VkT`$, with $`i=1,2,3`$ being the hydrophilic-hydrophobic, Coulomb, and charge-neutral interactions respectively, $`𝒮_lv_hS_l/Vk`$, $`𝒮_sv_hS_s/Vk+\varphi _s\mathrm{ln}z_s`$, and $`𝒮_cv_hS_c/Vk+\rho _c\mathrm{ln}z_c`$,

$`\begin{array}{ccccc}& \hfill L_\alpha D& \hfill H_{II}D& \hfill L_\alpha H_{II}& \hfill d(L_\alpha H_{II})/dz_c\\ _1\hfill & \hfill 0.6904& \hfill 0.6059& \hfill 0.0845& \hfill 0.0152\\ _2\hfill & \hfill 0.0037& \hfill 0.0024& \hfill 0.0013& \hfill 0.0196\\ _3\hfill & \hfill 0.0378& \hfill 0.0268& \hfill 0.0110& \hfill 0.1082\\ 𝒮_l\hfill & \hfill 1.1373& \hfill 0.9110& \hfill 0.2263& \hfill 0.0257\\ 𝒮_s\hfill & \hfill 0.6365& \hfill 0.5034& \hfill 0.1331& \hfill 0.0215\\ 𝒮_c\hfill & \hfill 0.0030& \hfill 0.0020& \hfill 0.0010& \hfill 0.0153\end{array}`$

Phase diagram of the neutral lipid as a function of dimensionless temperature, $`T^{}1/2\chi N`$, and relative headgroup volume, $`1/(1+\gamma _t)`$. In addition to the disordered phase, $`D`$, there are normal and inverted body-centered cubic phases, $`bcc_I`$ and $`bcc_{II}`$, normal and inverted hexagonal phases, $`H_I`$ and $`H_{II}`$, and the lamellar phase $`L_\alpha `$.

Phase diagram of a neutral lipid with $`1/(1+\gamma _t)=0.24`$ in a solvent with $`\gamma _s=0.1`$ as a function of temperature, $`T/T_0`$, and fraction of solvent, $`\varphi _s/(\varphi _s)_0`$, with $`T_0`$ and $`\gamma _s(\varphi _s)_0`$ the temperature and volume fraction of solvent at the azeotrope.

Comparison of theoretically calculated and experimentally measured values of the lattice parameter $`D/D_0`$ of the $`H_{II}`$ phase at a temperature just above the azeotrope vs. weight fraction of water, $`\varphi _w^w`$. The lattice parameter at the azeotrope is denoted $`D_0`$.

Comparison of theoretically calculated and experimentally measured values of the lattice parameter, $`D(T)/D_0`$, vs. absolute temperature $`T`$ (a) for two different weight fractions of lipid, $`\varphi _l^w`$, and (b) along coexistence with excess water. The absolute temperature of the azeotrope is denoted $`T_0`$.

Lattice parameter of the $`H_{II}`$ and $`L_\alpha `$ phases along their mutual coexistence as a function of absolute temperature. The absolute temperature of the azeotrope is denoted $`T_0`$.

Transition temperature between $`L_\alpha `$ and $`H_{II}`$ phases for the same lipid in Figure 2, but now with a charged headgroup. The dimensionless strength of the Coulomb interaction is $`\beta ^{}=1`$.The transition temperature is plotted as a function of the magnitude of the average charge density on the headgroups. The solvent fugacity is fixed at $`z_s=3`$.

Transition temperature between $`L_\alpha `$ and $`H_{II}`$ phases for the same lipid as in Figure 6, but now including the interaction between charges and neutral solvent of strength $`\lambda =0.1`$. The solvent fugacity is again fixed at $`z_s=3.`$

a) Volume fraction distribution in the $`L_\alpha `$ phase of the solvent, headgroups, and tails, in the above system at a counter ion chemical potential of $`z_c=0.1228`$ corresponding to the $`L_\alpha H_{II}`$ transition. The temperature is $`T^{}=0.04`$, and the lattice parameter of the lamellar phase is $`D_L/R_g=4.14`$; b) charge densities arising from the headgroups, the counter ions, and the total charge density in the $`L_\alpha `$ phase under the same conditions as in a).

a) Volume fraction distribution in the $`H_{II}`$ phase of the solvent, headgroups, and tails, in the above system at a counter ion chemical potential of $`z_c=0.1228`$ corresponding to the $`L_\alpha H_{II}`$ transition. The temperature is $`T^{}=0.04`$, and the lattice parameter of the inverted hexagonal phase is $`D_H/R_g=3.83`$. b) charge densities arising from the headgroups, the counter ions, and the total charge density in the $`H_{II}`$ phase under the same conditions as in a).