Membrane-mediated protein-protein interaction: A Monte Carlo study

We investigate membrane-mediated interactions between transmembrane proteins using coarse-grained models. We compare the effective potential of mean force (PMF) between two proteins, which are always aligned parallel to the z-axis of the simulation box, with those PMFs obtained for proteins with fluctuating orientations. The PMFs are dominated by an oscillatory packing-driven contribution and a smooth attractive hydrophobic mismatch contribution, which vanishes if the hydrophobic length of the protein matches the thickness of the membrane. If protein orientations are allowed to fluctuate, the oscillations are greatly reduced compared to proteins with fixed orientation. Furthermore, the hydrophobic mismatch interaction has a smaller range. Finally, we compare the two-dimensional thickness profiles around two proteins with the predictions from the elastic theory of two coupled monolayers, and find them to be in very good agreement.


I. INTRODUCTION
Biological membranes are essential components of all living organisms. Even though their basic building block is a fluid lipid bilayer, their properties and functionality crucially depend on the membrane proteins. Proteins function as catalysts, transport and store other molecules, provide mechanical support and immune protection, generate movement, transmit nerve impulses, and control growth and differentiation [1].
Unfortunately, the components of real biomembranes are too diverse and complex to obtain detailed and unambiguous information about interactions of transmembrane segments with a lipid bilayer [2]. Moreover, structural perturbations or transformations of the lipid bilayer are among the most difficult processes to probe experimentally [3]. 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 like protein aggregation and function. In this context, the use of coarse-grained models has become popular [4][5][6][7]. Even with today's computing resources, atomistic modeling of multicomponent lipid bilayers on length scales of several nanometer still remains a challenge. Moreover, coarsegrained model simulations give insight into generic properties and mechanisms in lipid membranes, which are difficult to extract from fully atomistic simulations of specific membranes.
In this paper we focus on lipid-mediated interactions between proteins or, more generally, membrane inclusions. We present simulational and theoretical results on fluid membranes containing model transmembrane proteins, i.e. proteins that span through the membrane. Many theoretical models have been developed to get a better understanding of the membrane-mediated inter- * Electronic address: friederike.schmid@uni-mainz.de actions between such proteins [8][9][10][11][12][13][14][15]. In these models, the inclusions were modeled as straight cylinders, which are aligned with the bilayer normal. Only recently, computational studies with coarse-grained models have been performed on the same problem, where both proteins in single-component bilayers [16][17][18] and the effect of cholesterol on protein-lipid and lipid-mediated protein-protein interactions have been studied [19]. From these studies, the following general picture has emerged: The lipidmediated protein interactions can be divided in shortrange and long-range contributions. The short-range contributions depend on the local structure of the lipid layer and reflect packing and layering effects. The longrange interactions result from the elastic distortion of the membrane due to the insertion of the proteins and can be tuned, by tuning the length of the hydrophobic section of the protein ("hydrophobic mismatch" interaction). If the hydrophobic length of the proteins matches the thickness of the membrane, they will vanish, otherwise they tend to be attractive. In the literature, a number of other mechanisms that would induce long-range lipid-mediated protein interactions have been discussed [20][21][22], but for straight cylindrical proteins, the hydrophobic mismatch interaction seems be dominant.
The work presented here is based on a coarse-grained lipid model [23,24], which reproduces the main phases [25][26][27] and elastic properties [18,27] of DPPC bilayers. In previous papers, this model was used to study the lipid-mediated interactions between infinitely strong straight cylinders with fixed orientation along the bilayer normal in stressfree [18] and stressed membranes [27]. The rigid restriction was motivated by the huge amount of related theoretical work (see above). In reality, however, proteins may tilt in the membrane, their orientation may fluctuate, and this affects the lipid-mediated protein interactions.
The present paper focuses on this effect. We study the lipid-mediated protein interactions, using two variants of a coarse-grained protein model. The first variant is the infinitely long straight cylinder studied earlier [18,27], and the second variant consists of a cylinder of finite length with full freedom to tilt. This allows to assess in detail the effect of orientation fluctuations on the protein-protein interactions. In addition, we characterize the profiles of membranes containing two proteins and compare them with the theoretical prediction of an elastic theory.
Our paper is organized as follows: In the next section, we introduce the simulation models and method and briefly comment on the elastic theory with which the simulation data are compared. The simulation results are presented in Section III. We summarize and conclude in Section IV.

