This Article Abstract Full Text (PDF) 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 Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Venturoli, M. Articles by Sperotto, M. M. Search for Related Content PubMed PubMed Citation Articles by Venturoli, M. Articles by Sperotto, M. M.

Simulation Studies of Protein-Induced Bilayer Deformations, and Lipid-Induced Protein Tilting, on a Mesoscopic Model for Lipid Bilayers with Embedded Proteins

Maddalena Venturoli ^* , Berend Smit ^* and Maria Maddalena Sperotto 

^* Department of Chemical Engineering, University of Amsterdam, Amsterdam, The Netherlands; Department of Chemistry, University College, London, United Kingdom; and Biocentrum, The Technical University of Denmark, Kgs. Lyngby, Denmark

Correspondence: Address reprint requests to Maria Maddalena Sperotto, E-mail: maria{at}cbs.dtu.dk.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MODEL AND SIMULATION METHOD RESULTS AND DISCUSSION CONCLUSION AND FUTURE... ACKNOWLEDGEMENTS REFERENCES

 
Biological membranes are complex and highly cooperative structures. To relate biomembrane structure to their biological function it is often necessary to consider simpler systems. Lipid bilayers composed of one or two lipid species, and with embedded proteins, provide a model system for biological membranes. Here we present a mesoscopic model for lipid bilayers with embedded proteins, which we have studied with the help of the dissipative particle dynamics simulation technique. Because hydrophobic matching is believed to be one of the main physical mechanisms regulating lipid-protein interactions in membranes, we considered proteins of different hydrophobic length (as well as different sizes). We studied the cooperative behavior of the lipid-protein system at mesoscopic time- and lengthscales. In particular, we correlated in a systematic way the protein-induced bilayer perturbation, and the lipid-induced protein tilt, with the hydrophobic mismatch (positive and negative) between the protein hydrophobic length and the pure lipid bilayer hydrophobic thickness. The protein-induced bilayer perturbation was quantified in terms of a coherence length, _P, of the lipid bilayer hydrophobic thickness profile around the protein. The dependence on temperature of _P, and the protein tilt-angle, were studied above the main-transition temperature of the pure system, i.e., in the fluid phase. We found that _P depends on mismatch, i.e., the higher the mismatch is, the longer _P becomes, at least for positive values of mismatch; a dependence on the protein size appears as well. In the case of large model proteins experiencing extreme mismatch conditions, in the region next to the so-called lipid annulus, there appears an undershooting (or overshooting) region where the bilayer hydrophobic thickness is locally lower (or higher) than in the unperturbed bilayer, depending on whether the protein hydrophobic length is longer (or shorter) than the pure lipid bilayer hydrophobic thickness. Proteins may tilt when embedded in a too-thin bilayer. Our simulation data suggest that, when the embedded protein has a small size, the main mechanism to compensate for a large hydrophobic mismatch is the tilt, whereas large proteins react to negative mismatch by causing an increase of the hydrophobic thickness of the nearby bilayer. Furthermore, for the case of small, peptidelike proteins, we found the same type of functional dependence of the protein tilt-angle on mismatch, as was recently detected by fluorescence spectroscopy measurements.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MODEL AND SIMULATION METHOD RESULTS AND DISCUSSION CONCLUSION AND FUTURE... ACKNOWLEDGEMENTS REFERENCES

 
Biological membranes are complex, organized, dynamic, and highly cooperative structures whose physical properties are important regulators of vital biological functions ranging from cytosis and nerve processes, to transport of energy and matter (Sackmann, 1995). To relate the structure and dynamics of biomembranes to their biological function (the ultimate goal of biomembrane science), it is often necessary to consider simpler systems. Lipid bilayers composed of one or two lipid species with embedded proteins, or natural or artificial peptides, provide a model system for biological membranes. Understanding the physics of such simplified soft-condensed matter systems can yield insight into biological membrane functions. Therefore these systems are extensively investigated, both experimentally and theoretically.

The hydrophobic matching between the lipid bilayer hydrophobic thickness and the hydrophobic length of integral membrane proteins has been proposed as one of the main physical mechanisms that regulate the lipid-protein interaction in biomembranes (Mouritsen and Blom, 1984; Sackmann, 1984; Mouritsen and Sperotto, 1993; Gil et al., 1998; Killian, 1998; Dumas et al., 1999). The energy cost of exposing polar moieties, from either hydrocarbon chains or protein residues, is so high that the hydrophobic part of the lipid bilayer should match the hydrophobic domain of membrane proteins. The results from a number of investigations have indeed pointed out the relevance of the hydrophobic matching in relation to the lipid-protein interactions, hence to membrane organization and biological function. It is now known that hydrophobic matching is used in cell membrane organization: the membranes of the Golgi have different thicknesses; along their secretory pathway, proteins that pass through the Golgi undergo changes of their hydrophobic length to match the membrane hydrophobic thickness of the Golgi (Munro, 1995,1998; Bretscher and Munro, 1993; Pelham and Munro, 1993). Hydrophobic matching seems also to play a role in sequestering proteins with long transmembrane regions (McIntosh et al., 2003) into sphingolipids-cholesterol biomembrane domains denoted rafts (Simons and Ikonen, 1997; Binder et al., 2003). The biological importance of rafts, and their involvement in Alzheimer's and prion diseases, is nowadays an intensively investigated subject (Fantini et al., 2002).

