This Article Free via Open Access: OA OA Abstract Full Text (PDF) Supplement OA All Versions of this Article: biophysj.107.120667v1 94/7/2580    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Google Scholar Articles by Khelashvili, G. Articles by Harries, D. PubMed PubMed Citation Articles by Khelashvili, G. Articles by Harries, D.

Protein Diffusion on Charged Membranes: A Dynamic Mean-Field Model Describes Time Evolution and Lipid Reorganization

George Khelashvili ^*, Harel Weinstein ^* and Daniel Harries 

^* Department of Physiology and Biophysics, Weill Medical College of Cornell University, New York, New York; and Department of Physical Chemistry and the Fritz Haber Research Center, The Hebrew University, Jerusalem, Israel

Correspondence: Address reprint requests to George Khelashvili, Tel.: 212-746-6539; E-mail: gek2009{at}med.cornell.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MODEL RESULTS AND DISCUSSION CONCLUDING REMARKS APPENDIX SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
As charged macromolecules adsorb and diffuse on cell membranes in a large variety of cell signaling processes, they can attract or repel oppositely charged lipids. This results in lateral membrane rearrangement and affects the dynamics of protein function. To address such processes quantitatively we introduce a dynamic mean-field scheme that allows self-consistent calculations of the equilibrium state of membrane-protein complexes after such lateral reorganization of the membrane components, and serves to probe kinetic details of the process. Applicable to membranes with heterogeneous compositions containing several types of lipids, this comprehensive method accounts for mobile salt ions and charged macromolecules in three dimensions, as well as for lateral demixing of charged and net-neutral lipids in the membrane plane. In our model, the mobility of membrane components is governed by the diffusion-like Cahn-Hilliard equation, while the local electrochemical potential is based on nonlinear Poisson-Boltzmann theory. We illustrate the method by applying it to the adsorption of the anionic polypeptide poly-Lysine on negatively charged lipid membranes composed of binary mixtures of neutral and monovalent lipids, or onto ternary mixtures of neutral, monovalent, and multivalent lipids. Consistent with previous calculations and experiments, our results show that at steady-state multivalent lipids (such as PIP₂), but not monovalent lipid (such as phosphatidylserine), will segregate near the adsorbing macromolecules. To address the corresponding diffusion of the adsorbing protein in the membrane plane, we couple lipid mobility with the propagation of the adsorbing protein through a dynamic Monte Carlo scheme. We find that due to their higher mobility dictated by the electrochemical potential, multivalent lipids such as PIP₂ more quickly segregate near oppositely charged proteins than do monovalent lipids, even though their diffusion constants may be similar. The segregation, in turn, slows protein diffusion, as lipids introduce an effective drag on the motion of the adsorbate. In contrast, monovalent lipids such as phosphatidylserine only weakly segregate, and the diffusions of protein and lipid remain largely uncorrelated.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MODEL RESULTS AND DISCUSSION CONCLUDING REMARKS APPENDIX SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Many proteins that peripherally adsorb on lipid membranes contain structured domains that target the protein to the bilayer; examples include the C2 (1), PH (2), FERM (3), and BAR (4) domains. Many of these domains act by specifically binding to a particular lipid species, like the PH domain that binds phosphatidylinositol 4,5-bisphosphate (or PIP₂) lipids. However, an apparently different type of targeting is achieved by numerous other proteins that contain natively unstructured clusters of basic residues, such as the well-studied examples of the GAP43, GTPase K-Ras, and MARCKS (5–9). The use of positively charged residues for targeting may come as no surprise, as cellular plasma membranes typically contain 20% anionic lipids. This affords a simple mechanism for protein-lipid binding that is essentially nonspecific, yet able to confine proteins to membrane interfaces.

This simple molecular picture has been challenged by recent theoretical and experimental evidence suggesting that the major anionic lipid component in many cells, phosphatidylserine (or phosphatidylglycerol), might not be the major participant in peripheral protein binding (10–13). Instead, the typically multivalent phosphoinositides, such as PIP₂ or even phosphatidylinositol 3,4,5-bisphosphate, are more likely implicated in segregation close to peripherally adsorbed proteins. This is interesting, because phosphoinositides are known to play an important regulatory role at the plasma membrane. Despite the fact that they constitute typically only 1% of membrane composition, these minority lipids can act at sites of regulation at least partly by electrostatic association with peripheral and embedded proteins (14). Concentrating PIP₂ at the site of protein adsorption is therefore a likely mechanism for local and specific recruitment. It has been suggested that segregated lipids can subsequently be released upon cellular changes, e.g., in Ca⁺² concentrations. This provides a way to control the amount of free PIP₂ in the membrane, and hence a mechanism for regulating PIP₂ known to participate in cellular signaling processes such as enzyme activation, endocytosis, and ion-channel activation (15).

To begin to understand why electrostatic targeting could primarily be achieved by polyvalent rather than the more abundant monovalent lipids, we must focus on the forces that underlie this protein-lipid interaction. Theory and experiments show that the attraction of positively charged protein domains to the oppositely charged membrane is due not only to Coulombic interaction, but is also entropically driven. This entropic gain is due to the release of counterions that were previously confined locally to the vicinity of the isolated protein or membrane by the requirement to preserve charge neutrality (16–19). Upon protein-membrane binding, these counterions are no longer required and are released to produce a translational entropy gain in the bulk solution, while the protein and membrane neutralize one another. At the same time, to allow maximal counterion release, charged lipids can migrate in the membrane plane toward the protein adsorption site to fully compensate charges on the protein, causing demixing (20). But this local lipid demixing comes at an entropic cost. The lower the membrane charge density, the higher the cost for the necessary lipid segregation and demixing upon protein binding. At equilibrium, the system has reached some compromise between maximal counterion release and minimal demixing.

Experiments have suggested that that PIP₂ preferentially segregates at sites of charged protein adsorption (10). This is reasonable because multivalent lipids should incur a smaller demixing penalty and larger counterion release entropy per segregated lipid, simply because each of them carries a larger charge. Recent theoretical studies predict that multivalent lipids should indeed segregate more than monovalent ones, and that the binding free energy to rigid macromolecules as well as to polyelectrolytes is significantly stronger for such lipids (11–13). But recognizing the dynamic nature of the adsorption problem raises the possibility that the kinetic energy of each adsorbing protein allows it to move so quickly on lipid membranes that some lipids rarely manage to segregate at all. Conversely, lipids may rearrange so quickly around an adsorbing protein that the protein appears stationary to them, creating a transient binding site and thus impeding the protein's motion in the membrane plane.