II. MODELS AND METHODS
In the following a brief summary of the lipid model and the protein models used in this work will be given. In our study, we vary both the hydrophobic mismatch of the proteins (i.e., the length of the hydrophobic section), and the hydrophobic strength. 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 [28]. Experimentally the hydrophobicity of proteins can therefore be tuned by changing amino acid residues with different hydrophobicity of the protein [2]. Alanine, e.g., is less hydrophobic than leucine [29]. Alternatively, changing the pH of the solvent and thereby making priorly uncharged side-chains charged, will also affect the hydrophobic interaction of the lipid bilayer and the proteins [30].

A. Lipid Model
Lipid molecules are represented by chains consisting of one head bead and six tail beads, and there are additional solvent beads [24]. Within the lipid chain, adjacent beads at a distance r interact via a finite extensible non-linear elastic (FENE) potential with the spring constant FENE , the equilibrium bond length r 0 , and the cutoff ∆r max . The angles θ between subsequent bonds in the lipid are subject to a stiffness potential with the stiffness parameter BA . Beads of type i and j which are not direct next neighbors in a chain interact via a truncated and shifted Lennard-Jones potential, The offset V c,ij is chosen such that V LJ (r/σ ij ) is continuous at the cutoff r c,ij . The parameter σ ij = (σ i + σ j )/2 is the arithmetic mean of the diameters σ i of the interaction partners, and r c,ij = 1 σ ij for all partners (ij) 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. We use the model parameters [24] σ h = 1.1 σ t , r 0 = 0.7 σ t , ∆r max = 0.2 σ t , FENE = 100 /σ 2 t , and BA = 4.7 . At the pressure P = 2.0 /σ 3 t , the model reproduces the main phases of phospholipids, i.e., a hightemperature 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 [25]. The energy and length scales can be mapped to SIunits [18] 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 properties of the membranes in the fluid state were then also found to be comparable to those of DPPC membranes [18].

Cylinder of infinite length
The first type of model protein considered here is the infinite straight cylinder originally introduced by West et al. [18]. (see Fig. 1, left), i.e., the protein cannot tilt. The axis of the cylinder is always aligned parallel to the z axis of the simulation box. Since this model protein spans the whole simulation box in the z direction and due to the periodic boundary conditions, the term "infinitely long cylinder" is justified. The interaction of this simple model proteins and the lipid or solvent beads has repulsive contributions, which are described by a radially shifted and truncated Lennard-Jones potential where r = x 2 + y 2 denotes the distance of the interaction partners in the xy plane, σ is given by σ = (σ t +σ i )/2 for interactions with beads of type i (i = h, t, and s for head, tail, and solvent beads, respectively), σ 0 = σ t , and V LJ has been defined above (Eq. (3)). The direct proteinprotein interactions have the same potential (Eq. (4)) with σ = σ t and σ 0 = 2σ t . In SI-units, the cylinders thus have a diameter of 2 nm, which roughly corresponds to the diameter of the 8-stranded β-barrel OmpA [31]. In addition, protein cylinders attract tail beads on a hydrophobic section of length L. This is described by an additional attractive potential that depends on the z distance between the tail bead and the protein center. The total potential reads with the attractive Lennard-Jones contribution and a weight function W P (z), which is unity on a stretch of length 2l = L − 2σ t and crosses smoothly over to zero over a distance of approximately σ t at both sides. Specifically, we use The hydrophobicity of the protein is tuned by the parameter pt .

Cylinder of finite length
The second type of protein is also a straight cylinder, where the hydrophobic section of length L is modeled analogously to the infinite cylinder. But now the cylinder is allowed to tilt away from the z direction of the simulation box (Fig. 1, right). The protein has finite length L and is capped at both ends by effectively repulsive ("hydrophilic") hemispheres. More precisely, the protein is parametrized by a line of length L, and the interactions between lipids and proteins are given by Eq. (5) with σ 0 = σ t as before, where r is now the shortest distance between the center of a lipid bead and the protein line. The interactions between two proteins are purely repulsive and given by (4) with σ 0 = 2 σ t and r the minimum distance between the two interacting protein lines.
Even though these cylinders are allowed to tilt, they basically stay aligned with the membrane normal in our L [σt] pt = 3 pt = 6 4 14.1 ± 2.1 12.4 ± 2 6 12.5 ± 0.8 11.9 ± 2.2 8 11.1 ± 1.2 7.74 ± 1.15 model bilayers. The distribution of solid angles assumes a maximum at tilt zero, and the average tilt angles α for the parameters used in the present work are less than 15 degrees, as shown in Table I. These tilt angles should be regarded as valid in the general case, since we are not aiming at modeling any specific protein. Our results are compatible with an earlier simulation study by Venturoli et al. [32], where the average tilt angle of proteins with both negative and positive mismatch of up to 40% was also only slightly affected by the mismatch. Only when a considerable positive mismatch of more than 70% was present, a significant increase in average tilt angle for model proteins of small aspect ratio (diameter to length) was observed.
Since a protein is "rough" on an atomic scale, representing a complex structure like a protein as a simple smooth cylindrical object may seem to be a rather crude approach. Nevertheless, our way of modeling the proteins can be justified by the fact that e.g. α-helices are packed with vanishing little free space within the helices [33] and are therefore fairly smooth on the scale of ∼ 10Å. There are no large cavities into which chains or even whole molecules would fit [34].

