Interactions of membranes with coarse-grain proteins: a comparison

We study the interactions between lipid bilayers and rigid transmembrane proteins by Monte Carlo simulations of generic coarse-grain models. Different popular protein models are considered and compared with each other, and key parameters such as the hydrophobicity and the hydrophobic mismatch are varied systematically. Furthermore, the properties of the membrane are manipulated by applying different tensions. The response of the membrane to the insertion of single proteins is found to be mostly generic and independent of the choice of the protein model. Likewise, the orientational distributions of single proteins depend mainly on the hydrophobic mismatch and the hydrophobicity of the proteins, and are otherwise similar for all protein models. Orientational distributions are generally found to be very broad, i.e. tilt angles fluctuate very much, in agreement with experimental findings. Weakly hydrophobic proteins respond to positive hydrophobic mismatch by tilting. Strongly hydrophobic (strongly bound) proteins distort the surrounding membrane and tend to remain upright. For proteins with intermediate hydrophobicity, the two mechanisms compete, and as a result, the tilt only sets in if the hydrophobic mismatch exceeds a threshold. Clusters of several strongly hydrophobic proteins with negative positive mismatch may nucleate raft-like structures in membranes. This effect is more pronounced for proteins with rough, structured surfaces.


Introduction
Biological membranes are a central component of all living creatures [1,2]. They are used to create compartments, to enclose substances and to control and regulate transport processes and signaling. All biomembranes have the same basic structure, i.e. a self-assembled phospholipid bilayer. Their individual properties and functionalities depend on the specific lipid content and, most importantly, on the associated set of membrane proteins. On average, proteins constitute 50% of the membrane mass, as in plasma membranes or exterior membranes. Relatively pure lipids with a low protein content of only about 18% are found around certain nerve fibers. In contrast, energy-transduction membranes have a very high content of proteins, typically around 75% [3].
Nowadays, our picture of biomembranes is based on the 'fluid mosaic model' of Singer and Nicolson [4,5], according to which the lipid membranes of biological cells are in a dynamic fluid state in which proteins are free to move around. The proteins themselves are basically found to be peripheral proteins that are only loosely attached to the membrane surface and integral membrane proteins that cannot be easily separated from the lipids. The latter type forms the major fraction of membrane proteins [6]. In their original 1972 model [4], Singer and Nicolson assumed that proteins and lipids diffuse freely and that the membrane structure on larger scales is largely homogeneous as a result. More recent experimental results suggest that the distribution of proteins in the membrane is in fact heterogeneous [5,7] and that some membrane proteins can be transiently trapped in certain membrane areas [8][9][10][11][12]. The mechanisms governing this lateral compartmentalization-whether it is driven by lipids or proteins-are still under debate. For example, the popular 'raft hypothesis' assumes that the driving mechanism is nanodomain formation in the underlying lipid matrix, which may be stabilized by lipid-protein interactions [13,14]. However, the nature of the rafts, and even the question of whether they really do exist in vivo, is still being discussed controversially [15].
The complexity of natural membranes makes it necessary to investigate simplified model bilayers which are composed of one or two different lipid species. Understanding the physical principles that govern the dynamic and structural behavior of these model membranes is the aim of a large number of both experimental as well as theoretical studies [16]. Unfortunately, structural perturbations or transformations of the lipid bilayer in the presence of proteins are among the most difficult processes to probe experimentally [17]. Thus, complementary theoretical approaches and computer simulations of membrane systems of well-defined compositions are necessary to elucidate the role of the lipid bilayer in processes such as protein aggregation and function. Depending on the length and time scales of interest, a variety of computational methods and models have been developed [16,18,19]. Microscopic studies in atomistic detail are restricted to relatively small system sizes, and the proteins are surrounded by only a relatively small ring of lipids in most simulations. To study the long-range influence of proteins on lipid bilayers, coarse-grained approaches are necessary [16,20]. These allow one to study generic aspects of lipid-protein interactions on larger scales, for example, the distortion of the lipid bilayers due to the presence of lipids [16,21]. This problem has also been considered by various theoretical approaches, ranging from molecular mean-field theories [22][23][24][25][26], to elastic continuum theories [27][28][29][30][31][32][33][34][35][36]. In all these studies, the integral proteins were represented by simplified objects, either as smooth cylinders (this representation is usually used in theoretical approaches) or as rigid objects made of beads (the most common approach in coarse-grained simulations [34,[37][38][39][40]). Such simplifications are based on the idea that generic aspects of membrane-protein interactions should not depend on the microscopic details of the protein structure. However, since even simplified objects, of course, have a microscopic structure (e.g. the structure of a smooth surface), it is not a priori clear which simulation results are generic and which are a specific property of the chosen model. So far, systematic comparisons of different protein models are lacking.
The present study attempts to close this gap. We compare the response of lipid bilayers to the insertion of different 'protein'-like inclusions, namely two variants of bead proteins and smooth spherocylinders. Key parameters such as the hydrophobic length of the inclusion and the strength of the hydrophobic interaction are varied systematically. Moreover, we also consider different fluid membrane states, i.e. tensionless membranes and membranes subject to a strong tension. Previous studies have shown that the internal structure of membranes changes under tension [41]. The monolayers are less well separated and one has a significant amount of interdigitation. Varying the tension, thus, allows us to assess the influence of the local lipid structure on the lipid-protein interactions.
Our study is based on a generic molecular membrane model [42,43], which has been shown to reproduce the main phases and phase transitions of single component phospholipid bilayers, including the high-temperature fluid phase L α , the low-temperature tilted gel phase L β and the intermediate ripple phase P β [43]. Moreover, the elastic properties of the membranes (such as the bending stiffness, area compressibility, etc) were found to be in semiquantitative agreement with those of bilayers made of dipalmitoylphosphatidylcholine (DPPC), which is one of the most abundant lipids in real biomembranes [44]. In previous work, we have investigated the membrane-mediated effective interactions between smooth cylindrical inclusions, either of infinite [41,44,45] or of finite length [45] in these model membranes. Here we will focus on the protein-lipid interactions and on the comparison of different protein models. We mainly consider single inclusions, but we will also discuss clusters of inclusions in the end.

Lipid bilayer
In this section, we will first briefly summarize our lipid model and then introduce the different protein models studied in this work. The lipid model was originally introduced in the context of Langmuir monolayers [46][47][48], and was shown to reproduce the generic phase behavior of fatty acid monolayers in the practically relevant region of the transition between the liquid expanded and the liquid condensed region [48]. Combined with a suitable, computationally cheap solvent model [49], it can also be used to study self-assembled lipid bilayers and their main phase transitions [42,43,50]. Lipids are represented by a linear chain of n tail beads (t) of diameter σ t , attached to one slightly larger head bead (h) with a diameter of σ h . They are immersed in solvent beads (s) of diameter σ s . Within lipid chains, beads are connected by a finite extensible nonlinear elastic (FENE) bond potential Here r 0 denotes the equilibrium length of the bond and r max the maximal deviation, i.e. subsequent beads within a chain cannot come closer than r 0 − r max or be pulled further apart than r 0 + r max . The angle θ between subsequent bonds in the lipid gives rise to a stiffness potential Beads that are not directly adjacent in the same chain interact through a truncated and shifted Lennard-Jones (LJ) potential. with The parameter σ i j = (σ i + σ j )/2 is the arithmetic mean of the diameters σ i of the interaction partners, and r c,i j = 1 σ i j for all partners (i j) except (tt) and (ss): r c,tt = 2 σ tt and r c,ss = 0. Hence tail beads attract one another, all other interactions are repulsive, and solvent beads do not interact at all with each other. The solvent outside the lipid bilayer behaves like an ideal gas fluid and has no internal structure. In the presence of solvent, the model lipid chains selfassemble spontaneously into bilayers. Specifically, we use the model parameters [42] n = 6 (i.e. the lipid chains comprise seven beads in total), σ h = 1.1 σ t , σ s = σ h , r 0 = 0.7 σ t , r max = 0.2 σ t , FENE = 100 /σ 2 t and BA = 4.7 , and we work at a nominal pressure of P = 2.0 /σ 3 t . As already mentioned in the introduction, the model reproduces the main phases of phospholipids, i.e. a high-temperature fluid L α phase at temperature k B T > k B T m ∼ 1.2 and a low-temperature tilted gel (L β ) with an intermediate modulated ripple (P β ) phase [43]. The energy and length scales can be mapped to SI units [44] by matching the bilayer thickness or, alternatively, the area per lipid and the temperature of the main transition to those of DPPC, giving 1 σ t ∼ 6 Å and 1 ∼ 0.36 × 10 −20 J.
The elastic material properties of the membranes in the fluid state were found to be comparable to those of DPPC membranes [44].

Protein models
All known integral membrane proteins are either composed of α-helices forming helical bundle proteins or β-strands, which produce β-barrel proteins. This has motivated generic models where integral proteins are represented by hydrophobic cylinders or cylinder assemblies.
The interaction strength between the hydrophobic core of the membrane and the hydrophobic part of the inclusion crucially influences the local perturbation of the bilayer. Since hydrophobic interaction does not arise from the binding of nonpolar molecules to each other, but from preventing polar solvent molecules from achieving optimal hydrogen binding, the strength of the hydrophobic interaction depends on the relative polarity of both the solute and the solvent [51]. Experimentally, the hydrophobicity of proteins can therefore be tuned by changing amino acid residues with different hydrophobicities of the protein [52]. As an example, alanine is less hydrophobic than leucine [53]. Alternatively, changing the pH of the solvent and thereby changing the ion dissociation on side chains will also affect the hydrophobic interaction of the lipid bilayer and the proteins [54].
A second critical parameter is the hydrophobic length of the protein, i.e. the length of the hydrophobic section on the protein, compared to the hydrophobic thickness of the membrane. In the case of 'hydrophobic matching', the hydrophobic length of any transmembrane domain matches the hydrophobic thickness of the bilayer [55]. If the thickness of the hydrophobic core of the unperturbed bilayer is larger than the hydrophobic length of the protein, one has 'negative mismatch', and if the hydrophobic length of the protein surmounts the bilayer thickness, one has 'positive mismatch'. Hydrophobic mismatch can have various consequences for proteins, such as tilt or conformational changes of transmembrane parts, lateral oligomerization or even failure of membrane insertion [32]. Furthermore, it may also influence the properties of the lipids [56,57]. The lipid chain order may be changing, the phase transition temperature can be modified, the formation of microdomains is possible or even nonlamellar structures may be induced.
Here we compare three different models for simple transmembrane proteins. All have a hydrophobic middle section, i.e. a section which attracts tail beads. The length of this section and the strength of the hydrophobic interactions can be tuned. A picture of the three models is shown in figure 1.
The first type of model protein is a smooth spherocylinder with a hydrophobic length L, capped at both ends by effectively repulsive ('hydrophilic') hemispheres [45] (figure 1, left). It is parameterized by a line of length L. The interactions between proteins and lipids or solvent beads have a repulsive contribution where r denotes the shortest distance between the protein line and the center of the bead, σ is given by σ = (σ t + σ i )/2 for interactions with beads of type i (i = h, t and s for the head, tail and solvent beads, respectively), σ 0 = σ t and V LJ has been defined above (equation (4)). In addition, tail beads are attracted to the straight inner section of the protein by an attractive potential that depends on the projection of the distance d between the centers of the tail bead and the protein onto the protein axis. The total potential for tail beads reads with the attractive Lennard-Jones contribution and a weight function W P (d), which is unity on a stretch of length 2 l = L − 2σ t and crosses over smoothly to zero over a distance of approximately σ t at both sides. Specifically, we use The hydrophobicity of the protein can be tuned by varying the parameter pt . In systems containing several proteins, the direct protein-protein interactions are purely repulsive and have the form of equation (5) with σ = σ t , σ 0 = 2σ t and r the minimum distance between the two protein lines. In our own previous work [41,44,45], we used a restricted variant of this model, where the cylinders had infinite length (albeit with a finite hydrophobic section L) and fixed orientation in the direction normal to the membrane. This corresponds to the situation most commonly considered in theories of protein-induced bilayer distortions. Here we focus on more realistic simulation models and study cylinders with full orientational freedom. Representing a complex structure like a protein by a compact cylindrical object might seem a rather crude approach. However, this can be justified by the fact that, e.g., α-helices are packed with very little free space within the helices [58] and are therefore fairly smooth on the scale of ∼10 Å. There are no large cavities into which chains or even whole molecules would fit [59].
The other two protein models have surfaces of varying roughness. The second type of model protein, which we will call 'rough bead protein' in the following, consists of hydrophobic (pT ) and hydrophilic (pH ) beads having the same diameter as the tail beads, i.e. σ pT = σ pH = 1 σ t . This way of modeling proteins is common in coarse-grained simulation studies of membrane-protein interactions [34,35,[37][38][39][40]. Specifically, the proteins are constructed as rigid stacks of discs separated by a distance d pT = σ t . Each disc consists of an outer ring of N p beads with the nearest-neighbor distance d pT , which is filled with further beads if necessary. The diameter D p of the protein is therefore given by One disc of hydrophilic beads is added at each end of the protein. Hence, the total length of a bead protein is L = L pT + 2L pH , where L pT sets the hydrophobic length (see figure 1, middle). The hydrophobicity of the protein can be tuned with the interaction parameter pT t governing the strength of the attractive interactions between hydrophobic protein and lipid beads. Beads within one protein are rigidly ordered and do not have any additional degrees of freedom. In comparable dissipative particle dynamics (DPD) simulations where protein beads were connected by springs, no appreciable internal deformation of the proteins was observed, except for a slight bending of very slim proteins [38]. Compared to the smooth spherocylinder, the surface structure of the rough bead proteins is rather corrugated. Our third model protein, which we denote by 'smooth bead protein' is constructed as an intermediate model with reduced corrugation. This is achieved by doubling the number of beads and reducing the minimum separation to d pT = 0.5 σ t . To obtain hydrophobic interaction strengths that are comparable to the rough bead protein model, the interaction parameter pT t has to be rescaled. The rescaled parameters will be marked as˜ in the following. As a first estimate, simple geometrical considerations suggest that a rescaled energy of˜ pT t = 0.25 should lead to similar behavior of the smooth bead protein as pT t = 1.0 for the rough bead protein.
In this paper, we will discuss proteins of diameter D p ∼ 3σ t , corresponding to the diameter of a β-helix such as, e.g., gramicidine [44]. Thus our rough bead proteins contain N p = 6 beads per disc, and the smooth bead proteins contain N p = 12 beads per disc (cf figure 1). We have also studied proteins with smaller diameter (down to D p ∼ 1σ t , corresponding to an α-helix) and larger diameter (up to D p ∼ 5σ t ). The results were qualitatively similar and will not be shown here. Details can be found in the theses [60,61].

Simulation method
The systems described above were studied by Monte Carlo simulations at constant pressure, constant temperature T = 1.3 /k B , constant pressure P = 2.0 /σ 3 t and constant surface tension (with = 0 unless stated otherwise) in a simulation box of variable size and shape. Following [62], we impose tension via an additional energy term − A in the Hamiltonian of the system, where A is the projected area of the bilayer onto the x y-plane. The noninteracting solvent particles, which probe the free volume and force the lipids to self-assemble, are not affected by this additional energy contribution. They ensure that the normal pressure P N is kept fixed at the required value. Thus, we are performing Monte Carlo simulations in the N P N T ensemble with the effective Hamiltonian where U is the interaction energy, V the volume of the simulation box, V 0 = σ 3 t our reference volume and N the total number of beads (cf [42]). We use three main types of Monte Carlo moves, namely (i) local moves of lipid or solvent beads, (ii) global moves which change the size or shape of the simulation box and involve rescaling of all particle coordinates [42] and (iii) protein rotation or translation moves [45]. The moves are proposed randomly and accepted according to a Metropolis criterion. In each Monte Carlo sweep (MCS), every bead is moved once on average and the protein is moved and rotated once. The computationally expensive global rescaling moves were attempted every 50th MCS.
The system sizes ranged from ∼780 lipids and ∼12 300 solvent beads for the thin proteins to ∼1640 lipids and ∼24 600 solvent beads for the thickest proteins. Typical run lengths were of the order of several million MCS with equilibration times of up to one million MCS.

Results: single proteins
In this section, we compare the interactions between lipid membranes and single proteins for our three different protein models. Since the direct interactions between the proteins and the lipid molecules are rather different in the three models, the direct quantitative comparison is nontrivial. To set the stage from a thermodynamic point of view, we first consider the free energies of protein insertion for the different models. Then we discuss the influence of the proteins on the surrounding lipid bilayer, and finally, the orientational distributions of the proteins.

Binding free energies
The effective binding energy of the proteins to the membrane can be determined from the Gibbs free energies of insertion. In our context, the quantity of interest is the difference G eff = G − G s between the Gibbs free energy of inserting a protein in a membrane and the Gibbs free energy of inserting a protein in a pure solvent. To determine G and G s , we use a variant of the Widom insertion method [63] and gradually insert the protein by modifying its interaction potentials with a parameter λ. At λ = 0.0 the interaction must vanish completely and at λ = 1.0 it reaches its full interaction strength. The difference in Gibbs energy G can then be calculated by thermodynamic integration, The derivative ∂ H/∂λ can be calculated analytically and its value is recorded similarly to the other observables during the simulation. For the smooth cylinders, we replaced the total interaction energy U ip between the protein and lipids or solvent by a rescaled energỹ At low values of λ, this potential is not sufficient to bind the proteins to the membrane, which results in sampling problems. Therefore, we used a restricted model where the protein cylinders had infinite length (but finite hydrophobic portion L) and thus stay in the membrane by construction [44]. Here, 'infinite length' means that the cylinder extends through the whole simulation box and is connected by the periodic boundary conditions, which implies that it cannot tilt.
In the case of the bead proteins, the interaction potentials are all of Lennard-Jones type and we can follow the approach of Beutler et al [64]: the modified interaction energy between the beads of the gradually inserted protein k and the (unaltered) beads of the lipids and the solvent i at distance r ki reads where ki,λ,cutoff is set such that the potential is continuous at the cutoff radius r ki,c . In the case of purely repulsive soft-core potentials, the cutoff radius is also shifted to r ki,λ,min = 6 1 − α(1 − λ)σ ki (14) in order to account for the shift of the local minimum of the modified Lennard-Jones potential. The values of n and α LJ can be tuned such that the proteins remain in the membrane for all values of λ. Good results were obtained with n = 1 and α LJ = 0.40. Figure 2 shows the results on the binding energies for different proteins with hydrophobic length L = 6σ t . These proteins are 'hydrophobically matched'; therefore deformation of the bilayer is minimal (see the next subsection) and the binding energy results mainly from the competition between the interaction energy and the entropy loss associated with conformational changes in the lipid bilayer. The comparison of binding energies enables one to relate the interaction parameters in the different models with each other. Figure 2 shows that smooth bead proteins with hydrophobic interaction value˜ pT t have roughly the same binding energies as rough bead proteins with hydrophobic interaction value pT t = 3˜ pT t and cylinders with hydrophobic interaction value pt = 10˜ pT t . This relation will be further supported below when inspecting quantities such as the bilayer thickness or the director fields around cylindrical proteins and bead proteins, respectively.
The influence of the binding energies on the hydrophobic length of the protein is examined in figure 3 for the case of rough bead proteins. For a small hydrophobic strength pT t = 0.5, the binding energy is positive, G eff > 0; thus the protein does not bind at thermodynamic equilibrium. Nevertheless, it was possible to prepare metastable states where such weakly hydrophobic proteins span through a membrane, and these states mostly remained metastable during the whole simulation run. Only proteins with a small hydrophobic length L pT = 4σ t were occasionally expelled out of the membrane. The transmembrane state is stabilized by a high kinetic barrier, created by the fact that the hydrophilic caps have to traverse the hydrophobic core of the membrane during an expulsion process. A related observation was made by Illya and Deserno [65] in a recent study of peptide-induced pore formation. They report that at a certain peptide-lipid attraction, proteins initially placed above the bilayer bound to the upper monolayer, but did not insert. If the same peptide was initially placed inside the bilayer with its long axis parallel to the bilayer normal, it remains in the bilayer. The 'binding energies' in the nonbinding regime do not depend on the hydrophobic length L pT of the protein. As soon as the protein binds, however ( G eff < 0), the binding energy decreases with increasing L pT . Figure 3 also shows binding energies in membranes under the tension = 2 /σ 2 t . The trends are similar and even the values are comparable.

Bilayer distortion close to proteins
Next we investigate the effect of the protein on the surrounding lipid bilayer. One particularly pronounced phenomenon is the distortion of the membrane thickness in the vicinity of the proteins. This is shown for the different protein models in figure 4. Close to weakly hydrophobic proteins with positive binding energies ( G eff > 0), the membrane thickness is reduced, compared to the bulk: the protein effectively repels the lipids. Strongly hydrophobic proteins with G eff 0 locally compress or expand the membrane depending on the sign of hydrophobic mismatch: the membrane thickness adjusts to the hydrophobic length of the protein in the vicinity of the protein, and relaxes at larger distances. The thickness profiles were determined as the mean distance between opposing head beads at a given distance r from the proteins. Their shapes can be fitted nicely with an elastic theory originally developed by Safran and coworkers [27-29, 35, 44], shown as solid lines in figure 4. This theory treats the bilayer as a system of coupled elastic monolayers, each having a mean thickness t 0 , a bending rigidity k c /2, an area compressibility k A /2 and a spontaneous curvature c 0 . Furthermore, an additional parameter ζ enters the theory, which is related to the derivative of the spontaneous curvature with respect to the lipid area. The parameters k c , k A and ζ have been determined independently for our model bilayers from the fluctuation spectrum, both for tensionless membranes [44] and for membranes under tension [41]. They are given in table 1. The general form of a radially symmetric monolayer thickness profile (r ) = t (r ) − t 0 reads as [29] (r ) = A 1 [41].
where J 0 (x) and Y 0 (x) are the zeroth-order Bessel functions of the first and second kind, and the parameters α ± are complex numbers in stable membranes, α ± = α r ± iα i (with real and positive α r , α i ). Since the profiles must not grow exponentially at infinity, the coefficients The remaining complex coefficient A depends on the boundary conditions at the surface r = R of the protein. For a given surface distortion (R) =: t R and surface curvature ∇ 2 r | R =: t R at the radius r = R (with ∇ 2 r = (1/r )∂ r r ∂ r ), it reads .
The parameters t R and t R are the fit parameters in the theoretical curves of figure 4; they characterize the surface of the protein. More precisely, the theory predicts [35,44] t R + 2ζ where t R = ∂ r | R , k G is the Gaussian curvature andc 0 subsumes the effect of the spontaneous curvature c 0 as well as possible free-energy contributions from local distortions of other quantities that couple to the membrane thickness [44]. Unfortunately, neither c 0 nor k G can be determined from the fluctuation spectrum. They can be estimated from the moments of the pressure profiles, but the estimate is not very reliable, especially in the case of k G . We will therefore use the parameters t R andC as defined in equation (18) to characterize the effect of the protein surface on the thickness profiles. In figure 5, we plot the parameterC against the bilayer distortion at the surface. For weakly hydrophobic proteins (leftmost curves), there is no clear dependence. For strongly hydrophobic proteins, however, the curves approach one common line for all protein models. We conclude that the main characteristics of the bilayer thickness profile in the vicinity of strongly bound proteins do not depend on the particularities of the protein model. It is worth noting that the amplitude of the bilayer deformations for the proteins studied here is comparable to that of proteins with fixed upright orientation (as determined in [44]). This is related to the fact that strongly hydrophobic proteins exhibit very little tilt (cf section 3.3). Another bilayer property that is influenced by a transmembrane protein is the orientation of the lipids. Even though lipid tails are fluctuating and permanently changing their conformation within the hydrophobic core of the bilayer, they have average orientations [66], which may be shifted in the presence of a protein. The concept of lipid director has been applied in theories for the lipid-mediated interaction free energies between hydrophobic surfaces [26] and membranemediated interaction between two cylindrical inclusions in a symmetric lipid bilayer [32]. Figure 6 shows radial profiles of the average lipid tilt direction for the spherocylinder protein model and for the rough bead protein model. The profiles for positive and negative mismatch are qualitatively different. Close to positively mismatched, strongly hydrophobic proteins, the tilt profile is nonmonotonic: as one approaches the protein, the lipids first tilt towards the protein, then they become straight and may even slightly tilt away from the protein.
Lipids near negatively mismatched tilt away from the protein at all distances, but the tilt profiles also become nonmonotonic for strongly hydrophobic proteins, such that the tilt exhibits a maximum at a distance of around r ∼ 4σ t . This complex behavior is found for both types of protein models. The tilt profiles for the two protein models are almost identical (figure 6).

Orientational distribution of proteins
Next we examine the orientation distributions of our model proteins in membranes. The orientation of proteins is believed to have a significant influence on their functionality, e.g. in the context of pore formation [67]. Recent coarse-grained simulations have suggested that the cross-angle distributions of packed helix complexes are mostly determined by the tilt angle of individual helices [68]. One important driving force leading to tilt is hydrophobic mismatch [ 38,69,70]. Proteins tilt in order to alleviate the free-energy costs associated with membrane deformations. This is predicted by theoretical considerations [71] as well as molecular dynamics simulations [38,72,73].
Experimental tilt measurements, e.g. by nuclear magnetic resonance (NMR) techniques, have in some cases supported this view [69,74]; in other cases the reported tilt angles were surprisingly small compared to theoretical expectations [75][76][77]. This has been explained by large orientation fluctuations, which complicate the interpretation of NMR data, especially if peptides are highly mobile [73,78,79]. Higher tilt angles are obtained if such fluctuations are taken into account in the analysis [80,81]. Furthermore, tilt can be influenced by the anchoring residues flanking the hydrophobic transmembrane domains, which have their own preferred orientation at the hydrophobic/hydrophilic interface [70,82,83] and might prevent tilting through a variety of mechanisms [76].
Our simulations reveal yet another factor that controls tilt in proteins: the hydrophobicity of the protein. This is demonstrated in figure 7, which shows average tilt angles of different protein models at different hydrophobicities as a function of the relative hydrophobic mismatch (L − 2t 0 )/2t 0 . (We recall that L is the hydrophobic length of the proteins and 2t 0 the hydrophobic thickness of the membrane.) At negative hydrophobic mismatch or for hydrophobically matched proteins, the average tilt angle takes values around α ∼ 10 0 , which is in the same range as experimental values [81]. For positively mismatched proteins, the behavior depends markedly on the hydrophobicity of the proteins. In figure 7, the hydrophobicity of the 'weakly hydrophobic' proteins is so small that the binding energy G eff > 0, i.e. the bound state is thermodynamically metastable only. The tilt angle of such metastably bound proteins increases with increasing positive hydrophobic mismatch as expected. In contrast, strongly hydrophobic proteins with large binding energies, G eff 0, exhibit average tilt angles which are almost independent of the hydrophobic mismatch and tend to be smaller than those of the corresponding weakly hydrophobic proteins. This result is unexpected and seemingly at Cylinder, weakly hydrophobic (ε pt =1ε) Cylinder, strongly hydrophobic (ε pt =6 ε) Beads, weakly hydrophobic (ε pTt =0.5 ε) Beads, strongly hydrophobic (ε pTt =2 ε) variance with experimental findings ofÖzdirekcan et al [76], who reported a slight increase of tilt with hydrophobicity. However, it should be noted that the tilt angles reported in this study were generally very small, between 5 • and 10 • , and the results might be affected by the abovementioned difficulties of analyzing NMR data for highly mobile peptides.
To analyze this unexpected phenomenon in more detail, we show in figures 8 and 9 a selection of the corresponding orientational distributions. In general, the distributions are very broad, which is consistent with the experimental picture that orientations fluctuate strongly [78,79,81]. If one increases the hydrophobic length L, the orientational distribution of weakly hydrophobic proteins initially broadens, i.e. more proteins have higher tilt angles. This results in an increased average tilt angle. However, the maximum of the distribution is still found at tilt angle zero. The proteins fluctuate strongly, but their mean position is straight ( figure 8). If one increases the positive hydrophobic mismatch even further by applying tension, thus reducing the membrane thickness, a second effect comes into play: the orientation distribution not only broadens further, but also develops a peak at nonzero angle α. We conclude that weakly hydrophobic proteins adjust to the membrane thickness by tilting as expected ( figure 9).
In the case of strongly hydrophobic proteins, the situation is different. The orientation distribution does not broaden with increasing L; instead it sharpens around tilt angle zero, such that longer proteins are on average less tilted than shorter proteins. This behavior is found both for smooth spherocylinder proteins and for bead proteins; hence it seems to be generic. It can be associated with the deformations that strongly hydrophobic proteins induce in the surrounding lipid membrane: the stretched lipids bound to the protein surface stabilize the upright orientation. They order the surrounding lipids, leading to the formation of a dynamic complex consisting of the protein and a lipid shell, which is preferably oriented normal to the  membrane. A similar effect has been reported for certain model WALP peptides with flanking tryptophane residues [84]. In that case, the anchoring residues were held responsible for the lipid stretching. Apparently, proteins which are capable of inducing significant membrane deformations in their vicinity by whatever mechanism-either due to their hydrophobicity or due to strongly anchoring residues-will tend to form condensed protein-lipid complexes instead of tilting. Real proteins are in a thermodynamically bound state, but the binding energies are not necessarily very large. Figure 10 shows the behavior of the average tilt for such more realistic proteins as a function of relative hydrophobic mismatch. In this case, the two mechanisms of adjusting to the mismatch compete: as long as the hydrophobic mismatch is small, the protein responds by (slightly) deforming the membrane (cf figure 4) and the surrounding lipids stabilize an upright orientation. If the hydrophobic mismatch becomes large, the protein tilts. Due to the competition, the onset of tilt is not identical with the point where the hydrophobic mismatch turns positive. As a result, the protein behaves as if its hydrophobic length were effectively reduced (figure 10).
To conclude, our simulations suggest that the hydrophobicity of proteins or, more generally, their ability to induce strong membrane deformations plays an important role in determining the tilt. Tilting competes with the formation of dynamic complexes consisting of proteins and their surrounding lipid shells. This second mechanism will be more important if the proteins are more strongly bound in the membrane. We expect that it will also gain importance with increasing protein radius. Most of the systematic experimental studies cited above were based on transmembrane proteins with α-helical structure. In our simulations, we have considered proteins with the radius of β-helices (comparable to gramicidine), which is about three times as thick. Coarse-grained simulations indicate that thicker proteins tilt less than thin proteins [38]. Thus experimental investigations on the interplay of membrane deformation and tilting are presumably more promising if one uses experimental model proteins based on β-barrels. Figure 11. Evolution of the COM angles φ i and the COM distances l i during 3 × 10 5 MCS for systems containing three moderately hydrophobic rough bead proteins with negative mismatch (L = 4σ t , left), no mismatch (L = 6σ t , middle) and positive mismatch (L = 8σ t , right). The hydrophobicity is pT t = 1 . Typical configurational snapshots are shown below the graphs.

Outlook: protein clusters
The next step is to consider membranes that contain several proteins. We have studied the distortion profiles of membranes containing two proteins and, as before, we found no noticeable differences for the different protein models (data not shown, see also the thesis [61]). Membrane-mediated interactions between two proteins have been studied for the spherocylinder model by us [44,45] and for bead protein models by other authors [85,86], and the general features are similar. However, it is well known that membrane-mediated protein-protein interactions are not pairwise additive [87]. Many-body effects are important even at low densities. Therefore, we will conclude with a brief discussion of many-body effects. We have studied the time evolution of membranes containing three proteins, which were initially set up on an equilateral triangle with mutual distance r ∼ 6.5σ t . After equilibration of the system, it was monitored over 3 × 10 6 MCS. Quantities used to characterize the orientation of the three proteins with respect to each other are their distance l i (i = 1, 2, 3) from their common center of mass (COM) and the angles φ i between the two vectors l i and l i+1 pointing from the COM to the center of two adjacent proteins (with l 4 ≡ l 1 ). In the following they will be called COM distances and COM angles, respectively. The angles usually do not give exactly 360 • when added up, since we are not measuring the projection of the angles onto the x y-plane, but the full angle in space. Figure 11 shows the corresponding time evolution for rough bead proteins with moderate hydrophobicity. Here we observe a new effect of hydrophobic mismatch which is not yet present in systems containing only one or two proteins. Negatively mismatched proteins may nucleate an ordered lipid state in their direct environment which in turn pins the proteins at their positions, such that they effectively freeze. A similar, albeit weaker, effect is observed for positively mismatched proteins. In contrast, hydrophobically matched rough bead proteins remain mobile and diffuse around slowly during the simulation run.
The onset of a similar effect can be observed for smooth spherocylinder proteins, but only at much higher hydrophobic strengths. For example, if one raises the hydrophobic strength by a factor of 2, clusters of negatively mismatched proteins begin to freeze, but the corresponding positively mismatched cylinders remain mobile (figure 12).

Summary
In summary, we have studied the influence of the protein on the membrane structure as well as the influence of the membrane on the protein orientations for different protein models. The protein models have in common that they represent rigid cylindrical structures with no internal degrees of freedom. They differ in the degree of roughness of the surface.
For single proteins, we found that the distortion of the lipid bilayer as well as the orientational distribution of the protein are mainly determined by a few generic key parameters such as the hydrophobic length or the strength of the hydrophobic interaction. Details of the protein structure do not matter. Regarding all the quantities considered here, the results for the different protein models with different surface corrugation could be related to each other at an almost quantitative level. The bilayer thickness profiles can be fitted reasonably well for all parameter values with an elastic theory, and the relation between the elastic parameters did not depend on the protein model. The lipid tilt profiles around the protein as well as the histograms of protein tilt angles were almost identical for the different protein types.
Regarding the protein tilt distribution, we observed an unexpected qualitative difference between weakly and strongly hydrophobic proteins: whereas weakly hydrophobic proteins exhibit strong tilt fluctuations in the membrane, strongly hydrophobic proteins form a complex with the surrounding lipids which keeps them in an upright position. Weakly hydrophobic proteins respond to hydrophobic mismatch by tilting. Strongly hydrophobic proteins remain untilted and distort the membrane. In the biologically most relevant case of moderately hydrophobic proteins that are weakly bound to the membrane (i.e. the hydrophobic strength is just about large enough that the insertion free energy is negative), the two mechanisms compete. Upon increasing the hydrophobic mismatch, the proteins first slightly distort the membrane and then tilt. As a result, their apparent hydrophobic length is reduced. This was observed in all protein models. Overall, the differences between the different protein models were negligible.
We conclude that the microscopic geometric structure has very little influence on the interaction of lipid bilayers with single proteins. The particular choice of the protein model is not very critical as long as one is mainly interested in generic features. However, this changes when looking at systems containing several proteins. In systems containing three hydrophobically mismatched rough bead proteins, we have observed trimerization through the nucleation of an ordered membrane domain. For smooth spherocylinder proteins, this effect was much weaker and only observed for strongly hydrophobic, negatively mismatched proteins. Thus the specific choice of the protein model becomes important if one wants to study complexes of lipids and several proteins. Translating this finding to experimental situations, we would expect the behavior and structure of lipid-protein complexes containing several proteins to be sensitive to details of the protein surface structure even in systems where the characteristics of single protein constituents (i.e. simple α-helices or β-barrels) are mostly determined by generic factors such as hydrophobic mismatch.
The observation that hydrophobically mismatched proteins can nucleate ordered domains with reduced mobility might be interesting in the context of the current raft discussion. We emphasize that these ordered structures have no thermodynamically stable membrane phase counterpart in the absence of proteins. As mentioned in the introduction, the role of proteins in raft formation (if they exist) is not yet clear. Our results indicate that multibody hydrophobic mismatch interactions might provide a possible mechanism that stabilizes rafts. This mechanism should become even more efficient if an ordered lipid phase is close by in parameter space. Similar mechanisms could also stabilize 'liquid ordered' domains with enhanced thickness in mixed membranes even in parameter regions where the liquid ordered state does not represent a distinct thermodynamically stable phase in the pure lipid membrane [88]. Future studies of model membranes containing many proteins in varying lipid environments should shed light on these issues.