If we consider free lipid diffusion in the membrane plane under the influence of the field exerted by the oppositely charged protein, we might conclude that multivalent lipids are more mobile (move faster) than monovalent ones, even if their diffusion in the absence of the protein field is similar. This conclusion simply follows if we assume lipid diffusion to be directed by spatial (second-order) derivatives of the electrochemical potential that in turn are linearly proportional to the product of charge and local potential. Higher mobility near the charged macromolecule may allow multivalent lipid segregation that can happen on timescales fast enough to follow the protein random motion as it moves on the membrane surface. This not only reflects a stronger adsorption free energy for protein on PIP₂-containing membranes, but also highlights the important role of the minority multivalent lipids in lipid segregation. Note, however, that living cells often exhibit diffusion behavior of membrane components (lipids, proteins) that is much different than in vitro (21).

To quantify the combined kinetic effect of many lipid species interacting with peripheral proteins, it is essential to be able to calculate the steady state of adsorbing macromolecules in a way that will include all important degrees of freedom in a self-consistent way. Previous theoretical studies have shown that equilibrium distributions can be predicted from a self-consistent mean-field model based on the modified Poisson-Boltzmann (PB) equation (20,22). But to address the concerted action of protein adsorption and lipid segregation, we extend this strategy here by using a dynamic propagation method to efficiently derive steady-state configurations for membranes interacting with macromolecules. Our model is an application of a time-dependent self-consistent mean-field approach that has been used to address similar problems (for examples, see (23–26)). The numerical scheme we developed provides not only the adsorption free energy and charge density distribution on the membrane at steady state, but also can be used to gain additional dynamic information on the segregation process.

Our method uses an atomic level representation in three dimensions, and takes into account lateral reorganization and demixing of lipids during adsorption. Lipids are allowed to move in the membrane plane according to a diffusion-like Cahn-Hilliard (CH) equation (27), where segregation rates are in proportion to the Laplacian of their chemical potential. The local chemical potentials are derived from the free-energy functional that depends on local lipid component densities and are calculated using results from nonlinear PB theory.

The membranes we consider here are binary or ternary mixtures of neutral (zwitterionic) lipids, as well as negatively charged ones bearing one or more charges. At steady state, solutions of our two-dimensional CH equations in the membrane plane match the equilibrium solutions of the corresponding modified PB equation. We find, in agreement with previous studies (11–13), that static positively charged polypeptides and proteins positioned near net-acidic membranes of biologically relevant composition preferentially recruit to the adsorption zone multivalent lipids such as PIP₂, rather than monovalent ones.

We then extend our model by using a dynamic Monte Carlo scheme to consider lipid dynamics combined with adsorbate diffusion on the membranes. This dynamic model allows us to conclude that it is the composition of the membrane on which proteins are diffusing that determines whether lipids will be sequestered. In particular, we find that PIP₂ lipids can be expected to segregate near oppositely charged proteins and thereby slow down the diffusion of the protein. An important prediction of our model is that PIP₂ lipids will be able to diffuse in concert with the retarded adsorbed proteins, while lipids such as phosphatidylserine (PS) will only weakly segregate and, in this case, both protein and lipid diffusion will be largely uncorrelated.

The implied consequences of the detailed dynamic picture are clear: by virtue of their charge alone, PIP₂ lipids can be sequestered and retained for extended periods of time in the vicinity of oppositely charged peripheral proteins. This may serve as a possible mechanism for the formation of mobile lipid microdomains that diffuse slowly in the membrane plane, allowing both the lipid and the protein the time needed to act in concert.

	MODEL

TOP ABSTRACT INTRODUCTION MODEL RESULTS AND DISCUSSION CONCLUDING REMARKS APPENDIX SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Free-energy functional
Consider charged proteins and lipid bilayers immersed in an aqueous solution. The solution also contains a symmetric 1:1 electrolyte of concentration n₀, corresponding to a Debye length

(1)

Here, k_B is the Boltzmann constant, T is the temperature, e is the elementary charge, ₀ is the permeability of free space, and _w = 80 is the dielectric constant of the aqueous solution. At the physiological conditions modeled in our calculations, at T = 300 K and n₀ = 0.1 M electrolyte (monovalent salt) concentration, _D 10 Å.

We consider the limit of low surface density of adsorbing proteins, so that interactions between proteins are negligible. We represent the adsorbing protein in full-atomistic three-dimensional details, whereas the membrane is considered as a two-dimensional fluid, allowing us to treat lipid headgroup charges in the continuum representation as usual in regular solution theory (28). Hybridizing these atomistic and continuum representations provides a realistic mesoscale description of the electrostatic problem, while still allowing the entire system to be described at the same mean-field level.

Fig. 1 shows a unit simulation cell containing a single charged macromolecule (poly-Lysine 13 residues long) adsorbed on a membrane. Both the macromolecule and the lipid membrane are treated as a low dielectric material, with dielectric constant _m = 2 within adsorbate and membrane. Each atom on the adsorbing protein is assigned a radius and a partial charge. The oppositely charged membrane consists of mixtures of net-neutral, mono-, and multivalent lipids. Each membrane bilayer is composed of two adjoined lipid monolayers, forming a slab of thickness d. We assume that the membrane is stiff with respect to any weak-deforming forces exerted by the protein, and we therefore keep the membrane flat in all our calculations.

View larger version (9K): [in this window] [in a new window]   FIGURE 1  Schematic view of a simulated unit cell containing a protein (or peptide) adsorbed on a membrane. For illustration, we use the basic poly-peptide Lysine, 13 residues long (Lys13). The membrane is represented by a rectangular-shaped slab of thickness d. Charges on the lipid headgroups are represented by a continuous surface charge density. The distance of nearest approach between protein and membrane surfaces is h. In all our calculations the dielectric constant for the membrane interior as well as protein is m = 2, the dielectric constant of the aqueous environment is m = 80, and the Debye length in the electrolyte solution is D = 10 Å.