C. Simulation Method
We have carried out Monte Carlo simulations at constant pressure P = 2.0 /σ 3 t and constant temperature T with periodic boundary conditions in a simulation box of variable shape and size. Thus, we are performing Monte-Carlo simulations in an N P T ensemble with effective Hamiltonian where U is the interaction energy, V the volume of the simulation box, V 0 an arbitrary reference volume and N the total number of beads. In practice, three main types of Monte-Carlo moves were proposed and then accepted or rejected according to a Metropolis criterion, namely 1) translational local moves of the lipid beads, 2) protein translation and rotation moves, and 3) global moves which change the volume of the simulation box or its shape. Most of these moves have been discussed earlier [24]. Here we only sketch the algorithm for the new protein rotation moves. The new trial direction d p of the protein is generated in several steps. First a vector u, which lies randomly distributed on a unit sphere, is generated. Choosing a random vector on the surface of a unit sphere is efficiently done by applying an algorithm proposed by Marsaglia [35]: Two random variables r 1 and r 2 within the interval (−1, 1) are generated and ζ 2 = r 2 1 + r 2 2 is calculated. If ζ 2 > 1 the random numbers are discarded and a new pair r 1 , r 2 has to be generated. For ζ 2 < 1 the components of the random unit vector in Cartesian coordinates are given by If the new direction of the protein was chosen to be u, the acceptance rate of the move would be prohibitively small. Therefore, u is scaled down by ∆ tilt and the vector t = ∆ tilt u is added to the tip of the normalized direction d p of the protein. Finally, normalizing the resulting vector d p = d p + t to d p gives the trial direction of the protein used to calculate the energy change for the Metropolis criterion. By varying the scaling factor ∆ tilt during the "prerun" the acceptance rate of the tilt moves can be adjusted to the desired rate. Effectively, this means adjusting the opening angle of the cone around d p , which limits the maximum tilt angle θ p,max of the move.
During one Monte-Carlo step (MCS) there is on average one attempt to move each bead and protein, and one attempt to rotate the protein. Since the global moves require rescaling of all particle coordinates, which is rather expensive from a computational point of view, they are performed only every 50th MCS on average.

D. Elastic Theory
In the following, we briefly sketch the elastic theory with which we compare the simulation results [10,14,18], and which has proved to describe very well the properties of our model membranes in the fluid state [18]. The membrane is treated as a system of coupled monolayers, which fluctuate subject to the constraint that the volume of lipids is locally conserved. It is also considered to be almost flat, i.e., the positions of both monolayers can be parametrized by single-valued functions h i (x, y).
Here we are mainly interested in the thickness deformations, which we characterize by the locally smoothed deformation profiles Φ el (x, y) of monolayers about the unperturbed value t 0 ( i.e., the total monolayer thickness is t 0 + Φ el ). The free energy of monolayer thickness deformations is then written as [15] Parameter L Value kc where k c and k A are the bending and compressibility modulus of the bilayer, c 0 is the spontaneous curvature, ζ a parameter related to the spontaneous curvature, and k G the Gaussian curvature. All these parameters have been determined for our model membranes in earlier work from analyses of the membrane fluctuations and stress profiles [18]. We note that the terms c 0 and k G can be rewritten as pure surface terms. Protein inclusions distort the membranes by imposing a fixed thickness at their boundaries [10], Φ el = t R , and may introduce additional surface fields, which enter in the same way as the spontaneous curvature term [18] and thus effectively renormalize c 0 . Hence the effect of proteins on the membrane can be characterized by the two parameters t R and c 0 , the renormalized curvature. They have been determined for our model in earlier work [18] from the distortion profiles of single membranes about infinitely long straight cylinders. Table II summarizes the values of all elastic parameters for the simulation models studied in this work. The elastic free energy, Eq. (10), can be minimized analytically for membranes containing only a single cylindrical protein. For two proteins, it must be minimized numerically. To this end, the free energy integral was discretized in real space using a square grid of spatial discretization h = 0.25 σ t , a system size of 50 × 50 σ 2 t and a second-order difference scheme to evaluate the derivatives. The boundary condition was implemented by setting φ = t R inside the inclusion. The minimization was done via a steepest descent method, using a relaxation scheme [36]. This procedure gave deformation profiles with which we could compare the simulation results (see below).

