Bilayer-thickness-mediated interactions between integral membrane proteins

Hydrophobic thickness mismatch between integral membrane proteins and the surrounding lipid bilayer can produce lipid bilayer thickness deformations. Experiment and theory have shown that protein-induced lipid bilayer thickness deformations can yield energetically favorable bilayer-mediated interactions between integral membrane proteins, and large-scale organization of integral membrane proteins into protein clusters in cell membranes. Within the continuum elasticity theory of membranes, the energy cost of protein-induced bilayer thickness deformations can be captured by considering compression and expansion of the bilayer hydrophobic core, membrane tension, and bilayer bending, resulting in biharmonic equilibrium equations describing the shape of lipid bilayers for a given set of bilayer-protein boundary conditions. Here we develop a combined analytic and numerical methodology for the solution of the equilibrium elastic equations associated with protein-induced lipid bilayer deformations. Our methodology allows accurate prediction of thickness-mediated protein interactions for arbitrary protein symmetries at arbitrary protein separations and relative orientations. We provide exact analytic solutions for cylindrical integral membrane proteins with constant and varying hydrophobic thickness, and develop perturbative analytic solutions for noncylindrical protein shapes. We complement these analytic solutions, and assess their accuracy, by developing both finite element and finite difference numerical solution schemes. We provide error estimates of our numerical solution schemes and systematically assess their convergence properties. Taken together, the work presented here puts into place an analytic and numerical framework which allows calculation of bilayer-mediated elastic interactions between integral membrane proteins for the complicated protein shapes suggested by structural biology and at the small protein separations most relevant for the crowded membrane environments provided by living cells.

Figure 1

Overlapping lipid bilayer thickness deformation fields can yield bilayer-thickness-mediated interactions between membrane proteins. Bilayer thickness deformations $u$ due to interacting membrane proteins obtained using (a) the cylinder model, (b) the crown model, and (c) the clover-leaf model of integral membrane proteins (see Sec. 2 for further details). The thickness deformations in (a) and (b) were obtained through exact analytic minimization of the thickness deformation energy (see Sec. 3), and the thickness deformations in (c) were calculated using finite elements (see Sec. 4a). The clover-leaf shape in (c) provides a simple coarse-grained model [46, 47, 88, 89] of the observed closed-state structure of the pentameric mechanosensitive channel of large conductance [98] (protein structural data shown as ribbon diagrams; Protein Data Bank accession number 2OAR). The calculated lipid bilayer thickness deformations depend on lipid and protein properties, the protein center-to-center distance $d$, and, for the crown and clover-leaf models, the protein orientations ${\omega }_{1,2}$. The color scale ranges from $u_{\text{min}}=-0.6$ nm to $u_{\text{max}}=0.4$ nm.

Figure 2

Two-center bipolar coordinate system for two membrane proteins with circular cross sections of radii $R_{1,2}$ separated by a center-to-center distance $d$ along the $x$ axis. The $x\text{−}y$ plane corresponds to the reference plane used in the Monge representation of surfaces. Expressions for $(r_1,{\theta }_1)$ in terms of $(r_2,{\theta }_2)$ can be obtained by considering the coordinates of the red (left) point, and expressions for $(r_2,{\theta }_2)$ in terms of $(r_1,{\theta }_1)$ can be obtained by considering the coordinates of the blue (right) point.

Figure 3

Thickness deformations $u$ minimizing Eq. (3) obtained using our FE approach for two clover-leaf pentamers. The triangular mesh (black overlay) indicates the mesh used in the FE calculation and is generated using the frontal algorithm of the gmsh package [112]. All model parameters were chosen as described in Sec. 2.

Figure 4

Illustration of the hexagonal grid used for our FD calculations. (a) Nodal values of the discretized thickness deformation field $u_{i,j}$ are indexed by $i$ along the horizontal direction and by $j$ along the oblique direction. The 19-point stencil discretizing the biharmonic operator at position $(i,j)=(0,0)$ is found by multiplying $u_{i,j}$ and nearby nodal values of the thickness deformation field by the indicated coefficients, adding up these contributions, and dividing by $h^4$. (b) We choose the protein boundary points (blue circles) to correspond to the nodes closest to the exact protein boundary curve (gray clover-leaf shape; see Sec. 2b3), and impose the general boundary condition in Eq. (11) at these nodes. We use a layer of interior points (red squares), together with corresponding exterior points (see main text), to impose the general boundary condition in Eq. (12).