(2)

For simplicity, we assume here that all lipids have the same lateral area per headgroup, a (though the model can easily be extended to include different headgroup areas as in Andelman et al. (29)). Denoting the valency of the i^th lipid species by z_i, we define the local surface charge density as

(3)

The adsorbing macromolecule is first considered as fixed in space at a distance h from the membrane surface, where h is defined as the shortest distance between the van der Waals surfaces of atoms on the adsorbate and lipid headgroups (Fig. 1).

The free-energy functional for the system can be written as a sum of the electrostatic energy, salt ion translational entropy, and lipid mixing entropy in the membrane plane (20):

(4)

The system's electrostatic energy is given, as usual, by

(5)

Here = e/k_BT is the dimensionless (reduced) electrostatic potential, being the electrostatic potential. The integration of this term must be carried out over the entire space. The contribution from the translational entropy of mobile (salt) ions in solution is

(6)

where n₊ and n_– represent local concentrations of positive and negative mobile electrolyte ions, respectively, and n₀ is the electrolyte concentration in the bulk solution.

The two-dimensional mixing entropy of mobile lipid molecules contribution is given by

(7)

The sum extends over m lipid species, represents the average composition of the i^th lipid species, and the integral is taken over the membrane surface.

Functional minimization of F with respect to the mobile ion concentrations leads to the nonlinear PB equation (22,30–35):

(8)

Solving this equation yields the electrostatic potential . The additional minimization of F with respect to the lipid compositional degrees of freedom leads to a second differential equation on the bilayer boundary that should be solved simultaneously with Eq. 8 (20,22). Here, however, we do not solve this boundary equation. Instead, we use a dynamic propagation scheme that reaches the equilibrium lipid distribution in the long-time limit, as shown in the next section.

Lipid propagation in time
We now derive the dynamic equations that govern the time-dependent lipid rearrangement in our model. The steady-state solutions to these equations correspond to the lipid distribution in the membrane plane that minimizes the free-energy functional F with respect to all _i-values. The starting point is the continuity equations for all lipid species. The equation for the i^th lipid species has the form (36)

(9)

where is the local (lipid) current at position corresponding to a mole fraction of the i^th lipid species at time t, and represents the flow velocity of the i^th lipid species at time t and located at on the membrane surface.

We next relate the local current of the i^th lipid species and the local gradients of its electrochemical potential. The idea is to assume that gradients in the chemical potential determine the velocities of lipids migrating in the membrane plane, in the spirit of Diamant and Andelman (26), and references therein. We define the electrochemical potential within our model:

(10)

Here N_i is the number of lipids of species i, and µ°_i represents the standard chemical potential for the i^th lipid species that is independent of _i. Using Eqs. 2 and 4 in Eq. 10, we find (12)

(11)

We assume that the above expression is valid not only in the equilibrium state of the system, but can be used as well to obtain local instantaneous chemical potentials for lipids as the system evolves in time toward the steady state. Then, the generalized force exerted locally on the i^th lipid species is given by –µ_i (37). This force is balanced by the friction experienced by all lipid species. Hence, we can write the following set of balance equations (37),

(12)

where D_ij values are the so-called diffusivities, or the elements of the diffusion matrix. Based on the experimental findings by Golebiewska et al. (10), showing that the effective diffusion coefficients for uncharged, mono-, and multivalent acidic lipids are similar, we assume here that all diffusivities are equal and simplify Eq. 12 by setting D_ij D_lip for all lipid types. With this approximation, the only distinguishing characteristic between the net-neutral, mono-, and multivalent lipids in our model is their net headgroup charge.

Combining Eq. 12 with the definition of currents, the lipid incompressibility constraint, and with the flux neutrality condition leads to the relationship between currents and gradients in electrochemical potentials:

(13)

With the above expression, Eq. 9 simplifies to

(14)

which is the desired CH equation for the time-evolution of lipid compositions.

Equation 14 is solved for all lipid species i 1 self-consistently with the nonlinear PB Eq. 8, so that at each iteration the local surface charge density and the electrostatic potential gradient on a membrane surface are linked through the boundary condition,

(15)

with z₀ describing the position of the charged interface (Fig. 1). The diffusion of i = 1 follows from the material conservation constraint Eq. 2. We note that within the CH formalism, lipid fractions (i = 1, ..., m) are globally conserved fields and, for a given time instance t, they obey

(16)

where S_A is the area of a membrane surface A. We stress that, in our model, the lateral motion of lipids is entirely due to the presence of the adsorbed protein. We do not consider any fluctuations, i.e., thermal motions of lipids, which in principle can be introduced into the CH dynamics using a stochastic noise term (27). Describing lateral reorganization of lipid molecules by the CH equations implies that the motion of lipid molecules in the membrane plane obeys the laws of normal (regular) diffusion. In contrast to artificial membranes, lipids in cellular membranes may also undergo anomalous diffusion (see for example (38,21) and reference therein). Therefore, our model does not attempt to describe a priori any anomalous or hop-diffusion dynamics, neither for lipids nor for the adsorbed macromolecule. However, within our model such diffusion could still be the result of lipid-protein interactions.

Mobile protein: hybrid Cahn-Hilliard and dynamic Monte Carlo
To this point, we have assumed that the adsorbing macromolecule remains stationary at a distance h from the membrane surface and only let lipids rearrange laterally in the two-dimensional membrane plane. We can further extend our model to account for the diffusional motion of the protein in the membrane plane within a dynamic Monte Carlo (DMC) scheme (38–44). Within DMC, we start with some particular system configuration as the initial state, and then generate a sequence of other possible trial configurations. Acceptance of these moves represents a stochastic dynamic trajectory of the system.

Our goal is to couple the DMC scheme for the motion of the adsorbing macromolecule on the membrane surface with the CH formalism for lateral rearrangement of lipids (which we term the CHDMC method). Protein and lipids simultaneously diffuse, each with their own typical diffusion rate. The two corresponding timescales in our model are d²/D_lip, for the lateral diffusion motion of lipids, and d²/D_prot, for the diffusion of adsorbing macromolecule, where d is the lattice constant, D_lip is the lipid diffusion coefficient, and D_prot is the diffusion coefficient of the macromolecule on a homogeneous membrane. The ratio of the two diffusion coefficients D' = D_prot/D_lip determines how close the two relevant timescales are, and therefore also reflects the coupling strength of these different modes of motion. Simply stated—if protein and lipid diffuse at similar rates, their motion is expected to be highly correlated. However, if the protein diffuses much faster than lipid, lipid rearrangement will not achieve complete relaxation, and in the limit of D_prot >> D_lip, the motions of lipid and protein can be expected to be largely uncorrelated.