A. Lipid-mediated interactions between inclusions
We now turn to discussing our simulation results for the membrane-induced interactions between cylindrical inclusions in the bilayer. In this subsection, we focus on the effect of protein orientation fluctuations on the effective interactions between the model proteins, i.e., the PMF.
The radial distribution function g(r) as a function of the protein-protein distance r was obtained from simulation runs using the technique of successive umbrella sampling [37] combined with a reweighting procedure. As starting configurations we used equilibrated systems with 750 to 760 lipids and two simple transmembrane proteins of diameter 3 σ t . A first estimate of g(r) was obtained during pre-runs of length 2×10 6 MCS. Then, biased runs of 3 × 10 6 MCS were performed to improve the statistics of configurationally less frequent protein-protein distances. In the case of tiltable proteins, 10 × 10 6 MCS had to be performed to achieve comparable accuracy in the histogram data.
After removing the bias from these results and combining the overlapping distributions the effective potential w(r) = −k B T ln g(r) can be extracted. At small inclusion-inclusion distances direct interactions of proteins and, in the case of long proteins, additional depletion induced attraction due to the solvent particles come into play. Since our main interest lies in the study of lipid-mediated medium and long ranged interactions, these parts of the curves have been cut off. Our sampling procedure resulted in a statistical error of about ±0.3 within each umbrella window. The distance range between 4 and 10 σ t was typically covered by 7-8 windows (with larger windows at larger distances), which gives an accumulated error of ±2 in the first minimum of the PMFs in the worst case, and of at most ±1 at distance 6 σ.
The hydrophobicity of the protein, i.e., the interaction strength between the hydrophobic core of the membrane and the hydrophobic part of the inclusion, is tuned by the parameter pt . Only for sufficiently high values of pt do the model proteins locally distort the bilayers. In the following we will consider two cases, weakly hydrophobic proteins with pt = 3.0 and strongly hydrophobic proteins with pt = 6.0 . In the first case lipid-protein contact is approximately energetically equivalent to lipidlipid contact, whereas in the second case lipid-protein contact is highly preferred. Weakly hydrophobic proteins only induce small perturbations of the lipid environment, strongly hydrophobic proteins lead to strong deformations.
Furthermore, the hydrophobic length L has been varied from L = 4 σ t to L = 8 σ t . The hydrophobic thickness of the membrane is 6 σ t , hence L = 4 σ t corresponds to a negatively mismatched protein, and L = 8 σ t to a positively mismatched protein.
The simulation results are compiled in Figs. 2, 3, and 4. The left and right graphs in each figure correspond to straight and tiltable cylinders, respectively. One immediately notices the first effect of orientational fluctuations: Due to lipid packing in the vicinity of the inclusions the curves feature oscillations with a wavelength of approx- Potential of mean force of two isolated proteins as a function of the protein distance at medium (solid) and strong (dashed) hydrophobicity pt for proteins with negative hydrophobic mismatch (L = 4 σt). The proteins of the left graph were always aligned parallel to the z axis of the simulation box, whereas proteins in the right graph were allowed to tilt.  imately 1 σ t , i.e. the diameter of the lipid tail beads. The higher pt the more pronounced the lipid layering. This holds both for straight and tiltable proteins. However, the oscillations are significantly reduced in the latter case. Hence orientation fluctuations reduce the effect of lipid layering on the protein-protein interactions.
Next we discuss separately the PMFs for strongly and weakly hydrophobic proteins.

