Insertion of protein molecules into a lipid bilayer perturbs the packing of nearby lipid chains, since the protein-lipid and lipid-lipid interactions are generically different. Such local perturbation leads to entropy loss of the lipid chains and induces a short-range interaction between the proteins. By defining an order parameter for lipid chain orientation, Marčelja (1976) constructed a mean-field Hamiltonian for a model system of two hexagonal proteins embedded in a flat lipid bilayer and predicted a pure attraction between the proteins. By taking the same assumption that fluctuations of lipid orientation are suppressed in the protein vicinity, Schröder (1977) derived an expression for the attraction. The attraction arising from lipid orientational entropy was first verified by coarse-grained Monte Carlo simulations (Sintes and Baumgärtner, 1997), as will be discussed in Section 3.1. May and Ben-Shaul (2000) applied detailed molecular chain packing theory to calculate the interaction between two protein walls in a bilayer and found that the interaction is attractive at small separations and repulsive at intermediate separations. This interaction starts to level off at separation around the hydrophobic thickness of the proteins or bilayer. Nonmonotonic lipid-mediated interaction potentials between protein inclusions were also reported by Lagüe et al. (2001) based on integral equation theory for liquids.

Protein inclusions generally exhibit different hydrophobic thickness from the embedding lipid bilayer. Lipid chains surrounding the inclusions will stretch or compress in order to avoid or alleviate the exposure of hydrophobic regions of the proteins or lipids. This deformation represents another source of free energy cost that contributes to the bilayer-mediated short-range interaction. Owicki and McConnell (1979) presented a Landau-de Gennes free energy to account for lipid chain deformation and bilayer area change caused by hydrophobically mismatched inclusions, and obtained a short-range attraction between such two inclusions. Huang (1986) formulated a continuum theory to describe the bilayer deformations around a rigid inclusion by considering the free energies associated with monolayer bending, lipid chain compression and surface tension, and reported a nonmonotonic bilayer thickness profile around single inclusions. Dan et al. (1993) and Aranda-Espinoza et al. (1996) adopted this continuum theory to calculate the membrane-mediated interaction between two cylindrical inclusions with hydrophobic mismatch. The interaction looks qualitatively similar, in the case of vanishing spontaneous curvature of the monolayers, to that obtained in Ref (May and Ben-Shaul, 2000), and has a range of about two to three times of the bilayer thickness.

The continuum approach (Huang, 1986; Dan et al., 1993; Aranda-Espinoza et al., 1996) based on membrane elasticity has been widely used to investigate the role of membrane-mediated interactions in such membrane protein processes as formation of gramicidin ion channels (Huang, 1986; Bitbol et al., 2012; Bories et al., 2018), cooperative gating of mechanosensitive channel of large conductance (Ursell et al., 2007; Haselwandter and Phillips, 2013; Kahraman et al., 2016a) and assembly of chemoreceptor trimers (Haselwandter and Wingreen, 2014). Extensions of the membrane elasticity model have been made by including lipid tilt (Fournier, 1999; Bohinc et al., 2003), Gaussian curvature (Brannigan and Brown, 2007), gradient of bilayer thickness (Bitbol et al., 2012), and asymmetry in two monolayers due to noncylindrical shape of the inclusions (Argudo et al., 2017). We briefly describe this approach for a up-down symmetric and single-component lipid bilayer. As shown in Figures 1A,B, the thickness deformations of the lipid bilayer around a hydrophobically mismatched protein inclusion are characterized by the relative displacement u(x, y) of the upper monolayer with respect to the horizontal midplane. By a Taylor-expansion around the unperturbed flat state of the bilayer with thickness 2u ₀ and area per lipid Σ₀, the monolayer free energy can be expressed in terms of u, gradient of u (i.e., ∇u), mean curvature H ≈ ∇² u/2, and Gaussian curvature K ≈ det(∇∇u). The free energy of monolayer compression or stretching per projected area is f c = K A ( u / u 0 ) 2 / 4 with K _A the bilayer’s area compression modulus. The corresponding surface-tension term is f _t = σ[u/u ₀ + (∇u)²/2]/2 (Huang, 1986; Haselwandter and Phillips, 2013; Kahraman et al., 2016b) with σ the bilayer tension. The bending energy density assumes f b = [ κ ( ∇ 2 u ) 2 / 2 + κ c 0 ∇ 2 u + κ ( c 0 − c 0 ′ Σ 0 ) ( u / u 0 ) ∇ 2 u + κ ̄ det ( ∇ ∇ u ) ] / 2 (Dan et al., 1993; Aranda-Espinoza et al., 1996; Brannigan and Brown, 2007), where κ is the bilayer bending rigidity, c ₀ the monolayer spontaneous curvature, c 0 ′ = ( ∂ c 0 / ∂ Σ ) Σ 0 the change of c ₀ due to lipid area variation, and κ ̄ the bilayer Gaussian modulus. The total free energy of the perturbed monolayer per inclusion is then given by the functional F [ u ] = ∫ ∫ d x d y ( f c + f t + f b ) . Minimization of F [ u ] under appropriate boundary conditions (Nielsen et al., 1998; Nielsen and Andersen, 2000; Brannigan and Brown, 2006) determines the bilayer deformations around the inclusions.