We start the CHDMC simulation with a homogenous distribution of lipids on the membrane, and with the adsorbing macromolecule at distance h from the membrane surface. At the initial step, we allow the adsorbate center of mass to make a random displacement in the two-dimensional membrane plane. In accordance with the fluctuation-dissipation theorem, we treat the random move as a combination of two independent displacements, each of size (27), in two membrane plane directions. Here _G is a Gaussian random number with zero mean and unit variance, and the dimensionless time step t' is related by the CH Eq. 14 to the real-time step t through

(17)

The move is accepted with a Metropolis-like criterion with the usual transition probability of W = 1 if F_new F_old, and if F_new > F_old. Here F_old and F_new are the adsorption free energies of the "old" and trial states of the protein-membrane system, respectively. If a trial move is accepted, the macromolecule is advanced to the "new" position, and the CH equations for lipids are solved taking into consideration the new position of the adsorbate. If, on the other hand, the trial move is rejected, the protein remains at its previous position and lipids rearrange with respect to the old location of the adsorbate. Because time is set by the CH equations, the dimensionless time is updated by t', regardless of the outcome of the DMC move, and the next stochastic step for the protein is attempted. Note that the adsorbate DMC move sizes are not arbitrary, but rather are determined by t'. Therefore, in the context of the dynamic scheme, each rejected DMC move can be viewed as a time-interval within which, on average, the protein's position does not change appreciably due to favorable local interactions with underlying charged lipids. Clearly, such an assumption requires the discretization of the CH equations with sufficiently small t' as compared to all relevant timescales in the system. With our choice of t', we found a rejection rate of only 1.5–4.6% (depending on the value of D').

A simple limiting case for the CHDMC algorithm is when lipid distribution on the membrane remains homogenous during the entire simulation. For this case, proteins should perform free diffusion with their effective diffusion coefficient D_prot unaffected by the presence of the membrane. Any protein slowing seen with our model is because of energetic barriers to lateral motion in the membrane plane that arise when charged lipids segregate around the adsorbate. In principle, the adsorbate might be able to overcome these barriers by diffusing away from the membrane and subsequently readsorbing. However, we also show in Results and Discussion that diffusion perpendicular to the membrane plane has an energetic cost. For example, when a globular protein carrying a surface charge of 1e^– per 93 Ų drifts 1 Å away from a membrane with neutral to monovalent lipid ratio of 80:20, we find a loss of k_BT in adsorption free energy (see also Fig. 2 a). In contrast, when moving the same distance as a result of a single lateral DMC step in the same system, the protein will face a significantly lower free-energy barrier of, at most, 0.1 k_BT. Motion in the perpendicular direction (increasing h) is, therefore, far less likely than lateral motion, where barriers are smaller. This finding allows us to ignore variations in h during the DMC run, as these are associated with prohibitively large energies. Extension to three-dimensional diffusion can be easily incorporated if required in future applications.

View larger version (9K): [in this window] [in a new window]   FIGURE 2  Adsorption of a spherical macroion (Rp = 10 Å) onto binary (PC/PS) lipid membranes with (a) Steady-state (equilibrium) adsorption free energy (in kBT units) as a function of macroion-membrane separation calculated using our method (solid line) and by the method used in May et al. (20) (dashed line). (b) Time sequence of radial profiles of the local fraction of PS lipids, for macroion-membrane separation of h = 3 Å. The steady-state distribution of PS lipids from the calculations by May et al. is shown for comparison (dotted-dashed line).

Simulation details
We focus on two types of mixed lipid bilayers. The first is composed of binary mixtures (m = 2) of neutral (phosphatidylcholine, i.e., PC) and acidic monovalent (phosphatidylserine, i.e., PS) lipids. The second is a ternary mixture (m = 3) of uncharged (PC), monovalent (PS), and acidic multivalent (PIP₂, valency of –4 at neutral pH) components. The discretized version of the CH equation (Eq. 14) for the case of ternary mixtures is presented in the Appendix. For all calculations the lipid membrane was modeled as a low dielectric slab (_m = 2) of dimensions 256 Å x 256 Å x 10 Å. Assuming that a = 65 Ų is the area per lipid headgroup, the slab dimensions correspond to roughly 1008 lipid molecules in one membrane layer. We note that, in principle, the electrostatic properties on one side of the bilayer may have an impact on those on the other side. Under physiological conditions, when this coupling has been shown to be weak (30,45). We have verified the above inequality by performing calculations on membranes of different thickness, d 10 Å. We found that the value of the electrostatic potential on the slab surface did not change when we varied the slab thickness. Therefore, to simplify calculations we completely decouple electrostatically the two membrane interfaces, and treat the lipid slab as a leaflet of thickness d = 10 Å (Fig. 1).

We note that attempting to describe dielectric properties of the interior of a lipid membrane by a uniform dielectric constant may not always be appropriate. Since solvent molecules generally penetrate deep into the lipid headgroup region, the dielectric constant in this area can reach much higher values than assumed in our work. To test the effect of the _m value on the predictions from our model, we compared the results for a membrane slab of thickness d = 10 Å with uniform _m = 2 to those obtained for a membrane described as two fused slabs, one of d = 5 Å and _m = 20 (interfacial headgroup region), and the other d = 10 and _m = 2 (hydrocarbon tail region). The compared adsorption free energies were within 1.5% of each other. We concluded that treating the entire lipid membrane as a low dielectric media was indeed an adequate representation of the systems investigated here.

At each iteration of the CH Eq. 14, the nonlinear PB Eq. 8 was solved using a modified version of the publicly available open-source software APBS, Ver. 0.4.0 (46). The system was placed on a 256 Å x 256 Å x 256 Å cubic grid with grid spacing of 1 Å, and the nonlinear PB was discretized with the finite-difference method. The APBS software was modified to include periodic boundary conditions in the (xy) bilayer in-plane directions. The charge, ion accessibility, and dielectric maps were configured and supplied to APBS. After each dynamic step, these maps were updated and fed back to the PB solver.