Strongly hydrophobic proteins
A general feature that straight and tiltable protein have in common at pt = 6.0 is an enhanced tendency to aggregate for hydrophobically mismatched proteins, both for positive or negative mismatch (dashed curves in Figs. 2 and 4). The effective interactions between mismatched proteins have an attractive long-range contribution which is not detectable for hydrophobically matching proteins (Fig. 3). Compared to straight proteins, however, the range of this additional attraction is reduced for proteins with orientational fluctuations, and the attraction essentially vanishes within the error at distances beyond r ∼ 7 σ t .
The first and deepest attractive minimum of the PMFs slightly below r = 4 σ t lies at approximately the same energy for both types of models, i.e. between −17 and −15 for proteins of length L = 4 σ t and 6 σ t and about −40 for proteins of length L = 8 σ t . This minimum is related to lipid bridging and a consequence of the strong attraction between lipids and proteins. At distance r ∼ 4 σ t , one layer of lipids is in contact with both proteins, thus stabilizing the conformation at that distance. The lipid bridging energy increases with increasing hydrophobic contact area between lipid tails and proteins. Compared to hydrophobically matching proteins, it is thus reduced for negatively mismatched proteins and enhanced for positively mismatched proteins. In addition, both negatively and positively mismatched proteins have the above-mentioned attractive interactions due to the elastic deformation of the layer. Compared to hydrophobically matching proteins, the contributions of lipid bridging and hydrophobic mismatch are opposite for negatively mismatched proteins, and cumulative for positively mismatched proteins. As a result the first minimum has approximately the same energy for proteins with negative hydrophobic mismatch and matching proteins, whereas proteins with positive mismatch show the strongest attraction.

Medium hydrophobicity
For weakly hydrophobic proteins with pt = 3.0 , the effect of hydrophobic mismatch at short distances is much less pronounced than in the case of strongly hydrophobic proteins. This is especially true for proteins with fluctuating orientations. In the latter case, the first minimum has almost the same depth for all types of mismatch, and the PMFs for different L are not further apart than 3 (solid curves in the right graphs of Fig. 2, 3, and 4).
At intermediate distances, however, the interaction between proteins with fluctuating orientations is found to depend crucially on the type of mismatch: For positively mismatched proteins, it has a long range attractive contribution similar to that observed for strongly hydrophobic proteins. For negatively mismatched proteins, this contribution is still present at distances r < 7 σ t , but then it turns around and gives way to a repulsive interaction. The repulsive barrier is small, but significant within the error of ±0.5 at distance 8 σ t . The behavior for matching proteins is intermediate -the layering interactions are superimposed by an attractive contribution up to r ∼ 6 σ t , which disappears within the error at larger distances. We note that this behavior is in marked contrast both to the behavior of strongly hydrophobic proteins and of weakly hydrophobic, but straight inclusions. In the latter cases, one observes a purely attractive long-range hydrophobic mismatch interaction, which vanishes for hydrophobically matched proteins.
In fact, the elastic theory does predict a repulsive barrier in the interaction free energy of two transmembrane proteins for positive and negative hydrophobic mismatch [18]. The maximum of this small repulsive barrier should be found in the range of r ∼ 6 σ t to 8 σ t and can be associated with a peak, i.e., a soft mode, in the spectrum of thickness fluctuations of the pure bilayer. One may speculate that the repulsive shoulders observed here are caused by the same effect. However, our simulations results show no repulsive barrier for proteins of length L = 8.0 σ t . The attractive range for positively mismatched proteins rather extends up to a protein-protein separation of about 9 σ t and levels off without any observable repulsive region. Furthermore, the attractive behavior of hydrophobically matching proteins is not predicted by the elastic theory either.
Other comparable simulations on tiltable proteins have been performed by Schmidt et al. Both for negative and positive hydrophobic mismatch they observe a longrange, lipid-mediated attraction. In the case of zero mismatch this interaction is less attractive [38]. Further, their PMFs exhibit an oscillating "fine structure", which they attribute to the discreteness of the membrane.
In another study, de Meyer et al. have compared PMFs for proteins which were allowed to tilt, and proteins which were always aligned parallel to the z-axis [17]. Their results for proteins comparable to our capped cylindrical inclusions of diameter D p = 3 σ t can be summarized as follows: In the case of positive mismatch the long-range interaction between two proteins should be influenced by the degree of protein tilt. If the protein was not allowed to tilt, they found an attractive region for small distances, followed by a small repulsive barrier and a shallow and broad minimum at large protein separation. If tilting of proteins was allowed, only a small, but broad repulsive barrier at medium inter-protein dis-tances was found after the typical attractive range at small distances for both positive and negative mismatch. The shallow attractive region at larger distance has vanished. In our model of long cylindrical proteins without tilt a shallow repulsive barrier at medium distances was also present, but only at very weak hydrophobicity of pt = 1.0 [18]. Other effects due to lipid layering were much less pronounced in the study of de Meyer et al., since they were using soft, merely repulsive potentials. But at very small protein separations they did observe slight oscillations in the PMFs, which they assigned to the free energy needed to remove lipids in between the proteins.