Biological membranes have at their disposal a number of ways to compensate for hydrophobic mismatch (de Planque and Killian, 2003), which may be used individually or simultaneously. These ways may imply changes of the membrane structure and dynamics on a microscopic, as well as on a macroscopic scale, and therefore can affect the membrane biological function (Montecucco et al., 1982; Johansson et al., 1981; In't Veld et al., 1991; Lee, 1998). To adjust to hydrophobic mismatch a membrane protein may cause a change of the lipid bilayer hydrophobic thickness in its vicinity. Experimental studies on reconstituted systems show that the range of the perturbation induced by proteins on the membrane thickness varies considerably from system to system (Jost et al., 1973; Hesketh et al., 1976; Jost and Hayes Griffith, 1980; Rehorek et al., 1985, Piknová et al., 1993; Harroun et al., 1999; Bryl and Yoshihara, 2001). A lipid sorting at the lipid-protein interface may also occur, where the protein prefers, on a statistical basis, to be associated with the type of lipid that best matches its hydrophobic surface (Dumas et al., 1997; Lehtonen and Kinnunen, 1997; Fahsel et al., 2002; Fernandes et al., 2003). Another way that a protein may have to adapt to a too-thin lipid bilayer is to tilt (Glaubitz et al., 2000; Killian, 1998; Sharpe et al., 2002; van der Wel et al., 2002; Koehorst et al., 2004; Strandberg et al., 2004). In addition to the protein as a whole, the individual helices of which a protein might be composed may also experience a tilt; and there is indeed some experimental evidence that the latter phenomenon may occur in channel proteins (Lee, 2003), and that a change in tilt-angle of the individual helices could be the cause of a change in protein activity. Long and single-spanning membrane proteins might also bend to adapt to a too-thin bilayer. Spectroscopic measurements on phospholipid bilayers, with embedded poly(leucine-alanine) -helices, suggest that the conformation of long peptides deviates from a straight helical end-to-end conformation (Harzer and Bechinger, 2000; Strandberg et al., 2004). A protein may also undergo structural changes to adapt to a mismatched lipid bilayer. Spectroscopy measurements indicate that, indeed, long hydrophobic polyleucine peptides might distort in the C- and N-terminus to reduce their hydrophobic length and thus match the thickness of the lipid bilayer in the gel-phase (Liu et al., 2002). Lipid-mediated protein aggregation could also occur to reduce the stress caused by hydrophobic mismatch (Mall et al., 2001; Fernandes et al., 2003; Harroun et al., 1999). Domains (Binder et al., 2003) may thus form whose functional properties differ from those of the bulk, i.e., the unperturbed bilayer (Tocanne, 1992; Tocanne et al., 1994; Thomson et al., 1995). When the degree of mismatch is too high to be counterbalanced by the adaptations just described, the proteins might partition between an in-plane and a transmembrane orientation, or even avoid incorporation in the membrane (Harzer and Bechinger, 2000; de Planque et al., 2001; Ridder et al., 2002). The phenomena just mentioned refer to local microscopic changes related to mismatch adjustment. Perturbations of the membrane on the macroscopic scale may also occur; these can range from in-plane protein segregation and crystallization, and gel-fluid phase separation (Mouritsen, 1998; Gil et al., 1998; Dumas et al., 1997; Morein et al., 2002; Fahsel et al., 2002), to changes of the three-dimensional structure of the membrane. The formation of nonbilayer phases upon protein incorporation in lipid bilayers is an example of the latter type of phenomena (Killian, 1992; Epand, 1998).

In an effort to elucidate the effects caused at the molecular level by the lipid-protein hydrophobic mismatch, and even their possible implications for the formation of biologically relevant domainlike structures such as rafts, a number of theoretical studies have been done with the help of different types of theoretical models (Sperotto and Mouritsen, 1991; Fattal and Ben-Shaul, 1993; Mouritsen et al., 1996; Gil and Ipsen, 1997; Belohorcová et al., 1997, 2000; Sintes and Bäumgartner, 1998; Dan and Safran, 1998; Nielsen et al., 1998; May, 2000; Duque et al., 2002; Petrache et al., 2000b, 2002; Jensen et al., 2001, 2002; Shen et al., 1997; Bohinc et al., 2003; Jensen and Mouritsen, 2004). One of the quantities that has drawn considerable attention in recent years is the extension of the domain size, which is determined by the coherence length of the spatial fluctuations occurring in the system. Such fluctuations, which depend on the thermodynamic state of the system, can be induced, as well as harvested, by proteins. In the past, computer simulations have been made on a lattice model to measure the extent of the perturbation induced by a protein on the surrounding lipid bilayer (Sperotto and Mouritsen, 1991). The results from these simulations indicated that the extension of the perturbation depends on factors such as the degree of hydrophobic mismatch, the size of the protein (i.e., the curvature of the protein hydrophobic surface in contact with the lipid hydrocarbon chains), and on the temperature of the investigated system. Also, it was found that, away from the protein, the perturbation decays in an exponential manner, and can therefore by characterized by a decay length, _P. The value _P is a measure of the size of small-scale inhomogeneities (i.e., domains) experienced by proteins when embedded in the lipid bilayer. In a sense, _P is also a measure of the extension of the range over which the lipid-mediated interaction between proteins may operate. Results from model studies of a phenomenological interfacial model for proteinlike objects in a bilayerlike system suggest that, under well-defined thermodynamic conditions, the protein-induced perturbation may propagate without decay over a number of lipid shells around the protein (the number of lipid shells being dependent, among other factors, on the size of the protein), may extend over long ranges, and might eventually establish a thermodynamic phase (Gil and Ipsen, 1997; Gil et al., 1998). The phase of the multilayered region that the protein prefers to be surrounded with is thus said to wet the protein (Gil and Ipsen, 1997; Gil et al., 1998). The disadvantage of using models such as the lattice models (Sperotto and Mouritsen, 1991) or the phenomenological models (Gil and Ipsen, 1997; Gil et al., 1998; Dan and Safran, 1998; Nielsen et al., 1998) is that these models do not allow for tilting of the model proteins as a whole. Therefore, with the help of these models, one cannot make predictions about the physical hydrophobic-mismatch condition that induce a protein to tilt in the lipid bilayer, other than, at the same time, by inducing a bilayer deformation in the vicinity of the protein. With the help of a microscopic model, Duque et al., (2002) have indeed studied how hydrophobic mismatch affects the way in which the inclination of transmembrane helices changes as a function of their hydrophobic length. Their self-consistent calculations predicted that peptides, whose hydrophobic length is less than that of the hydrophobic bilayer thickness, insert perpendicular to the bilayer—whereas peptides with a longer hydrophobic length than the bilayer hydrophobic thickness insert into the bilayer in a tilted manner, and with an angle with the bilayer-normal that increases with increasing mismatch.

The types of models described above are relatively crude, in the sense that either they cannot be used to investigate the physical conditions that cause protein tilting, or they do not take into full account the three-dimensional molecular structure of the bilayer. Simulations on more realistic models, such as all-atom models for lipid bilayers with embedded proteins, have confirmed that, at least within a time of the order of the nanoseconds, a mismatched protein can induce a deformation of the lipid bilayer structure (Chiu et al., 1999; Petrache et al., 2000b, 2002; Jensen et al., 2001), and that the deformation is of the exponential type (Jensen and Mouritsen, 2004). The same type of studies have also shown that tilting may also occur for membrane peptides (Belohorcová et al., 1997; Shen et al., 1997); however, to reduce a possible hydrophobic mismatch, synthetic peptides might instead prefer to deform the lipid bilayer, rather than undergo tilting (Petrache et al., 2002). Incidentally, the results from these studies indicated that the helical-peptides experience a slight bend in the middle of the helix.

Regardless of the huge body of experimental and theoretical studies on lipid bilayers with embedded proteins, issues like the range of the protein-induced lipid bilayer perturbation, its dependence on protein size, and the simultaneous occurrence of protein tilting (or even bending) to adjust for hydrophobic mismatch, are still a matter of debate. In this article we want to focus on these issues, which we investigate by means of a mesoscopic model for lipid bilayers with embedded proteins and a relatively new simulation method. Before introducing the model and the method, we would like to sketch the historical background that brought scientists to the development of mesoscopic models to study physical phenomena in biomembranes.

Because of the many degrees of freedom involved, the processes that take place even in model biomembranes occur over a wide range of time- and lengthscales (König and Sackmann, 1996). To model membranes, it is thus necessary to decide, a priori, the level of description of the system (i.e., to deliberately neglect those details unimportant to the process one wants to investigate). Often, this necessity follows the fact that some theoretical methods are limited in their applicability by the long computational times needed to calculate statistical quantities. The drawbacks of those otherwise relatively noncomputationally-time-demanding phenomenological, lattice, or interfacial models have outlined the necessity of more realistic models to investigate the effect of proteins on the bilayer structure and dynamics. Molecular dynamics (MD) simulation methods on all-atom models have been used to study the self-assembly of phospholipids into bilayers (Marrink et al., 2001) as well as the structure, dynamics, and interactions of individual membrane peptides or proteins with the lipid bilayer (Shen et al., 1997; Belohorcová et al., 1997, 2000; Jensen et al., 2001; Jensen and Mouritsen, 2004). MD simulations can provide detailed information about the phenomena that occur in biomembrane systems, although at the nanoscopic level and on a nanosecond timescale. Many membrane processes happen at mesoscopic length- and timescales, however—that is, above 1–1000 nm, and 1–1000 ns, respectively—and involve the collective nature of the system. This is the case for phenomena related to the gel-fluid phase transition and phase separation, the formation of domains on the mesoscopic scale, or the transition from a bilayer to a nonbilayer phase. Even though the speed of numerical computation is increasing very rapidly, it will be some time before it will be possible, by MD on realistic all-atom models, to predict the cooperative behavior of biosystems at mesoscopic timescales. Numerical studies of these phenomena require a considerable simplification of the model. These simplifications can be made by using a system of particles, or beads, in which each particle represents a complex molecular component of the system whose details are not important to the process under investigation. These models with simplified interactions between the beads are called coarse-grain (CG) or mesoscopic models. In recent years, CG models have been developed to study the phase equilibria of biomembrane-like systems at the mesoscopic level, and both MD and Monte Carlo (MC) simulation methods were used on such models (Goetz and Lipowsky, 1998). With the minimal modeling approach it was possible to simulate the self-assembly of phospholipids into various phases, both in the absence and presence of such biologically relevant molecules as anesthetics and alkanes (Shelley et al., 2001a,b). It was also possible to study the lipid-mediated range of attraction between two proteins embedded in a lipid bilayer (Sintes and Bäumgartner, 1998).

Despite the advantages that arise by minimal modeling in connection with simulation methods like MD and MC, the possibility to study processes that involve the collective behavior of the system is still limited. To try to overcome this limitation, the use of a faster simulation technique, i.e., dissipative particle dynamics (DPD), on CG models has thus been considered. The DPD-on-CG-model approach can be seen as a middle course between the approach based on pseudo-three-dimensional models (such as lattice and interfacial models) and the one based on all-atom models. The DPD method was originally developed to simulate complex fluids, such as surfactant and polymer melts, at the mesoscopic level. It was then adopted to study mesoscopic models for pure lipid bilayer systems (Venturoli and Smit, 1999), as well as lipid bilayers containing impurities such as alcohols (Kranenburg and Smit, 2004; Kranenburg et al., 2004b). The results from the simulation studies demonstrated that with the DPD-CG approach one was able to reproduce the structural and thermodynamic properties resulting from cooperative behavior of the lipid bilayer system (Kranenburg et al., 2003a,b).

We have adopted the DPD simulation method to study the behavior of a mesoscopic model for lipid bilayers with embedded proteins. The model was derived from the mesoscopic model for pure lipid bilayers (Venturoli and Smit, 1999) and was extended to account for the presence of proteins. We studied systems at low protein/lipid ratios.

The aim of the work presented in this article was to understand whether, and to what extent due to hydrophobic mismatch, and via the cooperative nature of the system, a protein may prefer to tilt (with respect to the normal of the bilayer plane), rather than to induce a bilayer deformation without (or even with) tilting. Therefore we have attempted to make a systematic correlation of the protein-induced perturbation and lipid-induced protein tilting with hydrophobic mismatch.

The article is structured as follows. First we describe the mesoscopic model for lipid bilayers with embedded proteins, and present the DPD simulation method. We then present the model parameters, the statistical ensemble used for the simulations, and the methods of calculation of statistical quantities. In Results and Discussion, data for both pure lipid bilayers and those with embedded proteins are shown and discussed. Whenever possible, we have validated the model by comparing the results obtained from our model study with existing theoretical and experimental data. The results from the simulation studies and the model predictions are summarized in Conclusion and Future Perspectives, together with possible future applications of the DPD-on-CG-model approach.

	MODEL AND SIMULATION METHOD

TOP ABSTRACT INTRODUCTION MODEL AND SIMULATION METHOD RESULTS AND DISCUSSION CONCLUSION AND FUTURE... ACKNOWLEDGEMENTS REFERENCES

 
Mesoscopic model
Within the mesoscopic approach, each molecule of the system (or groups of molecules) is coarse-grained by a set of beads. In the specific case of the lipid-protein system, we considered three types of beads: a waterlike bead, labeled w; a hydrophilic bead, labeled h, which models a part of the headgroup of either the lipid or the protein; and a hydrophobic bead, labeled either t_L or t_P, depending on whether it refers to a portion of the lipid hydrocarbon chain or a portion of the hydrophobic region of the protein, respectively. Each model-lipid is built by one or more headgroup-like h-beads connected to two tails of equal length. Each of these tails is formed by connecting with springs a chosen number of t_L-beads, depending on the type of lipids one wants to model. Fig. 1 a shows a schematic representation of a model-lipid. A w-water-bead represents three water molecules, and a t-bead represents three CH₂ groups (or one CH₂ plus one CH₃ group) of the lipid hydrocarbon chain (Kranenburg et al., 2004a). The systems that we have simulated are made of model-lipid chains having three headgroup beads and five beads in each chain; this corresponds to the case of an acyl chain with 14 carbon atoms, namely to a model for a dimyristoylphosphatidylcholine phospholipid, as illustrated in Fig. 2.

View larger version (58K): [in this window] [in a new window]   FIGURE 1  Schematic representation of a model lipid (a), and a model protein (NP = 43 and ) (b). A typical configuration of the assembled bilayer with a model protein embedded (as results from the simulations) is shown in the snapshot (c). The drawing in d shows the part of the system to which the following quantities refer: the pure lipid bilayer hydrophobic thickness, ; the perturbed lipid bilayer hydrophobic thickness, dL(r); the protein hydrophobic length, dP; the tilted-protein hydrophobic length, ; and the tilt-angle, tilt.

View larger version (18K): [in this window] [in a new window]   FIGURE 2  The atomistic representation of DMPC and its corresponding coarse-grained model. Hydrophilic head-beads are indicated in shading and hydrophobic tail-beads in open representation.

Dissipative particle dynamics
We studied the mesoscopic model with the help of the dissipative particle dynamics (DPD) simulation method (Hoogerbrugge and Koelman, 1992; Warren, 1998; Jury et al., 1999). The DPD method was originally based on the idea of simulating the fluid hydrodynamics of systems composed of particles, or beads, in analogy with the way the Navier-Stokes equations reproduce the motion of a real fluid. Each bead, which represents the center of mass of a small droplet of the fluid, moves according to Newton's equation of motion, and interacts according to simplified force laws. The beads interact with each other via conservative, random, and dissipative forces of the pairwise-additive type. The total force, f_i, acting on bead i, is thus expressed as a sum over all other beads, j, which are within a certain cutoff radius R_c from bead i,

(1)

The first term in Eq. 1 refers to a force of conservative type. This comprises two contributions, one related to interactions between beads not bound together, and the other related to interactions between beads that are linked together. The former contribution is chosen in such a way to model a soft-repulsive potential,

(2)

where the coefficients a_ij > 0 represent the maximum repulsion strength, r_ij = r_i – r_j is the distance between bead-particles i and j, and R_c is the cutoff radius, which gives the extent of the interaction range. The conservative force can have also an elastic contribution, which derives from the harmonic force used to tie two consecutive beads in the chains of either the lipid or the protein. This contribution is expressed as

(3)

where K_r is the elastic constant, and r_eq is the equilibrium value of r_ij. To control the chain flexibility, an extra bond-bending force between consecutive bonds is added,

(4)

(5)

where K is the bending constant, is the angle between two consecutive bonds, and _o is the equilibrium angle. The other two forces in Eq. 1 are a drag force (F^D) and a random force (F^R), which are expressed as

(6)

where v_ij = v_i – v_j is the velocity difference between particles i and j, is the friction coefficient, and is the noise amplitude. The quantity _ij is a random number, which is chosen from a uniform random distribution, and in an independent manner for each pair of particles. The chosen functional dependence on r_ij of the conservative force, F^C, permits us to use larger integration time-steps than are usually allowed by the MD simulation technique (which has to do with computationally demanding forces of the Lennard-Jones hard-core type). Also, the combined effect of the two forces, the dissipative and the random, acts as a thermostat—which conserves the (angular) momentum and thus provides the correct hydrodynamics to the system, at least for sufficiently long timescales and large system sizes.

Español and Warren (1995) have shown that the equilibrium distribution of the system is the Gibbs-Boltzmann distribution, if the weight functions and coefficients of the drag and random forces satisfy

(7)

(8)

where k_B is the Boltzmann constant and T is the temperature. Furthermore, all the forces assume the same functional dependence on the interparticle distance r_ij (as the conservative force does) if the weight function w^R(r) is of the type

(9)

Model parameters
The repulsion parameter (see Eq. 2) related to the interaction between the water beads, a_ww, was derived by fitting the calculated value of the compressibility of water, at room temperature, to the experimental one (Groot and Warren, 1997). In principle, this fitting procedure may be applied at any temperature, and may thus result in temperature-dependent a_ij parameters. To deal with temperature-dependent parameters would make the interpretation of the simulation data very difficult. Therefore, we make an approximation in which we assume that the parameters a_ij are not temperature-dependent. These parameters have been chosen so as to reproduce the structural and thermodynamic behavior of the pure system, i.e., of a pure DMPC bilayer. Because a direct mapping between the atomic level information and the model parameters is not always possible, it is worth mentioning that these are effective parameters and reflect this limitation. Therefore, it will not always be possible to make a direct comparison between the properties of the model system and the properties of the reconstituted system.

The values of the parameters referring to the lipid-lipid and lipid-water interaction have been chosen equal to those used for the pure lipid bilayer model (Kranenburg et al., 2003a). To model the amphiphilic nature of the lipids, the repulsion parameters a_ij (Eq. 2) between two beads (whether hydrophilic or hydrophobic) were chosen to be smaller than the repulsion parameters between two beads of which one is hydrophilic and the other hydrophobic. The numerical values of the interaction parameters between different bead types are given in Table 1.

Concerning the elastic contribution to the interaction energy (see Eq. 3), at an overall bead density of 3 (Groot and Warren, 1997), the resulting equilibrium distance is equal to r_eq = 0.7. To determine the spring constant, K_r, for the lipid chain, we required that 98% of the bond-distance distribution be within one R_c. A value of K_r = 100 was found to satisfy this requirement.

The values of the parameters related to the bond-bending force (Eq. 4) for the model lipids were derived from MD simulations on an all-atom model for a DMPC lipid bilayer (Kranenburg et al., 2004a). The resulting values for the bending constant and the equilibrium angle in the lipid tails are K = 6 and _o = 180°, respectively. For the bond-bending potential between the head-bead connected to the lipid tails and the first beads in the tails (beads 3, 4, and 9 in Fig. 2), values of K = 3 and _o = 90° were found to reproduce the correct configurational distribution, and structure, of the all-atom model for a DMPC lipid molecule.

Compared to the lipid hydrocarbon chains, the hydrophobic part of membrane proteins can be considered fairly rigid; therefore the value of the bending constant in the protein chains was set equal to K = 100, i.e., an order-of-magnitude larger than that used for the lipid chains. The equilibrium spring-distance and spring constant between the beads in the protein chains were chosen equal to the values used for the lipid chains, i.e., r_eq = 0.7 and K_r = 100, respectively.

Length-, time-, and temperature-scales
Usually, within the DPD approach, one makes use of reduced units for the mass, length, and energy (Groot and Warren, 1997; Groot and Rabone, 2001). The DPD unit of length is the cutoff radius, R_c, the unit of mass is the mass, m, of a bead (where all the beads in the system have equal mass), and the unit of energy is k_BT. Therefore the temperature is expressed in reduced units as well. Before presenting our results, we would like to spend a few words to give an estimate of the typical time- and lengthscales (expressed in terms of physical units) involved when one uses the DPD simulation method.

For the lengthscale, one can say that the level of coarse-graining (i.e., the number, N_m, of atoms, or molecules, represented by a DPD bead), is the renormalization factor for the mapping of the reduced units of length onto physical units. To estimate the value of the cutoff radius, R_c, one can reason as follows (Groot and Rabone 2001): If a DPD bead corresponds to N_m water molecules, then a cube of volume represents N_m water molecules, where is the density, i.e., the number of DPD beads per cubic R_c. Assuming that a water molecule has approximately a volume of 30 ų, one obtains

(10)

If a bead density equal to = 3 (Groot and Warren, 1997) is chosen, the cutoff radius is then equal to

(11)

As previously detailed in Mesoscopic Model, the results discussed below refer to a coarse-graining with N_m = 3. The choice of N_m = 3 results in a DPD bead having the volume of three water molecules, i.e., 90 ų. Therefore, from Eq. 11 one obtains R_c = 6.46 Å. This value may then be used to convert the reduced units of length into Ångstrøms.

In the literature, different mapping criteria have been adopted to derive the physical unit of time, . All these criteria were based on the mapping of the experimental value of the diffusion constants of one of the components of the system onto the value obtained from the simulations. For example, Groot and Rabone (2001) considered the self-diffusion constant of water, whereas Groot (2000) used the diffusion constant of a surfactant micelle. In both cases, a value of the integration timestep, t = 0.06, was used, which gave t 5 ps and t 25 ps, respectively. Both of these values show that DPD simulations allow for a timestep that is at least three orders-of-magnitude longer than that used in atomistic MD simulations, which is typically of the order of a few femtoseconds.

To give an estimate of the values of the reduced temperatures in terms of physical temperatures, we have mapped the reduced temperatures, T*, onto physical temperatures, T, according to the linear relation

(12)

The values of the coefficients a and b were found by solving the system of linear equations obtained by substituting in Eq. 12 the reduced and physical values of the main- and pre-transition temperatures, 24°C and 13.5°C (Koynova and Caffrey, 1998), respectively, for a pure DMPC phospholipid bilayer. The resulting values are a = 133°C and b = –33°C.

Statistical ensemble and surface tension
It has been suggested (Jähnig, 1996) that unconstrained, self-assembled bilayers are at their free-energy minimum, characterized by having a zero value of the surface tension. Nevertheless, it is still a matter of debate which value of the surface tension should be used in molecular simulations. For both self-assembled and preassembled membranes, a fixed number of lipid molecules and a fixed area combined with periodic boundary conditions are generally used in MD simulations. Periodic boundary conditions minimize the effects due to the finite size of the bilayer, but the fixed size of the simulation box imposes a constraint on the bilayer, which might also result in a finite surface tension. Although the constraint on the fixed area can be released by performing simulations at constant pressure or constant surface tension (Chiu et al., 1995, Zhang et al., 1995), there is still the question of which value of the surface tension should be used to reproduce the area per lipid of a simulated lipid bilayer. In their molecular dynamics simulations, Feller and Pastor (1996, 1999) observed that a tensionless state did not reproduce the value of the area per lipid derived from experiments. They argued that this is because the typical undulations and out-of-plane fluctuations of a macroscopic membrane cannot develop in a patch of a membrane, whose size is similar to that considered during MD simulations. They concluded that a positive surface tension (stretching) must be imposed on the system to compensate for the suppressed undulations, and to be able to reproduce the value of the area per lipid calculated from experiments. However, more recently, Marrink and Mark (2001) investigated the system size-dependence of the surface tension in large membrane patches ranging from 200 to 1800 lipids, simulated for times up to 40 ns. These authors found that simulations at zero surface tension correctly reproduce the experimental surface areas for an unstressed membrane. Goetz et al. (1998) suggested performing simulations for the exact area at which the interfacial tension is 0, and to determine this area iteratively. We have adopted a different approach, in which we mimic the experimental condition by simulating a system in which we impose a value of the surface tension. To be able to impose a given value of the surface tension on the model system, we have adopted a hybrid scheme based on both the DPD and the Monte Carlo (MC) simulation methods. The DPD method was used to evolve the positions of the beads; the Newton's equations of motion were integrated by adopting a modified version of the velocity Verlet algorithm (Groot and Warren, 1997). The MC method was used to impose a given surface tension on the bilayer. This was done by changing the bilayer projected area on the plane perpendicular to the bilayer normal, A, by an amount, A, and, at the same time, by changing the height of the simulation box to ensure that the total volume of the system remains constant (Venturoli and Smit, 1999), and therefore no work is done against the external pressure. The MC acceptance probability, P_acc, was expressed as

(13)

where U and U' define the energies before and after the change of the box sizes, respectively, and where ß = 1/k_BT. To obtain the tensionless state of the bilayer, was set to zero in Eq. 13.

Before collecting the data used to estimate the statistical quantities of interest, we have first equilibrated each bilayer system for 20,000 DPD-MC cycles. In each cycle it was chosen, with a probability of 70%, whether to perform a number of DPD steps, or to attempt to change the box aspect-ratio according to the imposed value of the surface tension, = 0. After equilibration, data were collected over 50,000 DPD-MC cycles, at = 0. The statistical averages of the quantities of interest (see Method of Calculation of Statistical Quantities, below) were then made over configurations, which were separated from one another by 50 DPD steps. On average, 10,000 independent configurations were considered for the calculation of each of the statistical averages.

Method of calculation of statistical quantities
We have studied the physical properties of the model system both in the absence and in the presence of the proteins. The pure lipid bilayer hydrophobic thickness, was estimated by calculating the difference between the average position along the bilayer normal (i.e., the z direction, if one considers the bilayer parallel to the x,y plane) of the tail-beads attached to the headgroup (beads 4 and 9, as illustrated in Fig. 2) of the lipids in one (top) monolayer, and of the lipids in the opposite (bottom) monolayer,

(14)

where z_t is the z position of either bead 4 or 9 (Fig. 2) of the lipid. The overline indicates an average over the two chains for each lipid, and over the total number of lipids in each monolayer. The difference between the two terms in the above expression is further averaged over the number of the ensemble configurations.

To study the effect of a protein on the surrounding bilayer structure, we have calculated the lipid-bilayer hydrophobic thickness, d_L(r), as a function of the radial distance r from the protein hydrophobic surface, that is, at the interface with the lipid hydrocarbon chains, as schematically illustrated in Fig. 1 d. The method of calculation of d_L(r) resembles the one used to calculate as illustrated in Fig. 3. For each configuration, we have first calculated the circularly averaged value of the positions along the bilayer normal of the tail-beads attached to the headgroup of the lipids within each circular sector k (k = 1,2,3,...) at distance r = kr from the protein surface. The bin size r was chosen to be of the order of the diameter of the lipid projected area on the bilayer plane. This was done for both monolayers of the bilayer. The instantaneous value of the bilayer hydrophobic thickness at distance r from the protein surface is then given by the difference of these two values. To obtain d_L(r), this difference has been further averaged over all the sampled configurations, as

(15)

View larger version (23K): [in this window] [in a new window]   FIGURE 3  Schematic drawing to illustrate the method of calculation of dL(r), which is described in detail in the text. The protein is represented by a shaded cylinder. The figure shows the case when the protein is parallel to the bilayer normal (a), and the case when the protein is tilted with respect to the bilayer normal (b).

The behavior of d_L(r) allowed us to access the extension of the protein-mediated perturbation on the bilayer. Based on a previous theoretical finding (Sperotto and Mouritsen, 1991), we first assumed that the perturbation induced by the protein on the surrounding lipids is of an exponential type. We have then verified this assumption later by analyzing the deviation of the functional form of the calculated d_L(r) from the one assumed. If the behavior of d_L(r) is exponential, the protein-induced perturbation can be expressed in terms of a typical coherence length, i.e., the decay length _P,

(16)

where is the mean hydrophobic thickness of the unperturbed pure lipid bilayer, and d_P is the protein hydrophobic length. The above equation expresses the fact that away from the protein surface, and at distances at least of the order of _P, the perturbed d_L(r) decays to the bulk value namely the value corresponding to that of the pure lipid system at the considered temperature, at least if no finite-size effects occur. In principle, by knowing d_L(r), d_P, and and by using Eq. 16, one can estimate _P. In our case we have determined the value of _P by best-fitting with Eq. 16 the values d_L(r) resulting from the simulations, where _P and are the fitting parameters. For the resulting value of the parameter obtained by the best-fitting, we have verified that this is equal, within statistical accuracy, to the value of the lipid bilayer hydrophobic thickness in the bulk, which we directly calculated from the simulations. Since the proteins can be subjected to tilt, the input parameter for d_P we used is not the actual hydrophobic length of the model-protein, but instead an effective length, The value is defined as the projection onto the normal of the bilayer plane of the protein hydrophobic length directly obtained from the simulations: where ^tilt is the tilt-angle (see Fig. 1 d). The degree of tilting of a protein with respect to the bilayer normal was computed by considering, for each bead-chain forming the protein, the vector that connects the position of the two hydrophobic beads bound to the protein hydrophilic beads (i.e., close to the lipid-water interface), one located in one monolayer of the bilayer, and the other in the opposite monolayer. The tilt-angle, ^tilt, is then defined as the average value, over all the chains forming the protein, of the angle between this vector and the bilayer normal.

In some cases, to facilitate the interpretation of the data obtained from the simulations, it was necessary to know the degree of order/disorder of the lipid chains in the vicinity of the protein, and eventually to compare it with that of the pure lipid bilayer, i.e., in the bulk, away from the protein-induced perturbation. Therefore, we have calculated the value of the lipid chain order parameter, S(r), which is defined as

(17)

with

(18)

where _S is the angle between the orientation of the vector r_ij = r_j – r_i (r_ij = |r_ij|) along two consecutive lipid chain beads, i,j, and the bilayer normal unit vector, . S has the value of 1 if r_ij is on average parallel to the bilayer normal, 0 if the orientation is random, and –0.5 if the bond is on average parallel to the bilayer plane. The value S(r) has been independently calculated for each of the two monolayers of the bilayer, as well as averaged over all the bonds of the lipid chains at distance r from the surface of the protein.

The details regarding the number of configurations used to estimate the statistical quantities were mentioned at the end of Statistical Ensemble and Surface Tension, above.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION MODEL AND SIMULATION METHOD RESULTS AND DISCUSSION CONCLUSION AND FUTURE... ACKNOWLEDGEMENTS REFERENCES

 
In this section, we present the results from the simulations of the pure-DMPC lipid bilayer model system and the bilayer with embedded proteins. We have focused on the low protein-concentration regime, where the correlation between different proteins can be neglected. Therefore we considered bilayers with just a single embedded protein. To investigate the dependence on mismatch and protein size of the extension of the lipid bilayer perturbation around an embedded protein, we first studied the behavior of the system at a constant temperature, well above the pure lipid bilayer main-transition, or melting, temperature. Because one of the ways to change the hydrophobic mismatch is by changing temperature, we then studied the temperature-dependence of _P and ^tilt in the temperature range above the melting temperature of the pure system, i.e., in the fluid phase. We did this for a number of lipid-protein model-systems.

For convenience, in the following sections, together with the given values of the reduced temperatures, T*, we have also added in brackets the corresponding approximated values in °C, estimated using Eq. 12. Because, as already stated, a direct mapping between the atomic level information and the model parameters is not possible, the values of the temperatures derived from the temperature-mapping previously described will also reflect this limitation. Therefore, in the following, the temperatures given in °C should be considered simply as general guidelines to which thermodynamic phase a system is in, at a given reduced temperature.

The results presented here refer to lipid bilayers composed of 900 lipids and 25 water beads per lipid, resulting in fully hydrated bilayers. The choice of 25 water beads per lipid was sufficient to ensure that the hydrophilic parts of the model protein would be fully hydrated even when the protein is not subjected to tilt. We have made calculations for smaller system sizes, and we have found that the chosen system size was sufficient to avoid finite-size effects, at least in the temperature range close or above the main-transition temperature of the pure bilayer system.

Pure lipid bilayer
Fig. 4 shows the phase behavior of the pure lipid bilayer hydrophobic thickness, as a function of reduced temperature T*. The system undergoes a main transition at a reduced melting temperature T*_m = 0.425, which is calculated from the inflection point of (T*). This value of the melting temperature corresponds to the main-transition temperature of DMPC, which is 24°C (Koynova and Caffrey, 1998). Above T*_m the lipid chains are in the melted state, hence the low value of and the system is in the so-called L, or fluid phase. The snapshot in Fig. 4 (bottom, right) shows a typical configuration of the system in the fluid phase. For sake of clarity, in this, as well as in all snapshots shown in this article, the water molecules are not shown. At very low temperatures the system is in the so-called L_ß' or gel phase, which is characterized by having ordered chains, hence the high value of In this phase the lipid chains are tilted with respect to the bilayer normal. A typical configuration at this temperature can be seen in the snapshot in Fig. 4 (bottom, left).

View larger version (29K): [in this window] [in a new window]   FIGURE 4  The pure lipid bilayer hydrophobic thickness, as a function of reduced temperature T* (top). The main-transition temperature of the pure system is at the reduced temperature T* = 0.425. Also shown are snapshots of typical configurations of the system simulated at reduced temperatures: T* < 0.35, corresponding to gel-phase, or Lß' (bottom, left); 0.425 > T* 0.35, corresponding to the ripplelike striated phase, which here we denote as Pß' phase (bottom, center); and T* > 0.425, corresponding to the L or fluid phase (bottom, right). The lipid headgroups are represented by black lines, the lipid tails by green lines, and the end-segments of the lipid tails are shown in yellow.

L

Using the scaling relation in Eq. 12, and the value of R_c = 6.46 Å for the conversion factor for the unit of length (see Length-, Time-, and Temperature-Scales, above), we can now compare the values of the bilayer hydrophobic thickness, and the area per lipid, obtained from our simulations with those referring to the fully hydrated DMPC bilayers, which are derived from experiments. Both sets of values are shown in Table 2. The values obtained from the simulations are in good quantitative agreement with the experimental data. Small deviations from the experimental values are only observed in the case of the area per lipid at high temperature (65°C), and in the case of the bilayer hydrophobic thickness in the gel phase (10°C).

View this table: [in this window] [in a new window]   TABLE 2  Values obtained from the simulations and from experiments of the pure bilayer hydrophobic thickness, and area per lipid, AL, of DMPC

Fig. 5 shows the calculated bilayer hydrophobic thickness profile, d_L(r), as a function of the distance r from the protein surface. The data refer to the different values of d_P, resulting in different values of hydrophobic mismatch d, ranging from d = –8 to 28 Å, and to the three protein sizes, which correspond to N_P = 4, 7, and 43. Because the probability of finding a lipid molecule in the lipid-shell closest to the protein (namely at the position r = nr with n = 1) is much lower than in the other lipid-shells (with n > 1), the data collected at a distance r = r from the protein surface have not been considered for the statistics. One can clearly see from the curves in Fig. 5 that the protein induces a perturbation of the lipid bilayer in its vicinity. The perturbation decays in a manner that depends on the hydrophobic mismatch, and on protein size. If the protein hydrophobic length is smaller than the unperturbed bilayer hydrophobic thickness (d_P < ), i.e., at negative mismatch d < 0 (open symbols in Fig. 5), to reduce exposure of the protein hydrophobic region to the water environment, the lipids around the protein shrink to match the protein hydrophobic surface. By choosing the peptide hydrophobic length to approximately match the value of the hydrophobic thickness of the unperturbed lipid bilayer, i.e., d 0, one can clearly see from Fig. 5 (crosses) that the perturbation induced by the protein on the surrounding lipids becomes negligible. Instead, when the chosen protein is such that d_P > i.e., at positive mismatch d > 0 (Fig. 5, solid symbols), the lipids in the vicinity of the protein, to match the protein hydrophobic surface, stretch and become more gel-like than the bulk lipids far away from the protein.

View larger version (21K): [in this window] [in a new window]   FIGURE 5  Lipid bilayer hydrophobic thickness profiles, dL(r), as a function of the distance r from the protein surface, for different hydrophobic mismatch, , and for the three protein sizes corresponding to (a) NP = 4, (b) NP = 7, and (c) NP = 43. The data refer to results from simulations made at the reduced temperature T* = 0.7, which is well above the main-transition temperature of the pure system, T*m = 0.425. The calculated value of the pure lipid bilayer hydrophobic thickness at T* = 0.7 is The symbols for d refer to the following estimated values of the protein hydrophobic length: the open circle to (negative d); the open triangle to (negative d); the plus symbol to (d 0); the solid triangle to (positive d); the solid circle to (positive d); and the solid diamonds to (positive d).

View larger version (25K): [in this window] [in a new window]   FIGURE 6  Calculated values of dL(r) (open circles) and fitted values using Eq. 16 (solid line) as a function of the distance r from the protein surface. The data refer to the three protein sizes NP = 4 (a,b), NP = 7 (c,d), and NP = 43 (e,f), and to two values of the protein hydrophobic length dP = 14 Å (d = –10 Å) and 41 Å (d = 17 Å). Also shown is the level value of the pure lipid bilayer thickness, (dashed line), the measured protein hydrophobic length dP (shaded area), and the effective protein hydrophobic length (open area), which is defined as the projection of dP onto the normal to the bilayer plane. The data refer to simulations at the reduced temperature T* = 0.7.

View this table: [in this window] [in a new window]   TABLE 3  Values of the decay length, P; the pure lipid bilayer hydrophobic thickness (both derived from fitting the thickness profiles dL(r) by using Eq. 16); the protein hydrophobic length, dP; and the effective protein hydrophobic length, (both calculated from the simulations) given for different values of hydrophobic mismatch, d, and for the three protein sizes corresponding to NP = 4, 7, and 43