INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORETICAL DEVELOPMENTS RESULTS AND DISCUSSION CONCLUSION REFERENCES

The association between transmembrane -helices is a process of basic importance for understanding protein structure and stability as well as protein-protein interactions in biological membranes. The interactions that drive such an association are highly specific in some cases, but relatively nonspecific in others (Lemmon et al., 1992; Lemmon and Engelman, 1994). In particular, nonspecific protein-lipid interactions arise from perturbations of lipid structure by the proteins themselves. The best documented nonspecific protein-lipid interactions are hydrophobic mismatch interactions between the hydrophobic length of the proteins and the hydrocarbon thickness of the surrounding lipids (Mouritsen and Bloom, 1993; Nielsen et al., 1998; Killian, 1998; May and Ben-Shaul, 1999; Harroun et al., 1999a,b), but there are also nonspecific lipid-packing effects due to hydrophobic interactions between proteins and lipids, which will be discussed below. Experimental investigations of nonspecific interactions of this kind are difficult to perform because of the length and time scales involved, although some experimental results have been obtained using spin-label electron paramagnetic resonance (Marsh and Horváth, 1998) and solid-state nuclear magnetic resonance (Watts, 1998). However, such interactions are quite amenable to a theoretical analysis, as discussed by Gil et al. (1998). The motivation of this paper is therefore to study nonspecific lipid-mediated protein-protein interactions, using a theoretical analysis involving integral equation formalisms for liquids.

Theoretical investigations of nonspecific lipid-mediated protein-protein interactions other than hydrophobic mismatch were initiated by Marcelja (1976), who proposed a mean field model of a lipid bilayer based on order parameters related to lipid chain conformational states. His theory described the effects resulting from a nonspecific interaction between integral membrane proteins and the surrounding lipids. The model assumed that the most important change in lipid structure was restricted to the annulus of those lipid chains that are in direct contact with the protein. However, at temperatures above the main gel-liquid crystal phase transition, it was found that the disturbance caused by the protein inclusion extended to its second or the third neighboring chains. Marcelja showed that the change in lipid order gives rise to an indirect lipid-mediated interaction between membrane integral proteins, leading to a monotonically attractive potential between two proteins embedded in a lipid bilayer in the liquid crystalline phase with a free energy well of 1-3 k_BT at contact. A similar idea was developed by Sabra et al. (1998), on the basis of a lattice model where two different lipid species were included: annular lipids interacting strongly with the proteins, and neutral lipids interacting weakly with the proteins. The Monte Carlo simulations performed with this model showed that lipid-mediated two-dimensional (2D) arrays of membrane proteins only form when there are annular lipids present in the bilayer. Such arrays have been observed experimentally (Kuhlbrandt, 1992).

Different mean-field theories have been proposed by Schroder (1977), Owicki et al. (1978), Owicki and McConnell (1979), and Pearson et al. (1984). In these theories, the state of the lipid bilayer is characterized by several spatially inhomogeneous coarse-grained "order parameters" that are directly related to fluctuations in the lateral density of the lipid chains. The equation for the spatial variation of the order parameter field is derived from a Landau-de Gennes free energy functional with a limited gradient expansion (de Gennes, 1974). Interaction strengths and correlation lengths for the pure membrane are described in terms of phenomenological parameters. Similar ideas were used in more recent theoretical work to incorporate the influence of membrane stretching, bending moduli, and spontaneous curvature (Goulian et al., 1993; Kralchevsky et al., 1995; Aranda-Espinoza et al., 1996; Kim et al., 1998).

Mean-field theoretical treatments of this type generally predict an attractive lipid-mediated protein-protein interaction (except for that of Kim et al., 1998) that decays monotonically as a function of the distance between two proteins embedded in a lipid bilayer in the liquid crystalline phase, although the strength of the interaction was often found to depend on the bilayer phase (Goulian et al., 1993; Aranda-Espinoza et al., 1996). Given that the influence of a protein inclusion is to perturb the natural liquid crystalline state of the membrane, the indirect interaction between membrane proteins was obtained under the assumption that fluctuations are suppressed in the vicinity of the proteins. The lipid-mediated protein-protein interaction is then caused by the overlap of the lipid annuli surrounding the proteins, and its range was found to increase with increasing correlation length. Typically, the magnitude of this interaction is on the order of 1-3 k_BT at protein contact.

Sintes and Baumgärtner (1997a,b) examined the problem of lipid-mediated protein-protein interactions by using Monte Carlo computer simulations based on a simple model of the lipid bilayer. The model represented the bilayer with 2 × 500 lipid molecules, where each molecule was modeled by a flexible chain composed of five monomers. The proteins were modeled by two hard transbilayer cylinders. The simulations gave a depletion-induced attraction between proteins lying closer than one lipid diameter and a fluctuation-induced attraction for larger interprotein displacements that has a correlation length of about three lipid diameters. However, the main limitation of their approach was again that the description of the bilayer lipids was necessarily simplified to decrease the computational cost.