The Cahn-Hilliard equations were discretized on a 256 Å x 256 Å square lattice with 1 Å grid spacing. Convergence of the CH equations to equilibrium was checked by confirming that within numerical uncertainty the lipid electro-chemical potentials for all lipid species at steady state are uniform across the membrane surface. In all simulations, we chose a real time step of 200 ps.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION MODEL RESULTS AND DISCUSSION CONCLUDING REMARKS APPENDIX SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
As a first test of the method, we considered the adsorption of a uniformly charged spherical macroion onto an oppositely charged lipid membrane, consisting of binary mixtures of monovalent (PS) and neutral lipids (PC). The equilibrium adsorption free energy of charged macroions on binary membranes, and the steady-state distribution of lipid molecules around the macroion, have previously been derived using the theoretical model of May et al. (20). To validate our numerical methodology, we repeated the calculation performed by May et al. (20), with the expectation that our dynamic minimization scheme of the free energy described in Model should yield a steady-state solution identical, within numerical uncertainty (see Simulation Details), to theirs. The free energies were obtained in May et al. by using the alternative method of solving Eq. 8 with a boundary condition that accounts for lipid mobility. We then considered ternary lipid mixtures of monovalent, multivalent, and neutral lipids, and again applied our model to a uniformly charged spherical macroion (i.e., protein) adsorbing onto such membranes.

Next, we present the model results for the basic poly-peptide (Lysine-13, i.e., Lys¹³) adsorbing onto similar net acidic membranes consisting of binary and ternary mixtures of lipids.

Finally, we present results from a set of simulations in which the adsorbing macromolecule is allowed to diffuse. We compare diffusion rates of the membrane-bound and free macromolecules, and show that lipid composition plays a crucial role in regulating diffusion properties of peripheral, bound proteins.

Macroion adsorption on mixed neutral and monovalent lipid membranes
Fig. 2 a shows binding (adsorption) free energies (in k_BT units) as a function of membrane-macroion separation from our model and as derived from the free-energy functional minimization scheme implemented by May et al. (20). We consider a lipid membrane composition of Assuming an area per lipid headgroup of 65 Ų for both PS and PC lipids, this composition corresponds to one negative charge per 325 Ų of membrane area. The adsorbing spherical macroion of R_p = 10 Å radius mimics a globular protein carrying a surface charge of 1e^– per 93 Ų—that is, 3.5 higher than the average membrane charge density. The macroion was placed at successively higher distances from the membrane surface ranging from 3 Å to 15 Å, and for each separation, lipids were allowed to evolve toward the steady state with the dynamic method described in the previous section.

Fig. 2 a demonstrates full agreement in binding free energies resulting from the two calculations for all membrane-macroion separations. Fig. 2 b shows the fraction of acidic lipids as a function of the radial distance r from the projected center of the ion onto the membrane surface. We show radial profiles for membrane-ion separation of h = 3 Å, at different times starting from a protein that is adsorbed on a homogenously charged membrane. For comparison, we also plot the theoretically predicted equilibrium charge distribution calculated in May et al. (20).

By following the dynamic evolution of lipid rearrangement around the macroion, we can identify two timescales corresponding to two processes. The first corresponds to a major recruitment of charged lipids toward the interaction zone. This process, driven by the strong electrostatic interactions between oppositely charged macroion and charged lipids, and subsequent counterion release, occurs on short timescales: during the first 100 ns, the fraction of PS lipids close to the ion increases 2.5-fold, whereas during the remaining part of the simulation only minor changes occur close to the macroion. As a consequence of this initial charged lipid sequestration, a deficiency of PS lipids is created starting from r = 20 Å away from the ion resulting in the formation of a depletion-well in the lipid distribution, where _PS(r) < At longer times, the second process of filling up the well begins, as PS lipids from the bulk start flowing into the area of lower electrochemical potential and replace neutral lipids in the process. Fig. 2 b shows the relatively long timescales for the second process, as the depleted zone still persists after 1 µs of simulation.

Comparison of the lipid distribution after 1 µs of simulation with that of the equilibrium state calculated in May et al. (20) also shows good agreement. The small discrepancy in the two plots at distances r = 20–40 Å is due to the slow filling-up process, and further evolution of the system yields even better agreement with the equilibrium plot. The corresponding free energy at this point is already within numerical accuracy. Close to the ion, our calculation predicts a 2.6% higher fraction of PS lipids, over the prediction by May et al. This minor difference can be attributed to the specific details of the two models and the membrane representation in each. In May et al. (20), the lipid membrane is a continuous planar surface of charge density _memb. This surface divides space in two with dielectric constant _m = 0 in the nonaqueous volume. The macroion dielectric constant was also taken there to be _m = 0. In contrast, in our model, the lipid bilayer is represented by a rectangular slab of finite thickness d and the desired local charge density is achieved by placing (partial) point charges 1 Å apart on the membrane surface. The dielectric constant inside the slab, as well as inside the adsorbing macroion, is taken as _m = 2 in our model—a value probably more representative of the dielectric properties of macromolecules such as proteins and lipids. This higher _m in the membrane and macroion may allow slightly more charged lipids to segregate (less repulsion from image charges).

To conclude, the set of simulations presented above demonstrates that the numerical solution of Eq. 14 converges at long times to the equilibrium state as it should. Thus, our numerical scheme for calculating free energies is fully validated and allows us to follow lipid diffusion in the membrane plane.

Multivalent versus monovalent lipid segregation
Fig. 3 details lipid segregation in the adsorption process of a spherical macroion on ternary mixtures of neutral (PC), monovalent (PS), and polyvalent (PIP₂) lipid membranes containing and The spherical macroion again has a radius of R_p = 10 Å and a surface charge of 1e^–/93 Ų.

View larger version (18K): [in this window] [in a new window]   FIGURE 3  Adsorption of spherical macroion (Rp = 10 Å) onto ternary (PC/PS/PIP2) lipid membrane with composition. (a) Normalized fraction * of PIP2 lipids (upper panel) and PS lipids (lower panel) as a function of the radial distance from the macroion, r, at different times from the initial macroion binding. (b) Time sequence for the electrochemical potentials (µ – µ°), reported in kBT units, of PIP2 (upper panel) and PS (lower panel) lipids upon macroion binding.

