INTRODUCTION

TOP ABSTRACT INTRODUCTION THE MODEL THEORY RESULTS REFERENCES

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, to inverted hexagonal, H_II, phases (Seddon, 1990), whereas 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., 1998). It is these factors that 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 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, 1985; Ben-Shaul et al., 1985; Fattal and Ben-Shaul, 1994). Results for the bilayer are in good agreement with molecular dynamic simulation (Tieleman et al., 1997). These methods have shown that, in a neutral, anhydrous system, the entropy of the lipid tails always favors the H_II over the L 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 that 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.

It would be extremely useful to have available a relatively simple and tractable model of lipids that was capable, at least, of describing the effect of their architecture upon their phase behavior. With this in hand, one could, inter alia, examine the various bicontinuous phases to determine their stability or metastability (Shyamsunder et al., 1988), and to explicate the reasons they facilitate the crystallization of membrane proteins (Landau and Rosenbusch, 1996). Further, one could explore mixtures of lamellar- and nonlamellar-forming lipids to determine the role that the latter play in lipid-protein interactions (Epand, 1998), membrane fusion (Markin et al., 1984; Siegel, 1993), and membrane function (Hui, 1997), all areas in which the importance of their presence has been indicated.

Toward this end, a model system of solvent and monoacyl lipid embedded in a continuous space was recently introduced. Its phase diagram was obtained by solving the mean-field theory exactly (Müller and Schick, 1998). It displayed both L 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 use a model of a lipid computationally more tractable than that used by Müller and Schick: one whose hydrocarbon tails are modeled as flexible chains rather than within the rotational isomeric states framework used earlier (Flory, 1969; Mattice and Suter, 1994). We first study the model with an uncharged headgroup. Its phase behavior, both with respect to variations in architecture and 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. We introduce counter ions into the system, include the Coulomb interaction between all charges, and also a short-ranged interaction between charges and neutral solvent, an interaction that models the thermally averaged interaction between charges and the dipole of water. As the charge on the headgroup is turned on, the L 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 that are attracted to the headgroup are also enlarged by their own waters of hydration. It is the totality of these waters that 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 the Theory section, we first 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 that determine the headgroup, tail, and solvent densities. One of these self-consistent conditions is simply the nonlinear Poisson-Boltzmann equation. We then 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 the Results section, we first 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 this subsection, 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 the next subsection, 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, and all short-ranged interactions between charges and solvent. We find that the L is indeed stabilized with respect to the H_II phase, in agreement with experiment (Hope and Cullis, 1980; Bezrukov et al., 1998).

	THE MODEL

TOP ABSTRACT INTRODUCTION THE MODEL THEORY RESULTS REFERENCES

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²/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, whereas the n_c counter ions have charge +e and negligible volume, v_c = 0. There are five dimensionless densities that totally specify the state of the system; the number density of the headgroups, _h, of the tail segments, _t, and of the solvent, _s, and the charge density of the headgroups, e_h, and of the counter ions, e_c. They can be written as

(1)

(2)

(3)

(4)

(5)

We have chosen v_h as a convenient volume to make all densities dimensionless. In the above, R_s,j is the position of the jth solvent particle, and R_c,i the position of the ith counter ion. The configuation of the lth lipid is described by a space curver r_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 lth lipid, eQ_h,l, is equal to e or 0 is p or 1  p, respectively. As we model the case in which charges can associate or dissociate 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 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, 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. This is to model the short-ranged, thermally averaged interaction between charges and the dipole of water, an attractive interaction that decreases like r⁴ and is of strength e²u²/6²kT, where u is the dipole moment of water (Israelachvili, 1985). Thus, the energy per unit volume of the system, E/V, can be written

(6)

where

(7)

is a dimensionless measure of the strength of the Coulomb interaction. The grand partition function (Matsen, 1995) of the system is

(8)

Here,  r_l denotes a functional integral over the possible configurations of the lth lipid and in which, in addition to the Boltzmann weight, the path is weighted by the factor [r_t,l(s); 0, 1], with

(9)

with  an unimportant normalization constant. The notation  Q_h,l denotes an integral over the probability distribution of the charge on the headgroup of the lth lipid. We have enforced an incompressibility constraint on the system with the aid of the delta function (1  _{h  ss  tt}), where _s = v_s/v_h, and _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 + _t).

The model is now completely defined. The solvent is specified by _s, its volume per particle relative to that of the headgroup, and the architecture of the lipid is characterized by _t. There are three interactions, hydrophobic-hydrophilic, charge-charge, and charge-solvent, whose strengths are given by , *, and , respectively. The external parameters are the temperature, conveniently specified in terms of a dimensionless temperature T*  (2N)¹, 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 section, we derive the self-consistent field theory for the model, first in real space, and then in Fourier space.

	THEORY

TOP ABSTRACT INTRODUCTION THE MODEL THEORY RESULTS REFERENCES

Real space