One important limitation of these earlier studies is the significant number of approximations that has to be introduced to construct tractable analytical theories. The lipid bilayer, a complex macromolecular assembly of amphiphilic phospholipid molecules, is thus described in terms of a necessarily limited phenomenological free energy functional. Such approximations contrast with the microscopic view of the structure and dynamics of biological membranes provided by molecular dynamics (MD) simulations based on detailed atomic models. In the last few years, studies of pure lipid bilayers (Egberts and Berendsen, 1988; Feller et al., 1997) and protein-membrane systems (Woolf and Roux, 1996; Shen et al., 1997; Tieleman and Berendsen, 1998) have demonstrated the feasibility and success of such detailed simulations (see also Merz and Roux, 1996, and references therein). In principle, free energy MD simulations and perturbation techniques could be used to calculate the solvation free energy of inclusions (Postma et al., 1982). However, such calculations are computationally prohibitive and cannot be used to address such aspects of membrane structure at the present time. For example, it is actually not feasible to examine long-range protein-protein interactions embedded in a lipid bilayer with MD simulations because of the very long time scales involved in the relaxation of the lipids. Progress on general questions concerning protein-protein interactions therefore requires alternative approaches.

Recently, an approach based on statistical mechanical theories involving integral equations was developed for the study of liquids (Hansen and McDonald, 1986; Chandler et al., 1986). It was used to examine the influence of lipid chains on protein-protein interactions (Lagüe et al., 1998). The theory was derived as a hypernetted chain (HNC) integral equation projected onto the 2D space of the lipid bilayer plane. The exact lateral density-density response function of the hydrocarbon core, calculated from the configurations of an MD simulation of a lipid bilayer (Feller et al., 1997), is used as an input to this theory. Response functions of this type are closely related to the bilayer structure factor, which can be extracted from low-angle x-ray or neutron scattering measurements (Hansen and McDonald, 1986). Such a theory, constructed on the basis of a density susceptibility response function used as an input to derive environment-mediated potential of mean force between impurities, is in the spirit of the Pratt-Chandler theory of the hydrophobic effect (Pratt and Chandler, 1977). To illustrate the approach, the lateral perturbations on the dipalmitoylphosphatidylcholine (DPPC) bilayer structure as well as the lipid-mediated protein-protein interaction were calculated for a 5-Å-radius cylinder corresponding to a -helical polyalanine molecule. This theory offers an intermediate approach, combining aspects of both mean-field theories and fully detailed atomic simulations.

In this paper we extend the HNC equation by exploring approximations based on the Percus-Yevick (PY) equation (Hansen and McDonald, 1986) to calculate several quantities related to lipid-mediated protein-protein interactions. We used the two approximations for comparison to see if any essential differences could be uncovered in the lipid density around the embedded proteins and in the lipid-mediated forces acting between two embedded proteins. In the next section we present the theoretical formulation of the problem, followed by the results for a DPPC bilayer. The paper is concluded with a brief summary and a discussion of future work.

	THEORETICAL DEVELOPMENTS

TOP ABSTRACT INTRODUCTION THEORETICAL DEVELOPMENTS RESULTS AND DISCUSSION CONCLUSION REFERENCES

Integral equation theories

The method employed here was originally presented by Lagüe et al. (1998). However, because several extensions have been added, we briefly review the theory. Isolated protein inclusions embedded in a uniform lipid bilayer in the liquid crystalline state are considered. It is assumed that the dominant perturbation affects only the lateral positions of the lipids. For the sake of simplicity, it is assumed that the protein inclusions are hard repulsive cylinders of radius  that interact only with the hydrocarbon chains, whereas the polar headgroups are not directly affected by the protein. The proteins are modeled as hard, vertical, straight cylinders. The carbons belonging to the lipid acyl chains interact with the protein inclusions via a repulsive potential, U(r)  U(x, y). For n protein inclusions, the total perturbation potential is

(1)

where r_i  (x_i, y_i) is the position of the ith inclusion in the membrane plane.

For the development of the theoretical analysis, we use the HNC and the PY approximations. These equations allow us to calculate the average density of the carbon atoms projected in the 2D membrane plane (r; r₁, ... , r_n), where r  (x, y). (The explicit dependence of the lipid density upon the position of the n protein inclusions will be omitted in the following for the sake of clarity.) We begin by writing an expression for the free energy density functional for nonuniform liquids in the HNC approximation (Hansen and McDonald, 1986; Chandler et al., 1986),

(2)