There exist two types of long-range interactions between protein inclusions mediated by the embedding fluid membranes: fluctuation-induced interactions and curvature-induced elastic interactions. As the name suggests, they originate, respectively, from modification of membrane fluctuations and from perturbation of the equilibrium membrane shape by the presence of protein inclusions. In theoretical considerations, the fluid membrane, on length scales larger than its thickness, is coarse-grained into a two-dimensional (2D) elastic surface governed by the Helfrich Hamiltonian (Helfrich, 1973) H el = ∫ [ κ ( 2 H − c 0 ) 2 / 2 + κ ̄ K + σ ] d A , where H is the membrane’s mean curvature, K the Gaussian curvature, c ₀ the spontaneous curvature due to bilayer asymmetry, κ the bending rigidity, κ ̄ the Gaussian modulus, and σ the lateral tension conjugate to the membrane area. Typical values for the physical properties of synthetic and biological membranes are κ ∼ 10k _B T ≈ 4 × 10⁻²⁰J, − κ ≤ κ ̄ ≤ − 0.7 κ (Deserno, 2015), and σ ∼ 1 μN/m (Simson et al., 1998). ℓ ≡ κ / σ ∼ 200 nm is the characteristic length scale over which the surface tension dominates over the bending energy. A membrane surface that exhibits small deformations is parameterized by using a displacement field with respect to a planar, spherical, or cylindrical reference surface; see Figure 1C. The Helfrich Hamiltonian is then a function of the displacement field. Proteins that are embedded in or attached to the membrane are treated by applying boundary conditions to the displacement field or via extra energy terms that enter the Hamiltonian of the system. The fluctuation-induced interaction can be estimated from cumulant expansion, whereas the curvature-induced elastic interaction is often derived by finding the minimum-energy shape of the membrane under appropriate boundary conditions and additional requirements that the inclusions are in mechanical equilibrium, i.e., force- and torque-free. In all the following, we will use the notations F fl and F el in energy units to distinguish the two types of long-range interactions, and cite the leading-order expressions for F fl and F el as derived in literature unless specified otherwise.

Monte Carlo (MC) (Sintes and Baumgärtner, 1997; West et al., 2009) and molecular dynamics (MD) (Venturoli et al., 2005; de Meyer et al., 2008; Schmidt et al., 2008; West et al., 2009) simulations based on coarse-grained models of lipids and proteins have been used to investigate membrane-mediated short-range interactions between two cylindrical protein inclusions with emphasis on the generic feature of the interactions. Simulation study of proteins without hydrophobic mismatch was first done by Sintes and Baumgärtner (1997), who considered two rigid cylinders (of diameter two to four times of lipid width σ) embedded in a lipid bilayer whose bending deformation was strongly suppressed. They found a depletion-induced attraction for inclusion separation R < σ and an oscillating interaction for σ < R < 6σ attributed to inhomogeneous distribution and orientational fluctuations of lipid chains around inclusions. In simulation studies with hydrophobic mismatch between the bilayer and protein inclusions, the bilayer thickness profile around single inclusions was found to exhibit similar nonmonotonic behavior (Venturoli et al., 2005; de Meyer et al., 2008; Schmidt et al., 2008; West et al., 2009), consistent with the elastic theory. The potential of mean force (PMF) between two rigid inclusions depends on inclusion size (de Meyer et al., 2008) and lipid-protein interaction (West et al., 2009). For protein inclusions that have no or weak affinity to lipid chains, the PMF is attractive at smaller separations and repulsive at intermediate separations (de Meyer et al., 2008; Schmidt et al., 2008), also in accordance with the elastic theory. For protein inclusions strongly attracting lipid chains, the PMF is highly oscillatory with a repulsion at close separations (West et al., 2009). Protein inclusions with hydrophobic thickness larger than that of the bilayer can be tilted (Venturoli et al., 2005; Klingelhoefer et al., 2009) or even bent (Venturoli et al., 2005) in the membrane in order to avoid exposure of their hydrophobic domains.