Evaluation of the partition function of Eq. 8 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 _h(r) which, from its definition in Eq. 1, depends on the coordinates of the headgroup, r_l(1/2). One introduces into the partition function the identity

(10)

in which _h(r) does not depend on any specific coordinates of one of the elements of the system, but is simply a function of r. The integration on W_h extends up the imaginary axis. Inserting such identities for the five densities _h, _t, _s, _h, and _c, and a similar identity for the delta function expressing the incompressibility condition, one rewrites the partition function, Eq. 8, as

(11)

where

(12)

is the partition function of a single lipid in external fields W_h, W_t, and U_h,

(13)

is the partition function of a single counter ion of unit positive charge in an external potential U_c, and

(14)

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_l  1/v_h, z_c  z_c/v_h, and z_s  z_s/v_h. The partition function, Eq. 11, can then be written in the form

(15)

with

(16)

No approximations have been made to this point. What has been accomplished is a rewriting of the partition function from a form, Eq. 8, in which all entities interact directly with one another, to a form, Eqs. 15 and 16, in which they interact indirectly with one another via fluctuating fields. Although the integrals in Eq. 15 over _h, _t, _s, P_h, P_c, and  could all be carried out, because 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 use the self-consistent field theory in which we replace the integral in Eq. 15 by its integrand evaluated at its extremum. The values of W_h, _h, etc., which satisfy the extremum conditions, will be denoted by the corresponding lower-case letters w_h, and _h, etc. The equations that determine them are six self-consistent equations for the six fields w_h, w_t, w_s, u_h, u_c, and . They are

(17)

(18)

(19)

(20)

(21)

(22)

Because the field  is easily eliminated, the six equations readily reduce to five. The simplicity of Eq. 21 reduces this, in practice, to a set of four equations. The five densities _h, _t, _s, _h, and _c are functionals of all of the above fields except , and, therefore, close the cycle of self-consistent equations:

(23)

(24)

(25)

(26)

(27)

(28)

(29)

The density _h(r) is simply the expectation value of _h(r) in the single lipid ensemble. Similar interpretations follow for the other densities. Note that one of the self-consistent equations, Eq. 20, is simply the nonlinear Poisson-Boltzmann equation, and u(r) the electric potential.

With the aid of the above equations, the mean-field free energy, _mf, which is the free energy function of Eq. 16 evaluated at the mean-field values of the densities and fields, can be put in the form

(30)

(31)

with E given by Eq. 6. The thermodynamic potential, , 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 1  p (Borukhov et al., 1998). The consequence of this averaging is that _l[w_h, w_t, u_h] of Eq. 12 becomes

(32)

where

(33)

(34)

Although this appears to introduce an unknown parameter p into the problem, the condition of charge neutrality,

(35)

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

(36)

which satisfies the equation

(37)

with initial condition

(38)

The partition function of the lipid is then

(39)

From this expression for the single-lipid partition function and Eqs. 23, 24, and 27, one obtains expressions for the local density of the lipid heads,

(40)

of the lipid tails,

(41)

and of the charge density on the lipid heads,

(42)

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. 17-22, can be taken to be

(43)

(44)

(45)

(46)

Note that we have chosen here to write the Poisson-Boltzmann equation, Eq. 20, in its local, rather than its integral form. When the four fields are known, the corresponding densities follow from Eqs. 26, 29, 40, 41, and 42.

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 that makes straightforward their solution for a phase of arbitrary space-group symmetry (Matsen and Schick, 1994).

Fourier space

We note that the fields, densities, and the end point distribution function depend only on one coordinate r. 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(r), i = 1, 2, 3, ... , each of which have the desired space group symmetry; e.g.,

(47)

(48)

Furthermore, we choose the f_i(r) to be eigenfunctions of the Laplacian

(49)

where D is a length scale for the phase. The functions for the lamellar phase are clear. They can be taken to be

(50)

(51)

Expressions for the unnormalized basis functions for other space-group symmetries can be found in X-ray tables (Henry and Lonsdale, 1969) because 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

(52)

(53)

(54)

(55)

To obtain the partition functions and densities, we proceed as follows. For any function G(r), we can define a symmetric matrix,

(56)

Note that (G)_1i = (G)_i1 = G_i, the coefficient of f_i(r) in the expansion of G(r). Matrices corresponding to functions of G(r), such as

(57)

are evaluated by making an orthogonal transformation, which diagonalizes (G)_ij. With this definition, Eqs. 26 and 29 yield the solvent density and counter ion charge density,

(58)

(59)

(60)

To obtain the remaining densities, we need the end-point distribution function. From Eq. 37, we obtain

(61)

(62)

with initial condition q_i(0) = _i,1. The solution of this equation is

(63)

From this, the remaining densities follow from Eqs. 40, 41, and 42:

(64)

(65)

(66)

with

(67)

The mean-field free energy, Eq. 31, takes the form

(68)

with the mean-field energy being given by

(69)

We have expressed the Coulomb energy as a product of the charge densities and electrostatic potential. Note that