(18)

Here a is the area per lipid headgroup, the integral is taken over the membrane area around the macroion, and is the lipid composition at the cell boundary (Fig. 3 a). With the above definition, _i measures the excess number of the i^th lipid species around the adsorbate, relative to the respective number in the homogenous membrane. From Eq. 18 we find that after 400 ns the concentration of PIP₂ in the first coordination shell is increased by 150%, whereas PS lipid content increased by only 17%. Furthermore, our calculations reveal that during the adsorption process, the macroion attracts –2.1e additional charges (compared to a homogenous membrane) from the acidic lipids, of which PIP₂ contributes –1.2e, and PS –0.9e charges. In light of the major role observed for PIP₂, despite its very low content in the membrane, we conclude that the binding free energy of the complex is minimized most efficiently by the macroion interacting primarily with PIP₂ lipid molecules.

The preferential increase in multivalent phospholipid concentration at the adsorption region over monovalent ones has been reported in several experiments (8,10,49) as well as in theoretical studies (11–13). In our model, the large difference in sequestration of PIP₂ versus PS lipids is mainly due to the lower entropic cost associated with recruiting a multivalent lipid. The multivalent lipid carries several charges (better electrostatic interaction) but still loses only the entropy of one free lipid when sequestered compared with four PS lipids that would have to be sequestered instead. In our model, this preference also has an additional dynamic aspect. Highly charged lipids tend to be more mobile and move faster under the influence of an external (protein) electric field. This is because the chemical potential is directly related to the product of the electrostatic potential, and lipids valency (Eq. 11), as are the derivatives of the electrochemical potential that determine lipid mobility. Therefore, PIP₂ lipids can be expected to move faster toward the interaction zone, and once sequestered, are expected to remain bound or localized to the protein with preference over PS lipids.

Due to the lower PIP₂ content, these lipids might not always segregate quickly enough to the protein's vicinity despite their high mobility, thus being unable to replace the majority of PS lipids that, initially, are already quite abundant. This is expressed in the later time spans of the simulation as a displacement of the initially segregated PS lipids by PIP₂ lipids that are thermodynamically favored.

Fig. 3 b follows the changes in local chemical potentials with time. Clearly, the tendency at long times is toward achieving uniform chemical potential throughout the membrane plane for each of the three lipid components (µ_PC follows from µ_PS, and Eq. 2). The figure also shows clearly that the gradients in initial chemical potentials are much larger for PIP₂ than PS, although its membrane content is much smaller. For both lipids, uniform chemical potential can only be achieved by major lipid recruitment to the interaction zone.

Poly-peptide adsorption on binary and ternary lipid membranes
A more realistic peptide-membrane interaction is modeled by the adsorption of Lys¹³ on an oppositely charged (mixed) membrane. For this calculation, Lys¹³ was represented in all-atom detail, without blocked ends. Atomic radii and partial charges for each peptide atom (including the C-terminus) were derived from the CHARMM force field (50). Following the findings of Ben-Tal et al. (51), we placed Lys¹³ in a flat conformation next to the membrane such that its plane was parallel to the membrane surface and the minimal distance between the peptide interface and the membrane was h = 3 Å, as shown in Fig. 1. Based purely on the electrostatic calculations, Ben-Tal et al. showed that the free energy of Lys⁵ binding to 2:1 PC/PS membrane was lowest for this particular configuration (51).

We consider two different lipid compositions with the same surface charge density: a binary mixture with and a ternary mixture with and Fig. 4 shows the charged lipid organization for the ternary mixture (Fig. 4, a and b) and binary mixture (Fig. 4 c) upon Lys¹³ binding. Fig. 4, a and b, show local lipid fractions and for the ternary system, and Fig. 4 c shows a similar plot for in the binary mixture. Snapshots are taken after 500 ns starting from the homogenous distribution, and the green shades (corresponding to * = 1) represent locations on the membrane where the lipid distributions are unaffected by the adsorbing peptide. Darker colors (* > 1) show areas with excess lipids, and lighter colors (* < 1) identify locations with deficiency in the corresponding lipid species.

View larger version (35K): [in this window] [in a new window]   FIGURE 4  Adsorption of Lys13 onto ternary (PC/PS/PIP2) lipid membrane with (a and b), and onto binary (PC/PS) lipid membrane with (c). (a) Normalized local fraction of PIP2 lipids in the ternary system. (b) Local PS lipid fractions in the ternary system after 0.5 µs. (c) Local PS lipid fraction in the binary mixture.

For comparison, Fig. 4 c shows that even in PIP₂-free membranes, the segregation of PS lipids around Lys¹³ is marginal. There are regions on the membrane where the fraction of PS has increased only approximately twofold. But mostly the plot shows very weak PS sequestration compared to that seen for PIP₂.

To better assess the degree of multivalent lipid segregation around the polypeptide, in Fig. 5 we plot the electrostatic potential isosurfaces for Lys¹³ adsorbed on a ternary (PC/PS/PIP₂) lipid membrane. Fig. 5, a and b, show side and top views, respectively, of the system in the initial configuration, with homogenous lipid distribution on the membrane, and Fig. 5, c and d, provide similar views of the final state of the system, after 500 ns. Van der Waals radii for Lys¹³ are colored in gray. Isocontours are for = –1.5 k_BT/e (–37.5 mV) (red surface), and = +1.5 k_BT/e (+37.5 mV) (blue mesh). For clarity, the lipid membrane is not shown. Fig. 5 reveals significant change in the electrostatic potential near Lys¹³ after 500 ns of dynamics. Comparing Fig. 5, a and b (or Fig. 5, c and d), illustrates the growth of the negative electrostatic potential isosurface close to Lys¹³ as PIP₂ lipids segregate around the polypeptide (Fig. 4 a). This increase is because accumulation of PIP₂ near Lys¹³ reduces gradients in the local electrochemical potential, which, in turn, is proportional to the electrostatic potential. Similar plots for Lys¹³-PC/PS membrane (data not shown) reveal much weaker electrostatic potential near the adsorbate, indicative of lower charge accumulation around Lys¹³.