where is the density of the hydrocarbon chains in the uniform 2D membrane plane and (r) = (r)   is the deviation from the uniform density (in other words, (r) represents the perturbation of the lipid density by the protein inclusion). The lipid-lipid direct correlation function C_m(|r  r'|) is defined in terms of _m(r), the equilibrium carbon-carbon density susceptibility of the uniform unperturbed membrane,

(3)

Here _m(r) is a response function related to the density fluctuations of carbon pairs in the unperturbed membrane (_m does not depend upon the position of the protein inclusions; see below), and _m¹(r) is the functional inverse of the density response function _m(r). According to the free energy variational principle (Chandler et al., 1986), the average density is obtained by minimization of the functional  with respect to the function (r). This leads to the 2D-HNC integral equation

(4)

which must be solved self-consistently. The integral equation can be rewritten in a form more suitable for numerical calculations as a pair of coupled equations,

(5)

and

(6)

where h(r)  (r)/ is the protein-lipid correlation function, and c(r) is the protein-lipid direct correlation function. Equation 6 is the well-known Ornstein-Zernike (OZ) equation (Hansen and McDonald, 1986) for an isolated impurity in an infinite bulk system. A 2D-PY equation is obtained by linearizing the exponential on the right-hard side of Eq. 5:

(7)

This equation can also be solved self-consistently when coupled with Eq. 6.

Lipid-mediated interaction energy between two protein inclusions

It is possible to obtain the lipid-mediated potential of mean force (PMF) between two protein inclusions directly from the OZ equation (Hansen and McDonald, 1986). The PMF between two cylindrical protein inclusions is then

(8)

where c(r) is the protein-lipid direct correlation function for a single isolated protein inclusion. The double convolution is calculated by using a Fourier transform. We call Eq. 8 the OZ route. This method assumes that the concentration of protein inclusions in the membrane is infinitely small.

Alternatively, it is possible to compute the PMF from the excess Helmholtz free energy of the system by solving for the lipid density around the two protein inclusions located explicitly at r₁ and r₂. The excess Helmholtz free energy  due to the two cylinders is obtained by substituting the self-consistent solution to the 2D-HNC integral equation (Eq. 4) into the free energy functional of Eq. 2. This equation, together with the OZ equation of Eq. 6, leads to the closed-form expression for the excess free energy of an arbitrary perturbation (see also Morita and Hiroike, 1960),

(9)

where the notation h(r; r₁, r₂) and c(r; r₁, r₂) indicates that the correlation functions depend parametrically upon r₁ and r₂. The PMF is the difference in free energy between the system of two cylinders at r₁ and r₂ minus the free energy when they are infinitely separated, (r) = (r)  2(). We call Eq. 9 the  route.

Finally, it is possible to compute the PMF from the reversible work needed to bring two protein inclusions from an infinite separation to a distance r (Kirkwood, 1935),

(10)

where F(r) is the mean radial force directed along the cylinder-cylinder axis acting on one cylinder. The system with two cylinders is shown schematically in Fig. 1. The lipid-mediated mean force acting on cylinder 1 at r₁, in the presence of a second cylinder 2 at r₂, is

(11)

where U(r₁, r₂) is, as above, the total perturbating potential, and ₁₂ is a unit vector oriented along the line joining the centers of the two proteins. The mean force is given by

(12)

where (r = , ) is the lipid density at the outer surface of the protein inclusion for a given angle . The angular integration is taken over the circumference of the protein inclusion (dr  r d dr). The identity

(13)

following from the assumption of infinitely repulsive cylinders of radius , was used. Equation 12 shows that the net average lipid-mediated force acting on two protein inclusions arises from the asymmetry in hydrocarbon density around the protein inclusions. By symmetry, there is no net force acting on a single isolated cylinder. We call Eq. 12 the F route. The F route can be used with both the 2D-HNC integral equation (Eq. 4) and the 2D-PY integral equation (Eq. 7).

View larger version (7K): [in this window] [in a new window]   FIGURE 1   Diagram showing how the mean forces between two protein inclusions, modeled as cylinders, were calculated. For distances greater than the protein diameter, there is a force arising from lipid density asymmetry around protein inclusions, as stated in the text. The resultant force is obtained by integrating the average density of hydrocarbon chains (r = , ) at the positive side of the protein radius around the protein inclusion, as given by Eq. 12.

The calculated PMF corresponds to the lipid-mediated interaction between two helices modeled as straight vertical cylinders. These cylinders have no tilting with respect to the membrane normal. This is clearly an approximation because it is well known that transmembrane helices are often associated with some angle to optimize the packing of their side chains; for example, the angle between the helices in the glycophorin dimer is around 40° (MacKenzie et al., 1997). Such a phenomenon, which arises from local helix-helix interactions, cannot be taken into account in the current theory, which is built upon a two-dimensional projection of all interaction and correlation onto the membrane plane. Extensions to a more complex theory to account for the three-dimensional spatial organization of the membrane are in progress to deal with such questions.

Pair correlation function of the unperturbed membrane

A central quantity in the present theory is the response function of the uniform unperturbed membrane, _m(r), defined in Eq. 3. Functions such as _m(r) play a central role in the response of the average structure of an equilibrium system to a small perturbation (Hansen and McDonald, 1986). The response function _m is related to lipid density-density fluctuations of carbon pairs in the unperturbed membrane at equilibrium,

(14)

where (r) is the density of the ensemble of carbon atoms comprising the lipid chains,

(15)

Here i is the index of the lipid molecules, and , which goes from 1 to n, is the index of the carbon atom along the lipid chains. In the uniform unperturbed system, the average of (r) is the average carbon density per unit area,

(16)

The average density is equal to 2 × n × the surface density of lipid molecule per leaflet, where n is the number of carbon atoms in the hydrophobic moiety of one DPPC molecule (the factor 2 appears because of the upper and lower leaflets of the bilayer).

Density-density fluctuations of carbon pairs can also be expressed in terms of the radial intramolecular and intermolecular pair correlation functions of the pure unperturbed membrane,

(17)

where the function S_m(r  r') represents the carbon-carbon intramolecular pair correlation within a given lipid molecule (i = j), while the function H_m(r  r') represents the carbon-carbon intermolecular pair correlation between distinct lipids (i  j), as displayed in Fig. 2. By symmetry, the pair correlation functions depend only on the distance r projected in the x, y plane, with r = .

View larger version (46K): [in this window] [in a new window]   FIGURE 2   (A) Snapshot of the lipid bilayer taken from the trajectory used to compute the pair correlation functions of the unperturbed membrane (Eq. 14). (B) A single lipid molecule, taken from the snapshot in A, with the distance d (with r = r + d) used to compute the intramolecular correlation function Sm(r) (Eq. 18). (C) Two lipid molecules, taken from the snapshot in A, with the distance d (again with r = r + d) used to compute the extramolecular correlation function Hm(r) (Eq. 19). In both B and C, only atoms represented by filled circles were used to calculate the pair correlation function of the unperturbed membrane.

In the present study, the pair correlation functions were calculated from 500 configurations generated by molecular dynamics simulations of a detailed atomic model of a pure DPPC bilayer at 323.15 K performed by Feller et al. (1997). The function S_m(r) was calculated as

(18)

where N^intra(r; r + r) is the total number of carbons from a given lipid found within the 2D annulus going from r to r + r centered around carbon , and a(r; r + r) = [(r + r)²  r²] is the area of the annulus. The function H_m(r) was calculated as

(19)

where N^inter(r; r + r) is the number of carbons from the other lipids found within the 2D annulus going from r to r + r centered around carbon  of a lipid molecule. At large values of r, intermolecular correlations vanish and the function H_m(r)  0.

All of the heavy atoms from the acyl chains and the glycerol backbone were counted in the calculation of the pair correlation function for a total of 39 particles per lipid. The polar headgroup was not included. Because the average cross-sectional area is 62.9 Å² per DPPC (Feller et al., 1997), the average carbon density per unit area is equal to 1.24 Å².

Computational details

The 2D-HNC equations (Eqs. 5 and 6) and the 2D-PY equations (Eqs. 7 and 6) were examined for two different systems. The case of a single isolated protein inclusion modeled as a hard repulsive cylinder was first studied. For this case, the PMF was computed using Eq. 8 when the 2D-HNC closure was used. Second, two identical cylindrical protein inclusions were examined at various separations. For this case, the lipid-mediated protein-protein free energy was calculated using Eq. 9, and the lipid-mediated mean force between the two proteins was calculated using Eq. 12 with the 2D-HNC closure, whereas only the lipid-mediated force between the two proteins was calculated using the same equation but with the 2D-PY closure. Three radii for the hard cylinder were chosen: a small radius of 2.5 Å, a medium radius of 5 Å, and a larger radius of 9 Å. The small cylinder corresponds closely to an aliphatic chain, the medium cylinder corresponds closely to a polyalanine -helix, and the larger cylinder corresponds to a small protein such as the gramicidin channel (Woolf and Roux, 1996).

The 2D-HNC equations (Eqs. 5 and 6) and the 2D-PY equations (Eqs. 7 and 6) were solved numerically by a method used previously to solve HNC integral equations (Beglov and Roux, 1995, 1997). It involves a mapping of all functions onto a 2D discrete grid, e.g., u(x, y)  u(i, j), h(x, y)  h(i, j), c(x, y)  c(i, j), and _m(x, y)  _m(i, j). Two grid dimensions were used, depending on the number of proteins inserted into the system. A discrete grid for N = 1024 × 1024 with a spacing d of 0.12 Å was used with single protein systems, and N = 2048 × 1024 with the same d spacing was used for two protein systems. The 2D convolution in Eq. 6 was calculated using a numerical 2D fast Fourier transform (FFT) procedure called FFTW (Frigo and Johnson, 1998). The convolution was calculated directly, without zero padding. This corresponds to a periodic system in the x and y directions. An iterative scheme with simple mixing was used to solve 2D-HNC and 2D-PY closures self-consistently. In this scheme, the mth iteration is obtained from

(20)

for the 2D-HNC closure, and the mth iteration is obtained from

(21)

for the 2D-PY closure. Approximately 50 iterations were necessary for convergence. The average force was calculated from 1024 points around a protein inclusion, and these points were obtained by linear interpolation from grid points. The numerical solution to the integral equation took less than 5 min to obtain on a 400 MHz Intel Pentium II.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THEORETICAL DEVELOPMENTS RESULTS AND DISCUSSION CONCLUSION REFERENCES

Lipid structure and response function

Fig. 2 A shows a snapshot of the lipid bilayer taken from the trajectory of molecular dynamics simulations of Feller et al. (1997), which was used to compute the pair correlation function of the unperturbed membrane. The membrane was in the liquid crystalline phase, in which the average overall orientation of the lipid chains is perpendicular to the bilayer plane. Fig. 3 gives the results for the carbon-carbon distribution function, _m(r), extracted from molecular dynamics simulations and used as input in the integral equation. The carbon-carbon intramolecular correlation function, S_m(r), involving carbons from the same lipid molecule as schematically shown in Fig. 2 B, has a large peak up to 3 Å and then exhibits a slow decay over a distance of 10-15 Å. The short-range contribution to the intramolecular correlation arises mainly from nearest-neighbor carbons along the acyl chains (i.e., carbon i with carbons i  1 and i + 1), but there are also contributions from intermediate range (second neighbor) and long-range correlations along the lipid chain. These include correlations between the last carbon atom of the aliphatic chain and the glycerol oxygen at the beginning of the aliphatic chain. The dominant peak is thus an indication of the significant amount of short-range and long-range order in the lipid chains perpendicular to the plane of the bilayer. The long-range contribution to the intramolecular correlation function, which extends to 15 Å, is due to carbons located in different acyl chains of a single lipid molecule.

View larger version (16K): [in this window] [in a new window]   FIGURE 3   Carbon-carbon density-density response function m(r) (solid line) as extracted from molecular dynamics simulations and used as input in the integral equation. The intramolecular Sm(r) (dashed line) and intermolecular, Hm(r) (dash-dotted line), pair correlation functions extracted from the molecular dynamics simulations of Feller et al. (1997) are shown. (to match the dimensions of Å2 of Sm(r), the function m/ and the function Hm(r) ×  are also shown).

The carbon-carbon intermolecular correlation function H_m(r), involving carbons from two different lipid molecules, as shown schematically in Fig. 2 C, has a strong negative contribution at short distances due to the lipid-lipid core repulsion. It is interesting to note that for r = 0 the intermolecular correlation function has a value higher than the corresponding value for a real liquid, i.e., a value of 0.44 as compared to a real liquid value of 1.0. This comes from the fact that overlap of carbons is possible when the 3D configuration of aliphatic chains is projected in 2D onto the membrane plane. The first peak around 5 Å is in good agreement with the carbon-carbon intermolecular pair correlation function of butane (Tobias et al., 1997), whereas a characteristic double peak appears in the region between 4 Å and 6 Å, which includes intramolecular correlations between terminal methyl groups for molecules in the trans conformation as well as intermolecular correlations. This result suggests a 5-Å distance between carbon atoms of two different aliphatic chains, corresponding to a radius of 2.5 Å for a single aliphatic chain. It is of interest to note that the negative contributions are much less important at distances greater than the lipid diameter. As a result, the response function _m(r) exhibits a strong peak for distances up to 3 Å, arising from the intramolecular correlations, and a second strong peak at a distance of 4-5 Å, arising from the intermolecular correlations. The response function decays in an oscillatory manner, with small positive peaks appearing around 9 and 14 Å.

Perturbed lipid density around protein inclusion

The observed structure of the correlation functions of a pure DPPC bilayer shown in Fig. 3 suggests that the lateral response to perturbations may be quite complex. To assess the response of the membrane quantitatively, the average density around three protein inclusions modeled as hard cylinders was examined: a small cylinder, corresponding to an aliphatic chain, has a radius of 2.5 Å; a medium cylinder, with a radius of 5 Å, corresponds to an -helical polyalanine molecule; and a large cylinder of 9 Å radius, representing a small peptide such as the gramicidin channel (Woolf and Roux, 1996).

The radial average lipid density around the cylinders as calculated using the HNC integral of Eq. 5 and the PY integral of Eq. 7 along with Eq. 6 of Section II, are shown in Fig. 4. It is observed that, in all cases, the perturbation of the membrane structure extends 20 Å from the edge of the cylinder, with strong oscillations in the correlation function separated by ~5 Å. Moreover, the 2D-HNC and 2D-PY theories give similar results, both qualitatively and quantitatively, except for a small deviation around 3 Å. We will therefore only discuss the results in terms of the different radii without referring to the specific approximation used.

For the 2.5-Å cylinder radius, shown in Fig. 4, the average density is higher than its bulk unperturbed value next to the cylinder up to 1.75 Å. The average lipid density is lower for the next region, between 1.75 Å and 21.5 Å, with one small peak between 5 Å and 6 Å and a second peak between 9.5 Å and 10.5 Å. The average lipid density is higher between 13.5 Å and 21 Å and finally relaxes to the bulk value at 21 Å from the edge of the cylinder. The region next to the cylinder can thus be regarded as a crowded layer with a lipid density higher than the uniform bulk value, and this region is followed by a depletion layer with a lipid density lower than the uniform bulk value. This trend is not observed in the case of the 5-Å cylinder radius, where the average lipid density is lower than its bulk value between 0 Å and 10 Å from the edge of the cylinder. The depletion layer is followed by a crowded region, spreading from 10 Å to 21 Å. The crowded layer of the 5-Å cylinder radius is more extended and has a higher density than for the crowded layer of the 2.5-Å cylinder radius. The oscillations in density of these two curves are completely in phase, i.e., when the 2.5-Å cylinder radius curve is at a maximum, the 5-Å cylinder radius curve is also at a maximum. Finally, a trend similar to that of the 5-Å cylinder is observed for the 9-Å cylinder radius. The depletion layer extends from 0 to 5 Å, and the crowded layer extends from 5 to 20 Å. The depletion layer around the 9-Å cylinder has a slightly lower lipid density and is thinner than that for the 5-Å cylinder. In contrast, the crowded layer of the 9-Å cylinder has a higher lipid density and is longer than that for both the 2.5-Å and 5-Å cylinders. In a general way, the lipid density next to the edge of a protein with a radius larger than 2.5 Å, in a DPPC bilayer, is lower than its bulk value, and this depletion layer is followed by a crowded region where the lipid density is higher than its bulk value. The lateral perturbation of the lipid density is within ~25 Å of the edge of the cylinder.

View larger version (17K): [in this window] [in a new window]   FIGURE 4   Radially averaged lipid density (r) around the cylinders as calculated from the HNC integral in Eq. 5 (solid lines) and from the PY integral in Eq. 7 (dashed lines) for 2.5-Å-, 5.0-Å-, and 9.0-Å-radius cylinders. The origin of the axis corresponds to the edge of cylinders, and the mean lipid density of the unperturbed membrane is represented by horizontal lines.

The depletion layer next to the protein inclusion edge increases the area available for lipid molecules in this region. To quantify the increase in area per molecule in this region, the average area per lipid molecule for lipid molecules within the layer of a given size next to the cylinder edge was calculated as follows. First, the number of carbon atoms is retrieved by integration of the density of the carbon atoms over the layer surface. The density curves are shown in Fig. 4 and are obtained with the 2D-HNC theory. Next, the number of lipid molecules is estimated by dividing the number of carbon atoms in the layer by the number of carbon atoms in a single lipid, and again dividing by 2 to get the number of lipids for only one leaflet of the membrane. In fact, this is an approximation because carbon atoms can belong to different lipid molecules. Finally, the surface area of the given layer around the protein inclusion is divided by the number of lipid molecules to give the area by lipid molecule. Results for different layer sizes and different protein inclusion radii are given in Table 1 (recall that the area per lipid molecule for an unperturbed membrane is 62.9 Å²/molecule). In a general way lipid molecules next to the cylinder edge have a greater area per molecule than lipid molecules in the unperturbed membrane. This observation is in accord with results of Husslein et al. (1998), who performed a MD simulation of a hydrated diphytanolphosphatidylcholine lipid bilayer containing an -helical bundle of four transmembrane domains of the influenza virus M2 protein. It was observed that the area per lipid molecule in the vicinity of the protein increases to 85 Å²/molecule, as compared to 74.6 Å²/molecule for an unperturbed membrane lipid bilayer under the same conditions. The presence of a depletion layer of lipids around a protein inclusion and the corresponding increase in cross-sectional area per lipid suggest that there is effectively a long-range repulsion between the lipids and the protein. It seems plausible that the origin of this repulsion is entropic. It has been observed that the lipid chains adopt more ordered configurations (higher carbon-deuterium order parameter) near transmembrane proteins (Woolf and Roux, 1996; Chui et al., 1999). A lipid molecule must reduce its disorder significantly to come close to a protein inclusion, which is entropically unfavorable. The reduced disorder is converted into an effective lipid-protein repulsion. As indicated by Table 1, the effect increases with the size of the protein.

Lipid-mediated forces between protein inclusions

In the context of integral equations, different routes can be taken to calculate the lipid-mediated interaction free energy. If the theories were exact, all of these different methods would yield the same result. However, because the theories are only approximations, there may be some differences. A comparison of the interaction free energy obtained by different routes thus provides a valuable assessment of the internal consistency of the theory. Three different methods were used to compute lipid-mediated interaction free energy (see above): the OZ route based on Eq. 8, the  route based on Eq. 9, and the F route based on Eq. 12. The  route and the F route require explicit consideration of the perturbation of the hydrocarbon density around two protein inclusions. For this reason, they are computationally more expensive than the OZ route, which requires only the correlation functions around a single isolated inclusion. However, the F route offers the possibility of analyzing in detail the origin of the PMF from the nonuniform hydrocarbon density around the protein inclusions.

The results for the lipid-mediated interaction free energy, as a function of the separation distance between two cylinders, are given in Fig. 5. The different approximations for the PMF yield similar results, indicating that the theories are internally consistent. The numerical noise observed in the PMF curves calculated with the  route (Eq. 9) arises from the finite grid combined with the harsh repulsive potential between lipid and cylinders. This is a well-known problem that is also observed in finite-difference Poisson-Boltzmann calculations. For the case of two 2.5-Å cylinders, there is a free energy barrier at a separation distance of 15 Å between the cylinder edges, followed by an attractive free energy well for distances less than 10 Å. The repulsive barrier is ~0.6k_BT and extends from 10 to 20 Å. At separations greater than 20 Å, the protein-protein potential is very small and oscillatory. Finally, at protein-protein contact, the magnitude of the lipid-mediated potential is approximately 2k_BT. For the case of two 5.0-Å cylinders a similar trend is observed, i.e., that there is a free energy barrier at a distance of 10 Å followed by an attractive free energy well for distances less than 5 Å. The repulsive barrier is ~4k_BT and extends from 5 to 20 Å. Again, at separations greater than 20 Å, the protein-protein potential is very small and oscillatory, and at protein-protein contact, the magnitude of the lipid-mediated potential is approximately 8k_BT. This value is in good agreement with, though somewhat more negative than, previous estimates in the literature (Marcelja, 1976; Schroder, 1977; Owicki et al., 1978; Owicki and McConnell, 1979; Pearson et al., 1984; Sintes and Baumgärtner, 1997b). For the case of two 9-Å cylinders there is a free energy barrier at a distance of 5 Å between the cylinder edges, and, in contrast to the case of two 2.5-Å and the two 5-Å cylinders, there is no attractive free energy well. In this case, the repulsive barrier is ~9k_BT and extends from protein-protein contact to 20 Å.

View larger version (23K): [in this window] [in a new window]   FIGURE 5   Lipid-mediated free interaction energy, in kBT units, between two hard repulsive cylinders of 2.5-Å radius (top), 5-Å radius (middle), and 9-Å radius (bottom) as a function of their separation distance between cylinder edges. Results obtained from the 2D-HNC closure: free energy difference (Eq. 9, solid line), integration of force (Eq. 12, dashed line), and the PMF (Eq. 8, long dashed line). Results obtained from the 2D-PY closure: integration of force (Eq. 12, dot-dashed line).

The lipid-mediated forces calculated on the basis of Eq. 12 as a function of the separation distance between two cylinders are shown in Fig. 6. Again, as for the PMF curves, the noise arises from the finite-grid effect in combination with the harsh repulsive potential between lipid and cylinders. In all cases there is a lipid-mediated force present up to as much as 22 Å between cylinder edges. Then, when two 2.5-Å protein inclusions at a great distance apart are brought closer, a small repulsive force arises when the separation distance is on the order of 20 Å. This force oscillates slightly around zero with a separation of 5 Å. This separation between oscillations corresponds approximately to half of a lipid diameter and is exactly the value of the distance between carbon atoms of two different aliphatic chains suggested by the extramolecular correlation function. The force becomes attractive at 5 Å and then repulsive at protein contact. For the case of two 5-Å cylinders, the effective interaction is again repulsive from ~20 Å to 10 Å and is then attractive from 10 Å up to contact between the cylinders, again with oscillations separated by 5 Å. For the case of two 9-Å cylinders the repulsion acts at a larger distance, from ~20 Å to 5 Å, and the range of attraction is from 5 Å up to cylinder contact. As for the 2.5-Å and 5-Å protein curves, the oscillations are again separated by 5 Å. The greatest difference between 2D-HNC and 2D-PY approximations occurs for the 9-Å cylinder curves. However, these two curves are qualitatively similar in shape.

View larger version (20K): [in this window] [in a new window]   FIGURE 6   Lipid-mediated forces between two cylinders as a function of their separation between cylinder edges. See the Fig. 4 legend for details.

As shown by Eqs. 11 and 12, the lipid-mediated forces between two protein inclusions arise from the asymmetry in the lipid density around the cylinders. To understand the origin of these forces as well as the asymmetry in lipid density, a series of contour plots of lipid aliphatic chain density around the protein cylinders was made. For each system of two protein cylinders of the same size three separation distance were chosen, giving a total of nine contour plots. The first separation distance corresponds to the distance where the interaction energy begins to increase in Fig. 5, i.e., a separation distance of 20 Å for every protein size. The second separation distance was chosen to see what happens when the interaction energy is at the maximum in Fig. 5, i.e., a separation of 15 Å for 2.5-Å-radius protein cylinders, 10 Å for 5-Å-radius protein cylinders, and 5 Å for 9-Å-radius protein cylinders. The last separation distance is 2 Å for each protein cylinder size, corresponding to the case where the two protein cylinders are very close to each other. The series of contour plots is given in Fig. 7. In this figure, we see that the depletion layer of the two 9-Å-radius protein cylinders is thinner than for the two 5-Å-radius protein cylinders, but it is the depletion layer of the two 2.5-Å-radius protein cylinders that is the thinnest. Furthermore, for the 5-Å-radius and the 9-Å-radius protein cylinders, the maximum of interaction energy corresponds with the beginning of an overlap of the depletion layers of each protein cylinder.

View larger version (138K): [in this window] [in a new window]   FIGURE 7   Contour plots of lipid aliphatic chain density around two protein inclusions at different distances. The cylinder radius (r, in Å) and the distance between protein edges (d, in Å) are shown at the top of each plot. Dark blue represents protein inclusions, light blue corresponds to a lower density than the bulk, red corresponds to a higher density than the bulk, and cyan corresponds to the bulk density.

The differences between the results for the two protein inclusions of different sizes shown in Fig. 5 for the lipid-mediated free interaction energy and in Fig. 6 for the lipid-mediated forces provide some insight into the origin of protein inclusion size effect in lipid-mediated effective interaction. This effect arises from a difference of lipid packing around the cylinders of different sizes. It is important to note that, for a given protein inclusion size, the point where the attraction between two inclusions begins in Fig. 6 corresponds to a free energy maximum in Fig. 5. As stated above, for the two 2.5-Å inclusions, the lipid-mediated free energy maximum is at 15 Å, at 10 Å for the two 5-Å inclusions, and at 5 Å for the two 9-Å inclusions. These distances correspond to the diameter of three aliphatic chains for the case of two 2.5-Å inclusions, two aliphatic chains for the case of two 5-Å inclusions, and only one aliphatic chain for the case of two 9-Å inclusions. Furthermore, for the two 5-Å inclusions and the two 9-Å inclusions, these distances correspond exactly to the complete overlap of the depleted regions of each inclusion (see Figs. 4 and 7, middle column). This suggests that the lipid molecule between the two 5-Å inclusions is bound to both proteins, and that this lipid would have to be entirely removed to bring these two cylinders closer. In contrast, for the two 9-Å inclusions, there is only a single aliphatic chain between these two inclusions before the attraction occurs. The curvature of the 5-Å protein inclusion is perhaps too pronounced for a single lipid to be able to place both its aliphatic chains along the edge of one cylinder of this size. However, the curvature of the 9-Å cylinder does not seem to be too pronounced, and a single lipid can then place its two aliphatic chains along the edge of one cylinder of this size. The suggested configurations of lipid packing around a 5-Å cylinder and a 9-Å cylinder are shown schematically in Fig. 8. For the two 2.5-Å inclusions, the free energy maximum being at 15 Å, the interpretation of the result is not obvious.

View larger version (14K): [in this window] [in a new window]   FIGURE 8   Two-dimensional top view of the lipid packing around a 5-Å-radius cylinder (A) and a 9-Å-radius cylinder (B). Proteins are represented by circles and lipid molecules are represented by ellipses. See text for discussion.

The lipid-mediated forces obtained with the two 5-Å cylinders are comparable to those obtained by Sintes and Baumgärtner (1997b) with different protein radii. They used Monte Carlo computer simulations based on a simple model of the lipid bilayer, where the bilayer is represented with 2 × 500 lipid molecules and where each lipid molecule was modeled by a flexible chain composed of five monomers. The proteins were modeled by two hard transbilayer cylinders. They observed that the lipid-mediated forces were almost independent of the protein radius. In contrast, we observe a size effect in the effective interactions between two cylinders. However, the presence of both an attractive and a repulsive part in the lipid-mediated protein-protein effective interaction is in qualitative agreement with their results.

	CONCLUSION

TOP ABSTRACT INTRODUCTION THEORETICAL DEVELOPMENTS RESULTS AND DISCUSSION CONCLUSION REFERENCES

The dependence of lipid-mediated interactions between protein inclusions on protein size was investigated using a theory for examining the structure of the hydrocarbon chains around protein inclusions embedded in a lipid bilayer. This theory, based on the HNC integral equation theory for liquids, was recently developed (Lagüe et al., 1998) and uses the exact lateral density-density response function of the hydrocarbon core as an input to the calculations. This response function was computed from configurations taken from the MD simulation of a pure DPPC lipid bilayer of Feller et al. (1997). In addition to the original theory, where the lipid density around protein inclusions and lipid-mediated forces acting between two protein inclusions were computed, a new physical quantity was calculated: the lipid-mediated free energy acting between two protein inclusions. The theory was also extended to the PY integral equation (Hansen and McDonald, 1986) for comparison with the results obtained using the HNC closure.

We first examined the structure of the DPPC bilayer in the neighborhood of a single protein modeled as a hard cylinder for three different sizes: a small cylinder of 2.5-Å radius, corresponding an aliphatic chain; a medium cylinder of 5-Å radius, corresponding an -helical polyalanine protein; and a large cylinder of 9-Å radius, representing a small protein such as the gramicidin channel (Woolf and Roux, 1996). The results showed that the average lipid order is perturbed over a distance of 20 Å from the edge of the protein. For the 5-Å an