All-atom simulations can provide atomistic details of the bilayer deformation around single protein inclusions. Kim et al. (2012) performed all-atom MD simulations to measure the thickness profile of different bilayers embedding a Gramicidin A channel, and found qualitative discrepancy from theoretical predictions, probably due to specific protein-lipid interactions that are not addressed in the simplified membrane elasticity model. Mondal et al. (2011) studied the energetics of lipid bilayer deformations around a noncylindrical protein inclusion like G-protein coupled receptors by taking the continuum theory where the membrane-protein boundary conditions were extracted from atomistic MD simulations. Argudo et al. (2017) revisited the atomistic MD simulations of a Gramicidin A channel embedded in a POPC bilayer and showed that the membrane deformations and tilt of the ionic channel can be quantitatively captured by a refined bilayer model that incorporates the chemistry and geometry of the protein inclusions. However, due to the computational cost, it remains challenging to measure the membrane-mediated short-range interactions between two protein inclusions from all-atom MD simulations.

To directly measure the fluctuation-induced interactions remains a computationally difficult task, since they are weak and often coupled with the curvature-induced elastic interactions. Weikl (2001) studied the fluctuation-induced aggregation of protein inclusions much more rigid than the fluid membrane via MC simulations, where the membrane was represented by a discretized 2D elastic sheet and protein inclusions occupy single vacant sites on the membrane surface. The systems with inclusions of κ _p and κ ̄ p two order-of-magnitude larger than those of the membrane were found to separate into inclusion-rich and inclusion-poor phases even in the absence of any direct protein-protein attraction, whereas the systems with no contrast in the Gaussian moduli ( κ ̄ p = κ ̄ ) exhibit the same critical point as if the membrane were completely flat without shape fluctuations. This finding points out that a difference in Gaussian moduli is necessary for a fluctuation-induced interaction between membrane inclusions, as mentioned in Section 2.2.1. Pezeshkian et al. (2016) reported from MD simulations the clustering of rigid pentagon shaped nanoparticles, coarse-grained model of bacterial Shiga toxin, on lipid membranes driven by the fluctuation-induced attraction. Very recently, Sadeghi and Noé (2021) extracted from MD simulations the membrane-mediated effective interactions between protein particles embedded in a fluid membrane modeled by particle-based elastic sheet. The interaction varies non-monotonically with interparticle separation and has a depth of about k _B T for different values of protein stiffness, which can not be accounted for by the sum of the two-body interactions F fl + F el ≈ − 6 k B T ( a / R ) 4 + 4 π κ ( α 1 2 + α 2 2 ) ( a / R ) 4 . It is not clear whether the discrepancy is model specific.

MD simulations with coarse-grained models of lipid membranes at different levels of resolution (Reynwar et al., 2007; Olinger et al., 2016; Xiong et al., 2017; Spangler et al., 2018; Noguchi, 2016; Noguchi and Fournier, 2017), MC simulations (Bahrami et al., 2012; Šarić and Cacciuto, 2012b,a; Vahid et al., 2017; Bahrami and Weikl, 2018; Bonazzi and Weikl, 2019) and numerical minimizations (Reynwar and Deserno, 2011; Schweitzer and Kozlov, 2015) based on mesoscopic elastic surface models showed that proteins or particles adhering to membranes experience curvature-induced interactions, which are strongly attractive in many cases and can drive particle assembly on the membranes. We first review the studies of spherical particles. Reynwar et al. (2007) computed directly from coarse-grained MD simulations the force between two capsids adhering strongly to a lipid bilayer, and obtained repulsive forces at small capsid separations followed by attractive ones at large separations as shown in Figure 4A, seemingly contradictory to the theoretical prediction of pure elastic repulsion for conical inclusions in Section 2.2.1. By using the Surface Evolver package (Brakke, 1992) to numerically minimize the energy of a membrane adhering to two spherical particles, Reynwar and Deserno (2011) confirmed later that the curvature-induced force is indeed repulsive for small contact angles (i.e., weak adhesion), whereas for contact angles larger than 90° (i.e., strong adhesion), the force changes from repulsive to attractive with increasing separation, consistent with their MD results in Ref (Reynwar et al., 2007).