B. Profiles around two proteins
After having discussed the effective lipid-mediated interactions between two transmembrane inclusions, we now turn to the closer investigation of structural lipid rearrangements around two model proteins at small and medium separation. In this section, we will consider strongly hydrophobic proteins with pt = 6 .
As we have seen in the previous section, one characteristics of our lipid model with its hard-core, Lennard-Jones type interactions is the presence of lipid packing phenomena, which are usually much less prominent in systems modeled with soft dissipative particle dynamics potentials [18]. Especially in the area surrounding proteins one can expect ordering of the lipids, which is rather distinct from the situation further away from inclusions in the lipid bulk, and which may affect the interaction with other inclusions. Here, two-dimensional profiles give useful insight into the structure of membrane-characterizing quantities. To determine them, we have simulated equilibrated bilayers containing two capped proteins, keeping the positions of the proteins fixed, but allowing for tilt fluctuations. As before, we compare systems with cylindrical proteins of varying hydrophobicity, hydrophobic mismatch and protein-protein distance.
The two-dimensional profiles can be separated into two regions: In the inner region between the two inclusions, the profiles reflect the combined effect of both proteins. In the outer region, the influence of one protein is geometrically screened, and the profiles mainly reproduce the perturbations induced by one single protein.

Areal tail bead density
First, the areal density of tail beads χ t around two proteins is investigated. In the following graphs, the values are always given with respect to one monolayer, i.e. the average number of tail beads obtained for each bin lying in the xy plane was divided by 2. Both length and height of a single bin is set to 0.5 σ t . Therefore each sampled cell has an area of 0.25 σ 2 t only. At this resolution interesting details of the lipid structure in the vicinity of the inclusions are revealed.
Not surprisingly, negative hydrophobic mismatch, which goes along with thinning of the surrounding bilayer, leads to a reduction of the areal tail density, whereas positive hydrophobic mismatch enhances the areal tail density along with the membrane thickness. The fine structure exhibits additional features. In the following, we will focus on hydrophobically matching proteins of length L = 6 σ t . The layering of the lipids around the proteins, which gave rise to the oscillatory behavior of the PMFs in the previous section, manifests itself clearly in concentric rings with enhanced (light) and reduced (dark) areal tail bead density (Fig. 5). The variations in areal density are usually in the order of 20% compared to the average value at larger distance from the inclusions. If the proteins are close to each other (Fig. 5, left), they pin the lipids entirely i.e., the areal tail density exhibits sharp local peaks at two well-defined positions between the two proteins.

Thickness profiles
Next, we inspect the thickness profiles around two cylindrical inclusions at fixed distances. Specifically, we consider short distances, slightly below a protein-protein separation of ∼ 4 σ t , medium distances at a separation of ∼ 6.5 σ t , and larger separations lying between 8 and 9 σ t . A regular grid of cells with an area of ∆x∆y = 1.0 σ 2 t was used to subdivide the system. Averaging of the membrane thickness was done within each quadratic bin thus obtained. If necessary, the systems were rotated such that the connecting line between the centers of the proteins was parallel to the x-axis for better comparison and visualization purposes. Light shading indicates thick membranes, dark shading thin membranes. For comparison, we have also calculated the same profiles with the  (Table II). The results are shown in Figs. 6, 7, and 8. Both for small and medium protein-protein separations the perturbations induced by the proteins overlap. This amplifies the reduction or increase of thickness which would be caused by a single protein of negative or positive hydrophobic mismatch, respectively. At large distances, the proteins are surrounded by separate perturbation shells. Fig. 6 shows the results for proteins with negative mismatch. For small and medium protein-protein distances the bilayer thickness is smallest in the region between the proteins (Fig. 6, top left and middle left), which is in remarkable agreement with the predictions of the elastic theory (Figs. 6, top right and middle right). For proteinprotein distances 5 − 7σ t , the elastic theory even predicts a certain amount of "undershooting", i.e., the membrane thickness on the axis between the two proteins is predicted to be smaller than at the protein surface. This is also found in simulations, albeit to a much lesser extent.
Proteins with positive hydrophobic mismatch (L = 8.0 σ t ) are surrounded by a clearly thickened bilayer (Fig. 7). As for negatively mismatched proteins, the membrane perturbations are enhanced in the region between the proteins. For protein distances 5 − 7 σ t , the elastic theory predicts "overshooting", i.e., the thickness is expected to be largest at the point right between the two proteins. However, the simulation data do not show this effect. In this respect, positively mismatched proteins show a different behavior than negatively mismatched proteins. Finally, the bilayer thickness between hydrophobically matching proteins (L = 6.0 σ t ) remains practically unaltered as expected, except in the very close vicinity of the proteins, where a small annulus of enhanced thickness develops (Fig. 8, left). Elastic theory predicts a slightly enhanced thickness in the vicinity of the proteins (Fig. 8,  right).