Figure 5

Lipid bilayer thickness deformations due to a single cylindrical membrane protein. (a) Bilayer thickness deformation profile $u$ vs radial coordinate $r$ obtained from the exact analytic solution in Eq. (19) for $d\to \infty$ [38], the FE approach, and the FD approach. The gray vertical line indicates the protein boundary. (b) Bilayer thickness deformation energy vs thickness mismatch $U$ obtained using analytic, FE, and FD approaches. The insets show (a) the difference in thickness deformation profile, $\mathrm{\Delta }u$, and (b) the percentage difference in thickness deformation energy, ${\eta }_G$ in Eq. (54), between the analytic solution and the corresponding results of FE (blue solid curves) and FD (green dashed curves) calculations, respectively. For the FE solution we used an average edge size of the FE mesh $\langle l_{\mathrm{edge}}\rangle =0.3$ nm, and for the FD solution we used a lattice spacing $h=0.05$ nm. We set $\tau =0$ for both panels. All model parameters were chosen as described in Sec. 2, and analytic and numerical solutions were obtained as discussed in Secs. 3 and 4.

Figure 6

Comparison of analytic and FE solutions for cylindrical membrane proteins. (a) Percentage error in thickness deformations in Eq. (45) (${\eta }_u$; orange dashed curve) and percentage error in curvature deformations in Eq. (46) (${\eta }_{{\mathbf{\nabla }}^{2}u}$; magenta solid curve) vs $\langle l_{\mathrm{edge}}\rangle$ for a single cylindrical membrane protein. Errors in the thickness mismatch and curvature deformations decay as $\langle l_{\mathrm{edge}}\rangle$ and ${\langle l_{\mathrm{edge}}\rangle }^2$, respectively. (b) Percentage difference between analytic and FE results for the thickness deformation energy per cylindrical membrane protein, ${\eta }_G$ in Eq. (54), vs $\langle l_{\mathrm{edge}}\rangle$ for a system consisting of two identical cylindrical membrane proteins in the noninteracting regime ($d=14$ nm; red solid curve) and in the strongly interacting regime ($d=5$ nm; blue dashed curve; $N=11$ for the analytic solution). We set $\tau =0$ for both panels. All model parameters were chosen as described in Sec. 2, and analytic and numerical solutions were obtained as discussed in Secs. 3 and 4.

Figure 7

Percentage difference between analytic and FD results for the thickness deformation energy per cylindrical membrane protein, ${\eta }_G$ in Eq. (54), vs lattice spacing $h$, for a system consisting of two identical cylindrical membrane proteins in the noninteracting regime ($d=14$ nm for the FD solution; red solid curve) and in the strongly interacting regime ($d=5$ nm; blue dashed curve; $N=11$ for the analytic solution). We set $\tau =0$. All model parameters were chosen as described in Sec. 2, and analytic and numerical solutions were obtained as discussed in Secs. 3 and 4.

Figure 8

Thickness deformation energy per protein, $G$, for two cylindrical membrane proteins vs center-to-center protein distance $d$, calculated analytically at $N=3,N=5,N=7$, and $N=11$ in Eq. (18), and numerically using our FE and FD solution procedures. The inset shows the percentage difference in thickness deformation energy, ${\eta }_G$ in Eq. (54), between the analytic solution at $N=11$ and the corresponding results of FE (blue dashed curve) and FD (orange solid curve) calculations. For the FE solution we used $\langle l_{\mathrm{edge}}\rangle \approx 0.27$ nm, and for the FD solution we used $h=0.05$ nm. We set $\tau =0$. All model parameters were chosen as described in Sec. 2, and analytic and numerical solutions were obtained as discussed in Secs. 3 and 4.

Figure 9