Using a triangulated surface model for vesicle membranes, Bahrami et al. (2012) determined from simulated annealing MC simulations the minimum-energy shape of vesicle membranes interacting with adhesive spherical particles, and discovered stable membrane tubules that wrap one row of two or three particles; see Figure 4B. As shown in Figure 4C, similar membrane tubular structures were observed in constant-temperature MC simulations by Šarić and Cacciuto (2012b), who also reported linear aggregation of spherical particles adsorbed on vesicle membranes (Šarić and Cacciuto, 2012a). MD simulations (Reynwar et al., 2007; Xiong et al., 2017; Spangler et al., 2018) with molecular models for flat lipid bilayers showed that linear aggregation of adsorbed spherical particles induces membrane tubulation and vesiculation. These studies point towards curvature-induced strong attractions between spherical particles adhering to fluid membranes. Vahid et al. (2017) further demonstrated that, in the case of quasi-spherical vesicles, the curvature-induced attraction between two adsorbed spherical particles becomes weaker as the vesicle gets bigger. In the case of quasi-ellipsoidal vesicles, the attraction for two particles placed along the major axis is different from that along the minor axis, and their relative magnitude depends on the ellipticity of the vesicles; see Figure 4D. This MC simulation result suggests that the background curvature of closed membranes affects the curvature-mediated attraction. In the extension of their previous study (Bahrami et al., 2012), Bahrami and Weikl (2018) reported that two spherical Janus particles with one side strongly adhering to vesicle membranes can attract or repel each other, depending on the area fraction of the adhesive side and on the shape (i.e., concave or convex) of the adhering membrane segments; see Figure 4E. In all these simulated systems (Reynwar et al., 2007; Olinger et al., 2016; Xiong et al., 2017; Spangler et al., 2018; Noguchi, 2016; Noguchi and Fournier, 2017; Bahrami et al., 2012; Šarić and Cacciuto, 2012b,a; Vahid et al., 2017; Bahrami and Weikl, 2018; Bonazzi and Weikl, 2019; Reynwar and Deserno, 2011; Schweitzer and Kozlov, 2015), particle-membrane adhesion energy, membrane bending energy, and possible constraint due to conservation of the volume enclosed by membranes determine together the optimal membrane shape and thus how the curvature-induced interaction varies with inter-particle distance.

Numerical studies of anisotropic protein inclusions or scaffolds interacting with fluid membranes have also revealed that the membrane-mediated attractions are important for the protein to assemble and to remodel the membranes. Using Surface Evolver to find the minimum-energy shape of membranes interacting with two rigid protein scaffolds, Schweitzer and Kozlov (2015) showed that for circular scaffolds with anisotropic curvature (saddle-like or ellipsoidal shape), or isotropically curved scaffolds of elongated shapes (noncircular footprint on the membrane), the curvature-induced interaction is repulsive at small inter-scaffold separations and attractive at large separations. Specifically, the curvature-mediated attraction between two BAR-domain-like scaffolds was found to be very strong; see Figure 5A. Using a coarse-grained molecular model of N-BAR domain, Simunovic et al. (2013) showed from MD simulations that the proteins assemble on flat or vesicle membranes at low concentrations and form a mesh of linear aggregates as shown in Figure 5B. Noguchi (2016), and Noguchi and Fournier (2017) studied the assembly of arc-shaped proteins on flat membranes with coarse-grained MD simulations, and found side-by-side alignment of the proteins around membrane tubules; see Figure 5C. The membrane-mediated side-by-side arrangement was also reported by Bonazzi and Weikl (2019) in MC simulations of arc-shaped particles remodeling an initially spherical vesicle modeled by a triangulated surface; see Figure 5D, where the rather loose arrangement of particles has been experimentally reported for N-BAR proteins interacting with membrane tubules (Daum et al., 2016). We note that, in addition to the membrane-mediated indirect interactions, direct protein-protein interactions may also play a part in the assembly of anisotropic proteins on membranes. For instance, the helical arrangement of N-BAR proteins around membrane nanotubes found in coarse-grained MD simulations (Simunovic et al., 2016) is very likely due to the direct attraction between the proteins. Using the same discrete model as in Ref (Weikl, 2001) and treating saddle-like inclusions as point-like constraints that impinge anisotropic curvature on the membrane, Dommersnes and Fournier (2002) simulated the assembly of those inclusions, assisted by the curvature-induced attraction, into regular arrays that shape the membrane into the experimentally observed egg-carton pattern as shown in Figure 5E.

It is worthy to mention the membrane-mediated interactions in the systems of cell adhesion that is mediated by the specific binding of membrane-anchored receptors and ligands. The receptor-ligand complexes constrain the local separations of the two adhering membranes and thus experience fluctuation-induced attractions (Bruinsma and Pincus, 1996; Krobath et al., 2009; Weikl, 2018). An important biological consequence of these membrane-mediated attractions is the cooperative binding of cell adhesion proteins, as corroborated by MD simulations (Hu et al., 2013; Hu et al., 2015; Xu et al., 2015) and experiments (Steinkühler et al., 2018).