View larger version (68K): [in this window] [in a new window]   FIGURE 5  Electrostatic potential isosurfaces for Lys13 adsorbed on a ternary (PC/PS/PIP2) lipid membrane. (a and b) Side and top views, respectively, of the system in the initial configuration. (c and d) Similar views for the final state of the system, after 500 ns. Lys13 van der Waals surfaces are colored in gray. The red surface represents = –1.5 kBT/e (–37.5 mV) equipotential contour, and the blue mesh depicts the = +1.5 kBT/e (+37.5 mV) equipotential contour. For clarity, the lipid membrane is not shown.

Diffusion of macroion on binary and ternary lipid membranes
Above we showed that stationary basic macromolecules will preferentially sequester PIP₂ lipids. Because this picture might change if the adsorbate is allowed to diffuse, we considered the effect of protein mobility. By introducing this degree of freedom we allow the system access to states around the mean field free-energy minimum, and with that we pose several subtly related questions: how are the macromolecule diffusion rates affected by the acidic lipids in the membrane, and how will different lipids influence the apparent protein diffusion rates? To address these questions, we performed a set of calculations in which a (model) spherical macroion (R_p = 10 Å radius and uniform surface charge of 1e^–/93 Ų kept at h = 3 Å from the membrane) was allowed to move concomitantly with lipid diffusion. By performing CHDMC simulations, as detailed in Model, we studied mixed membranes of two different lipid compositions: binary mixtures with and ternary mixtures with and

For protein mobility, we focus on two typical cases. In the first, the model protein has a diffusion constant much faster than that of lipids in the unperturbed (bare) membrane, D' D_prot/D_lip = 10, while in the second the diffusion constant is comparable to that of the lipids, D' = 2. As we show, these two scenarios lead to different lipid and protein diffusion characteristics.

Modeling a fast protein diffusing over binary versus ternary membranes
Using color grade, Fig. 6 c shows the local lipid fractions after 0.6 µs of dynamic time evolution. The model protein's trajectory for the entire time interval is shown as a black line that follows the protein projected center-of-mass. The dotted red line indicates the size of the macroion projected onto the membrane, with the arrow indicting the initial protein position. For clarity, we zoom in on the membrane surface region explored by the macroion. The entire trajectory can also be found as an animation file in Supplementary Material, file No. TM-10D.avi.

View larger version (61K): [in this window] [in a new window]   FIGURE 6  Diffusion of charged spherical macroion on mixed membranes. The panels show the local surface charge densities after 0.6 µs of simulations (color scale) and the entire macroion trajectories in that time (connected black lines) for binary (PC/PS) mixture, D' = 10 (a); for binary (PC/PS) mixture, D' = 2 (b); for ternary (PC/PS/PIP2) mixture, D' = 10 (c); and for ternary (PC/PS/PIP2) mixture, D' = 2 (d). The red-dashed circles on each panel represent the projected size of the macroion with black arrows indicating the starting position of macroion center of mass. For clarity, the figures zoom in on the relevant membrane surface region explored by the macroion.

The free diffusion of the macroion does not last too long, because PIP₂ lipids quickly segregate again around the new position. This segregation is due to the large forces acting on the PIP₂ lipids by the electrostatic field emanating from the adsorbate. Concomitantly, the local lipid composition in the region of the membrane just abandoned by the macroion is restored to that of the homogenous mixture. Essentially, the macroion diffuses and drags PIP₂ lipids along the way, while the PIP₂ units that are segregated retard the free diffusion of the protein. We found that PS segregation for ternary mixtures is very weak, in accordance with our previous findings (Fig. 4 b). Therefore, we show here only the local changes in PIP₂ lipids.