Thickness deformation energy per protein for two crown shapes, $G$, vs center-to-center protein distance $d$, calculated analytically at $N=11$ and numerically using FE and FD schemes, for the minus-minus configuration $[{\omega }_1=0^{\circ }$ and ${\omega }_2={36}^{\circ }$ in Eq. (15)], the plus-minus configuration $[{\omega }_1=0^{\circ }$ and ${\omega }_2=0^{\circ }$ in Eq. (15)], and a configuration intermediate between the minus-minus and plus-minus configurations $[{\omega }_1=0^{\circ }$ and ${\omega }_2={24}^{\circ }$ in Eq. (15)]. We only consider here FD solutions for the mirror-symmetric minus-minus protein configuration. The inset shows the percentage difference in thickness deformation energy, ${\eta }_G$ in Eq. (54), between the analytic result at $N=11$ and the corresponding results of the FE [blue curves; minus-minus (solid), intermediate (short dashed), and plus-minus (long dashed) configurations] and FD (orange triangles with ${\eta }_G\approx 1.2\%$ at $d=5$ nm) calculations. We used $\langle l_{\mathrm{edge}}\rangle \approx 0.27$ nm for the FE and $h=0.05$ nm for the FD solutions, and set $\tau =0$. All model parameters were chosen as described in Sec. 2, and analytic and numerical solutions were obtained as discussed in Secs. 3 and 4.

Figure 10

Comparison of analytic and FE solutions for crown shapes. (a) Percentage error in thickness deformations in Eq. (45) (${\eta }_u$; orange dashed curve) and percentage error in curvature deformations in Eq. (46) (${\eta }_{{\mathbf{\nabla }}^{2}u}$; magenta solid curve) vs $\langle l_{\mathrm{edge}}\rangle$ for two crown shapes in the minus-minus configuration in the interacting regime at $d=7$ nm. The error is evaluated using the corresponding analytic result at $N=11$. Errors in the thickness mismatch and curvature deformations decay as $\langle l_{\mathrm{edge}}\rangle$ and ${\langle l_{\mathrm{edge}}\rangle }^2$, respectively. (b) Percentage difference between analytic and FE results for the thickness deformation energy per protein, ${\eta }_G$ in Eq. (54), vs $\langle l_{\mathrm{edge}}\rangle$ for two crown shapes in the minus-minus configuration in the noninteracting regime ($d=14$ nm for the FE solution; red solid curve) and in the strongly interacting regime ($d=5$ nm; blue dashed curve; $N=11$ for the analytic solution). We set $\tau =1k_{B}T/{\mathrm{nm}}^2$. All model parameters were chosen as described in Sec. 2, and analytic and numerical solutions were obtained as discussed in Secs. 3 and 4.

Figure 11

Percentage difference between analytic and FD results for the thickness deformation energy per protein, ${\eta }_G$ in Eq. (54), vs lattice spacing $h$, for two crown shapes in the plus-plus configuration in the noninteracting regime ($d=14$ nm for the FD solution; red solid curve) and in the strongly interacting regime ($d=5$ nm; blue dashed curve; $N=11$ for the analytic solution). We set $\tau =0$. All model parameters were chosen as described in Sec. 2, and analytic and numerical solutions were obtained as discussed in Secs. 3 and 4.

Figure 12

Thickness deformation energy of clover-leaf shapes. (a) Thickness deformation energy per protein for two clover-leaf shapes, $G$, vs center-to-center protein distance $d$, calculated analytically at $N=11$ and numerically using FE and FD schemes, for the face-on configuration $[{\omega }_1=0^{\circ }$ and ${\omega }_2={36}^{\circ }$ in Eq. (16)] and the tip-on configuration $[{\omega }_1={36}^{\circ }$ and ${\omega }_2=0^{\circ }$ in Eq. (16)]. (b) Thickness-mediated interaction energy $G_{\mathrm{int}}$ vs $d$ obtained by subtracting the protein-induced thickness deformation energies in the noninteracting regime from the respective perturbative analytic, FE, and FD solutions in (a). The inset shows the difference in the thickness-mediated interaction energies obtained from the FE solution procedure, and the analytic and FD solution procedures. We use the same labeling conventions for (b) as for (a). For our numerical calculations, we used $\langle l_{\mathrm{edge}}\rangle \approx 0.26$ nm for the FE and $h=0.05$ nm for the FD solutions. We set $\tau =0$ for both panels. All model parameters were chosen as described in Sec. 2, and analytic and numerical solutions were obtained as discussed in Secs. 3 and 4.