IV. SUMMARY AND DISCUSSION
In summary, we have studied protein-protein interactions in model lipid membranes, using two types of model proteins: Infinitely long cylinders with fixed orientation, and finite spherocylinders with fluctuating orientation. In general, the interactions are dominated by packing effects and smooth attractive hydrophobic mismatch interactions. Besides the geometric mismatch the hydrophobicity of the proteins is a key quantity.
The comparison of the PMFs for the two protein types allowed to assess the effect of orientation fluctuations. It can be summarized as follows: The main and most prominent effect of orientation fluctuations is to reduce the lipid packing effects on the PMFs, i.e., the corresponding short range oscillatory behavior. Furthermore, they slightly reduce the range of the hydrophobic mismatch interactions. For weakly hydrophobic proteins, they may also affect qualitatively the shape of the potentials: In the case of negatively mismatched proteins, orientation fluctuations introduce a weak repulsive maximum in the PMFs which is not present for orientationally fixed proteins. This is consistent with earlier results by de Meyer et al. [17], using a different protein model with hydrophobic interactions that would be considered 'weak' in our context. The direct experimental measure of interaction potentials between membrane proteins is still a challenging task due to the required temporal and spatial resolution. The "protein" diameter in our study roughly corresponds to that of gramicidin, for which the effect of hydrophobic mismatch interactions was verified experimentally in 1999 by Harroun and coworkers [39]. The hydrophobic mismatch effect has also been investigated by systematic studies of synthetic model peptides [2,40]. In these studies, the evidence for mismatch-induced interactions was mainly based on the observation of protein clustering. Recently, spatially resolved interaction energies be-tween proteins of positive hydrophobic mismatch have been calculated for mobile ATP-synthase c-rings in membranes basically consisting of phosphatidyl glycerophosphate from the probability distribution of the centerto-center distance by Casuso et al. [41]. After a softcore short range repulsion attributed to the structural perturbation of the lipid/protein organization, a minimum could be observed in the interaction energy. The authors consider this attractive potential to be caused by membrane-mediated interactions due to the deformations of the lipid bilayer. The ATP-synthase c-rings of diameter 65Å ± 5Å would correspond to proteins of diameter 10 σ t in our systems. Therefore, direct comparison with our simulation data is not possible. Nevertheless, the experimental finding of protein-protein attraction by mismatch-induced deformation is in accordance with the resulting smooth attractive interaction in our simulations.
For the case of strong interactions, we have also studied the distribution of lipid tails in systems containing two proteins and found that two proteins can form complexes with surrounding lipids, pinning these to their surface. We speculate that such effects may also occur in biological environments. Furthermore, we have compared the thickness profiles in systems containing two proteins with the elastic theory of two coupled monolayers. We found remarkable agreement between the theoretical prediction and the simulation results, with only slight qualitative deviations at intermediate protein distances ∼ 7 σ t : here, the elastic theory predicts "over"-or "undershooting" effects which are much weaker or not discernible at all in the simulations. The theoretical prediction for the PMFs has not been shown here (see [18]), because it is difficult to compare with the simulation data due to the lipid packing effects. Regarding the smooth hydrophobic mismatch part, the theory is in rough qualitative agreement with the simulations, except that its predicts a repulsive barrier at distances ∼ 8 σ t which is not observed in the simulations. The discrepancy is possibly related to the slight deviations between theory and simulations in the thickness profiles mentioned above.