We contrast these conditions with the same rapid model protein (D' = 10) diffusing on a binary membrane containing monovalent (PS) lipids. The corresponding trajectory for the binary mixture after 0.6 µs is shown in Fig. 6 a, and also in Supplementary Material, file No. BM-10D.avi. The color grade here describes the local PS fraction Clearly, acidic (PS) lipids segregate around the macroion to a much lesser extent compared to the ternary mixture, resulting in low energetic barriers to adsorbate motion. Hence, the diffusion of the macroion here is less restricted compared to that seen for the ternary mixture.

Slow protein diffusing over binary versus ternary membranes
Diffusion of a slower model protein, with D' = 2, on the same binary and ternary membranes (Fig. 6, b and d, respectively, and also in Supplementary Material, files No. BM-2D.avi and No. TM-2D.avi) shows qualitatively similar behavior to that observed for D' = 10. However, due to the lower mobility of the macroion, the acidic lipids have more time to segregate near the adsorbate, and therefore segregate more strongly. The result is that the majority of macroion moves are restricted to the acidic lipid-rich patch that forms close to the protein. This is particularly noticeable for the ternary system, where the macroion practically never escapes to go beyond the circular patch formed by PIP₂ lipids, but rather diffuses within it. Whereas for the fast protein on ternary mixtures we observe creation and destruction of macroion/PIP₂ "binding sites", for the slower protein this lipid-protein "complex" stays intact for the entire trajectory. In a sense there are always PIP₂ lipids associated with the macroion as it diffuses on the membrane.

Diffusion analysis
We now turn to the quantitative analysis of these simulation results. The final snapshots in Fig. 6 show that, on the simulation timescales, the adsorbate explores a more extended region of space when it diffuses on binary, rather than ternary membrane mixtures. To quantify this finding, we calculated the mean-square-displacement (MSD) of the macroion as a function of time for the systems shown in Fig. 6. Fig. 7 shows the MSD for a given time lag t, obtained by averaging over all pairs of points t time-steps apart (52,53). For each system, we only analyzed the last 400 ns of the trajectories. In Fig. 7, we plot for the first 150 ns of the productive runs, as sampling becomes poor for longer time lags (52,53). If we assume that, in all cases, the relationship between the MSD and elapsed time is linear, we can derive the effective (observed) protein diffusion constant from a linear regression analysis of the data. For all simulations, the diffusion coefficient of membrane-bound macroion, D_bound, was calculated first, and then we obtained D'_bound = D_bound/D_lip, which describes the draglike effect that lipids can have on the protein. The results are summarized in the table inset shown in Fig. 7.

View larger version (9K): [in this window] [in a new window]   FIGURE 7  Mean-square-displacement (MSD) as a function of time for D' = 10 (upper panel) and D' = 2 (lower panel). The table inset shows the apparent macroion and lipid diffusion coefficient ratios for free and for membrane-bound macroion, the latter calculated from linear regression analysis of the MSD plots.

Why do multivalent lipids influence the adsorbate diffusion rates more than monovalent lipids? And what role does the inherent (or free) diffusion rate of the macroion (diffusion rate on homogeneous surfaces) play in this process? Assume that the adsorbate inherently travels with high diffusion rate, D' = 10. Then, only highly mobile lipids would be able to segregate around the protein. Our simulations of ternary mixtures show that PIP₂ can be fast enough. These multivalent lipids experience strong electrostatic forces from the adsorbate because of their strongly charged headgroups and therefore are quickly sequestered near the macroion. Once sequestered, PIP₂ lipids act to confine the motion of the adsorbate, as any protein motion leaves behind an exposed charged patch that is highly unfavorable. The protein thus finds it difficult to overcome and escape this electrostatic-well created by sequestered PIP₂ molecules (8,10,49), and the retarded macroion has a lower diffusion rate.

In contrast, monovalent (PS) lipids segregate only weakly around a quick macroion. Hence, the energetic barriers created by these sequestered acidic lipids are low and the adsorbate, due to its high mobility, finds it easy to consistently escape them. The preferential interaction coefficients show on average –1.2e membrane surface charges for the fast protein on ternary mixtures versus –1.0e in the binary mixture accumulated near the macroion. Furthermore, our calculations reveal that PIP₂ and PS lipids contribute with their headgroup charges almost equally in the ternary mixtures, –0.7e charge coming from PIP₂ and –0.5e from PS. Based on these findings and in light of substantial difference in mono- and multivalent lipid composition, we conclude that localized electronic charge on multivalent PIP₂ lipid plays a major role in regulating the macroion diffusion rate.

As the free macroion diffusion rate decreases, monovalent and multivalent lipids tend to affect adsorbate diffusion rates to a similar extent, as revealed from the MSD plots for D' = 2, Fig. 7. This is because the underlying acidic lipids are able to rearrange around such slow-diffusing adsorbates, resulting in lipid segregation rather similar to that observed in equilibrium for the immobile protein. Under these conditions, even the less mobile PS lipids can retard the adsorbate's diffusion. Conversely, since the protein is slow enough so that all lipids can relax for all protein steps, the relative drag experienced by the protein and its slowing down are smaller than for fast proteins.

From Figs. 6 and 7 it is evident that there is some level of confinement with possible binding and unbinding events of interacting lipid and proteins, and hence the MSD plots need not be interpreted as strictly linear. In fact, the higher the protein diffusion rate D', the more anomalous the observed macroion diffusion on the ternary membranes. In the limit of the unrealistically high D' = 50, our simulations (data not shown) show adsorbate motion on ternary membranes that includes extended local motions with occasional rapid diffusion, resembling a hop-diffusion-type mechanism (52,53). Importantly, even for such extreme D', multivalent lipids are still capable of affecting protein diffusion rate. On the other hand, for large D' the diffusion on the binary mixtures is still close to linear, and for D' = 50 we see practically no slowing-down of the macroion diffusing on a binary PC/PS membrane. Obviously, in the limit of protein and lipid motion will become largely uncorrelated regardless of whether the protein is diffusing on ternary or binary membrane.

Because the DMC method allows thermal fluctuations on the order of k_BT away from the mean-field free energy, it is interesting to follow these changes for different membrane-protein interactions. In Fig. 8, a and b, we plot the instantaneous adsorption free energies (in thermal energy units) given by Eq. 4 as the macroion diffuses on binary and ternary mixtures, respectively. In both panels, the horizontal lines (D' = 0) represent the calculated equilibrium mean-field binding free energies for the respective stationary macroion and membrane complexes (no protein motion is allowed), measured with respect to the unbound protein and bare membrane in solution. However, lipids are still mobile and are allowed to segregate around the stationary macroion.

View larger version (24K): [in this window] [in a new window]   FIGURE 8  Value of the instantaneous adsorption free-energy functional F (in kBT units) as a function of time in the CHDMC simulations for (a) D' = 2 systems and for (b) D' = 10. The horizontal lines show the calculated equilibrium adsorption free energies for the respective systems when the macroion is stationary (see text for details).

Note that the largest possible difference in binding free energies between any two states of the mobile protein/membrane complex is 1.5 k_BT, as expected for thermal fluctuations. In one state the membrane is homogenous, and the lipids are ideally mixed (adsorption free energy at t = 0 in Fig. 8, a and b), and in the other state lipids are completely relaxed around the immobile adsorbed protein (D' = 0 blue lines in Fig. 8, a and b). Because a stationary macroion sequesters charged lipids most efficiently, the equilibrium-binding free energies for the stationary macroion/membrane complexes in Fig. 8 provide a lower-bound for binding free energies of the macroion/membrane systems. The plots reveal that the adsorption free energies for the ternary mixtures are close to the lowest possible values at multiple timepoints in the simulation, indicating that at those instances the mobile protein sequestered multivalent lipids to an extent similar to a stationary protein. The picture is much different when the macroion diffuses on mixed binary membranes: the adsorption free energies never reach the limiting value, because lipids are too slow to relax locally around the protein.

	CONCLUDING REMARKS

TOP ABSTRACT INTRODUCTION MODEL RESULTS AND DISCUSSION CONCLUDING REMARKS APPENDIX SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
We have described a self-consistent dynamic mean-field model to study the process of adsorption of charged macromolecules onto oppositely charged lipid membranes. Extending earlier thermodynamic calculations (20), our model adds the dynamic aspects of lipid lateral reorganization and demixing when adsorbed charged macromolecules diffuse upon membranes. In particular, to our knowledge this is a first attempt to study the kinetic mechanisms of acidic monovalent and multivalent lipid segregation around model basic proteins.

The model allows us to solve the full three-dimensional electrostatic problem within the nonlinear Poisson-Boltzmann theory rapidly and efficiently, and at the same time to use th