Out-of-equilibrium mechanochemistry and self-organization of fluid membranes interacting with curved proteins

Denoting by ρ_l(θ^α, t) the lipid areal density, balance of mass of lipids requires $\partial {\rho }_{l}/\partial t+{\rho }_{l}\,\mathrm{tr}\,{\boldsymbol{d}}=0$. Note that in this equation $\partial {\rho }_{l}/\partial t$ coincides with the material time-derivative since we consider a Lagrangian parametrization. This equation ignores the area fraction occupied by protein insertions, which except for channels is much smaller than the area fraction occupied by the membrane protein at the periphery of the bilayer, see figure 1. At moderate tensions we can assume that lipids are inextensible, and hence their density constant. Consequently, lipid membrane inextensibility requires that

$$\begin{eqnarray}&&\mathrm{tr}\ {\boldsymbol{d}}=0.\end{eqnarray} \tag{ 4 }$$

Since proteins are a diffusive species, balance of mass of proteins can be expressed as

$$\begin{eqnarray}&&\displaystyle \frac{\partial \phi }{\partial t}+\phi \ \mathrm{tr}\ {\boldsymbol{d}}+{\rm{\nabla }}\cdot (\phi {\boldsymbol{w}})=0\quad \mathrm{on}\quad {\rm{\Gamma }},\end{eqnarray} \tag{ 5 }$$

where ${\boldsymbol{w}}$ is the diffusive velocity relative to the Lagrangian coordinates and we have ignored protein sorption. Again, because of our Lagrangian parametrization, $\partial \phi /\partial t$ coincides with the material time-derivative $\dot{\phi }$.

Zoom In Zoom Out Reset image size

To describe the bending energy of the membrane, we follow a classical Helfrich model [36, 37, 69], ignoring the Gaussian curvature terms for simplicity but including a spontaneous curvature that depends on protein density. Assuming a linear dependence of preferred curvature on protein density, a common form of the bending energy is

$$\begin{eqnarray}&&{{ \mathcal F }}_{b}={\int }_{{\rm{\Gamma }}}\displaystyle \frac{\kappa }{2}{\left(H-\bar{C}\phi \right)}^{2}\,{\rm{d}}{S},\end{eqnarray} \tag{ 6 }$$

where κ is the bending modulus and $\bar{C}$ encodes the curving strength of proteins. Other variants of this model have been proposed [40] and further discussed in section 7.

For the free energy of the proteins on the surface, we adopt a Flory–Huggins [38, 39] model accounting for the entropy of mixing and for protein–protein interactions

$$\begin{eqnarray}&&{{ \mathcal F }}_{p}={\int }_{{\rm{\Gamma }}}\left\{\displaystyle \frac{{k}_{B}T}{{a}_{p}}[\phi {\rm{log}}\phi +({\phi }_{m}-\phi ){\rm{log}}({\phi }_{m}-\phi )]+\displaystyle \frac{\chi }{2{a}_{p}}{\phi }^{2}+\displaystyle \frac{{\mu }_{0}}{{a}_{p}}\phi \right\}\,{\rm{d}}{S}.\end{eqnarray} \tag{ 7 }$$

In this equation, ${k}_{B}T$ is the thermal energy, a_p is the area on the membrane of a single protein so that ϕ/a_p is the number density, the term involving ϕ_m (a saturation area fraction) accounts for the entropy of uncovered spaces, χ determines the strength and sign (attractive or repulsive) of protein–protein interactions, and μ₀ is a reference chemical potential. Variants of this model have been used to understand the linear stability of fully mixed states [42, 43] or to examine protein sorting by curvature at fixed shape [18, 44, 47].

Depending on the parameters and the boundary conditions, the model can lead to phase separation. For instance, negative χ promotes demixing and coexistence of a protein-rich and depleted phases. To regularize phase boundaries, we consider a term penalizing the gradient, which can be interpreted as a non-local interaction of proteins

$$\begin{eqnarray}&&{{ \mathcal F }}_{\mathrm{nl}}={\int }_{{\rm{\Gamma }}}\displaystyle \frac{{\rm{\Lambda }}}{2{a}_{p}}| {\rm{\nabla }}\phi {| }^{2}\,{\rm{d}}{S},\end{eqnarray} \tag{ 8 }$$

Λ being a material parameter. The length scale of this interaction is of the order $\sim \sqrt{{\rm{\Lambda }}/| \chi | }$, which, for a planar membrane, determines the thickness of the interface between protein-rich and depleted domains. Without this term, there is no energetic penalty to domain boundaries and the equations may become ill-posed [70]. Total free energy of this system is then given by

$$\begin{eqnarray}&&{ \mathcal F }={{ \mathcal F }}_{b}+{{ \mathcal F }}_{p}+{{ \mathcal F }}_{\mathrm{nl}}={\int }_{{\rm{\Gamma }}}W(H,\phi ,| {\rm{\nabla }}\phi {| }^{2})\,{\rm{d}}{S},\end{eqnarray} \tag{ 9 }$$

where W represents the total energy density of the membrane-protein composite system. We note that this form of energy density is not the most general that can be conceived for diffusing proteins on an inextensible fluid membrane without orientational order. For instance, the energy could also depend on the scalar invariant ${\rm{\nabla }}\phi \cdot {\boldsymbol{k}}{\rm{\nabla }}\phi$ [71].

A central thermodynamic quantity that drives protein diffusion and sorption is the chemical potential μ [59], measuring the amount of work required to add a molecule at a particular membrane location. The chemical potential can be identified from the variational derivative of ${ \mathcal F }$ with respect to ϕ as

$$\begin{eqnarray}&&{\delta }_{\phi }{ \mathcal F }={\int }_{{\rm{\Gamma }}}\mu \displaystyle \frac{\delta \phi }{{a}_{p}}\,{\rm{d}}{S}={\int }_{{\rm{\Gamma }}}\left[{D}_{2}W\ \delta \phi +{D}_{3}W\ \delta (| {\rm{\nabla }}\phi {| }^{2})\right]\,{\rm{d}}{S},\end{eqnarray} \tag{ 10 }$$

where ${D}_{i}W$ denotes the partial derivative of W with respect to its ith argument. Because this variation is taken at fixed shape, $\delta (| {\rm{\nabla }}\phi {| }^{2})=2{\rm{\nabla }}\phi \cdot {\rm{\nabla }}\delta \phi$. Using the notation ${W}_{H}={D}_{1}W,{W}_{\phi }={D}_{2}W$, and ${W}_{{\rm{\nabla }}\phi }=2{D}_{3}W\ {\rm{\nabla }}\phi$, and integrating by parts, we obtain

$$\begin{eqnarray}&&{\delta }_{\phi }{ \mathcal F }={\int }_{{\rm{\Gamma }}}({W}_{\phi }-{\rm{\nabla }}\cdot {W}_{{\rm{\nabla }}\phi })\delta \phi \,{\rm{d}}{S}+{\int }_{\partial {\rm{\Gamma }}}(\delta \phi ){W}_{{\rm{\nabla }}\phi }\cdot {\boldsymbol{\nu }}\,{\rm{d}}{l}\end{eqnarray} \tag{ 11 }$$

with ${\boldsymbol{\nu }}$ representing the in-plane normal at the edge ∂Γ. Ignoring the boundary term, comparison of the two equations above allows us to identify the chemical potential from equation (9) as,

$$\begin{eqnarray}&&\mu ={a}_{p}({W}_{\phi }-{\rm{\nabla }}\cdot {W}_{{\rm{\nabla }}\phi })={\mu }_{0}-{a}_{p}\kappa \bar{C}(H-\bar{C}\phi )+{k}_{B}T\mathrm{log}\displaystyle \frac{\phi }{{\phi }_{m}-\phi }+\chi \phi -{\rm{\Lambda }}{\rm{\Delta }}\phi ,\end{eqnarray} \tag{ 12 }$$

where Δ denotes the surface Laplacian. This expression clearly shows the contributions to the chemical potential given by bending elasticity, entropy and protein self-interactions.

Having described the mechanisms through which membrane and proteins store energy, we now detail the dynamical modes through which the system dissipates energy. Lipid bilayers in a fluid phase behave like interfacial viscous Newtonian fluids [67, 72]. Here, we ignore dissipative forces resulting from inter-monolayer slippage [52, 53, 73]. Having assumed that the membrane is locally inextensible, $\mathrm{tr}\ {\boldsymbol{d}}=0$, the dissipation potential accounting for in-plane shear can be written as

$$\begin{eqnarray}&&{{ \mathcal D }}_{m}={\int }_{{\rm{\Gamma }}}\eta \ {\boldsymbol{d}}:{\boldsymbol{d}}\,{\rm{d}}{S},\end{eqnarray} \tag{ 13 }$$

where η is the in-plane shear viscosity of the membrane and the rate-of-deformation tensor ${\boldsymbol{d}}$ is given by equation (2) [67].

Protein transport is characterized by the collective protein velocity ${\boldsymbol{w}}$ relative to the lipids in the membrane, which has already appeared in the statement of balance of proteins, see equation (5). As a single protein moves relative to the lipids, it experiences a drag force given by $-\xi {\boldsymbol{w}}$, where ξ is the molecular drag coefficient. By superposition in a dilute approximation, a collection of proteins characterized by local number density ϕ/a_p experiences a drag force per unit area given by $-\xi \phi {\boldsymbol{w}}/{a}_{p}$. The associated dissipation potential can then be written as

$$\begin{eqnarray}&&{{ \mathcal D }}_{p}={\int }_{{\rm{\Gamma }}}\displaystyle \frac{\xi }{2{a}_{p}}\phi | {\boldsymbol{w}}{| }^{2}\,{\rm{d}}{S}.\end{eqnarray} \tag{ 14 }$$

For simplicity, we ignore the dissipation occurring in the bulk fluid. This approximation could break down if fast shape changes occurred over length-scales larger than the Saffman–Dellbrück length, of about 1–10 µm [67, 74]. In most situations, however, fast curvature generation by proteins leads to much smaller geometric features [19, 75, 76]. Thus, the total dissipation potential of the system is given by:

$$\begin{eqnarray}&&{ \mathcal D }[{\boldsymbol{v}},{v}_{n},{\boldsymbol{w}}]={{ \mathcal D }}_{m}+{{ \mathcal D }}_{p}.\end{eqnarray} \tag{ 15 }$$

Let us consider a membrane patch Γ, possibly with smooth boundary ∂Γ. In absence of body forces or sorption of proteins from the bulk, power can only be supplied to the system through edge tractions, moments or flux of proteins at a given chemical potential. Assuming that Γ is a material surface, and thus no lipids can flow through ∂Γ, we write the power input functional as

$$\begin{eqnarray}&&{{ \mathcal P }}_{\mathrm{ext}}={\int }_{\partial {\rm{\Gamma }}}\displaystyle \frac{{\mu }_{b}}{{a}_{p}}(\phi {\boldsymbol{w}})\cdot {\boldsymbol{\nu }}\,{\rm{d}}{l}-{\int }_{\partial {\rm{\Gamma }}}({F}_{\tau }{\boldsymbol{v}}\cdot {\boldsymbol{\tau }}+{F}_{\nu }{\boldsymbol{v}}\cdot {\boldsymbol{\nu }}+{F}_{n}{v}_{n})\,{\rm{d}}{l}+{\int }_{\partial {\rm{\Gamma }}}M{\boldsymbol{\nu }}\cdot \dot{{\boldsymbol{n}}}\,{\rm{d}}{l},\end{eqnarray} \tag{ 16 }$$

where μ_b is a fixed chemical potential for proteins at the boundary, e.g. maintained by a protein reservoir, ${\boldsymbol{\tau }}$ is a unit tangent vector along ∂Γ so that ${\boldsymbol{\nu }}={\boldsymbol{\tau }}\times {\boldsymbol{n}},{F}_{\tau }$, F_ν and F_n are traction components at the boundary, M is a bending moment per unit length, and $\dot{{\boldsymbol{n}}}$ represents the material time derivative of the surface normal. We can express the last integral in terms of our dynamical variables ${\boldsymbol{v}}$ and v_n using the relation

$$\begin{eqnarray}&&{\boldsymbol{\nu }}\cdot \dot{{\boldsymbol{n}}}=-{\rm{\nabla }}{v}_{n}\cdot {\boldsymbol{\nu }}-\tau ({\boldsymbol{v}}\cdot {\boldsymbol{\tau }})-{\kappa }_{\nu }({\boldsymbol{v}}\cdot {\boldsymbol{\nu }}),\end{eqnarray} \tag{ 17 }$$

where $\tau ={\boldsymbol{k}}{\boldsymbol{\nu }}\cdot {\boldsymbol{\tau }}$ is the geodesic torsion of the boundary curve and ${\kappa }_{\nu }={\boldsymbol{k}}{\boldsymbol{\nu }}\cdot {\boldsymbol{\nu }}$ its normal curvature [77, 78]. Each of these terms can be defined on Neumann parts of the boundary ∂Γ.

Onsager's recipe to derive the dissipative dynamics is to minimize the Rayleghian functional defined as

$$\begin{eqnarray}&&{ \mathcal R }=\dot{{ \mathcal F }}+{ \mathcal D }+{{ \mathcal P }}_{\mathrm{ext}},\end{eqnarray} \tag{ 18 }$$

where $\dot{{ \mathcal F }}$ is the rate of change of the energy [58, 59]. To enforce local inextensibility, see equation (4), and possibly the incompressibility of the fluid enclosed by the membrane, we introduce a Lagrange multiplier field σ (contributing to membrane tension) and a Lagrange multiplier p (pressure difference) to form the Lagrangian

$$\begin{eqnarray}&&{ \mathcal L }=\dot{{ \mathcal F }}+{ \mathcal D }+{{ \mathcal P }}_{\mathrm{ext}}+{\int }_{{\rm{\Gamma }}}\sigma \ \mathrm{tr}\ {\boldsymbol{d}}\,{\rm{d}}{S}-p{\int }_{{\rm{\Gamma }}}{v}_{n}\,{\rm{d}}{S}.\end{eqnarray} \tag{ 19 }$$

Since Γ is a material surface, the rate of change of energy can be written as

$$\begin{eqnarray}&&\dot{{ \mathcal F }}={\int }_{{\rm{\Gamma }}}(\dot{W}+W\mathrm{tr}{\boldsymbol{d}})\,{\rm{d}}{S},\end{eqnarray} \tag{ 20 }$$

where $\dot{W}$ is the material time-derivative of the energy density and is given by

$$\begin{eqnarray}&&\dot{W}={W}_{H}\dot{H}+{W}_{\phi }\dot{\phi }+{D}_{3}W\ \displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\left(| {\rm{\nabla }}\phi {| }^{2}\right).\end{eqnarray} \tag{ 21 }$$

Note carefully that, because $| {\rm{\nabla }}\phi {| }^{2}={g}^{\alpha \beta }{\phi }_{,\alpha }{\phi }_{,\beta }$ with ${\phi }_{,\alpha }=\partial \phi /\partial {\theta }^{\alpha }$ involves the inverse of the metric tensor, which depends on the membrane configuration, the time derivative in the last term involves not only $\dot{\phi }$ but also the velocity of the surface. Indeed,

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\left({g}^{\alpha \beta }{\phi }_{,\alpha }{\phi }_{,\beta }\right)=-2{\rm{\nabla }}\phi \cdot {\boldsymbol{d}}\cdot {\rm{\nabla }}\phi +2{\rm{\nabla }}\phi \cdot {\rm{\nabla }}\dot{\phi }.\end{eqnarray} \tag{ 22 }$$

Using the above equation, the relation $\dot{H}={\boldsymbol{v}}\cdot {\rm{\nabla }}H+{\rm{\Delta }}{v}_{n}+{v}_{n}({H}^{2}-2K)$ [67], balance of mass of proteins in equation (5) and the divergence theorem, we can express the rate of change of the energy by explicitly highlighting its dependence on the variables ${\boldsymbol{w}}$ and ${\boldsymbol{V}}={\boldsymbol{v}}+{v}_{n}{\boldsymbol{n}}$

$$\begin{eqnarray}\begin{array}{rcl}\dot{{ \mathcal F }} & = & {\displaystyle \int }_{{\rm{\Gamma }}}\displaystyle \frac{{\rm{\nabla }}\mu }{{a}_{p}}\cdot (\phi {\boldsymbol{w}})\,{\rm{d}}{S}\\ & & +{\displaystyle \int }_{{\rm{\Gamma }}}\left[{W}_{H}{\rm{\nabla }}H+{\rm{\nabla }}\left(\displaystyle \frac{\mu \phi }{{a}_{p}}\right)-{\rm{\nabla }}W+2{D}_{3}W{\rm{\nabla }}\cdot ({\rm{\nabla }}\phi \otimes {\rm{\nabla }}\phi )\right]\cdot {\boldsymbol{v}}\,{\rm{d}}{S}\\ & & +{\displaystyle \int }_{{\rm{\Gamma }}}\left[{\rm{\Delta }}{W}_{H}+{W}_{H}({H}^{2}-2K)+\displaystyle \frac{\mu }{{a}_{p}}\phi H-{WH}+2{D}_{3}W{\boldsymbol{k}}:({\rm{\nabla }}\phi \otimes {\rm{\nabla }}\phi )\right]{v}_{n}\,{\rm{d}}{S}\\ & & -{\displaystyle \int }_{\partial {\rm{\Gamma }}}\displaystyle \frac{\mu }{{a}_{p}}(\phi {\boldsymbol{w}})\cdot {\boldsymbol{\nu }}\,{\rm{d}}{l}-{\displaystyle \int }_{\partial {\rm{\Gamma }}}2{D}_{3}W({\rm{\nabla }}\phi \cdot {\boldsymbol{\nu }}){\rm{\nabla }}\cdot (\phi {\boldsymbol{w}})\,{\rm{d}}{l}+{\displaystyle \int }_{\partial {\rm{\Gamma }}}\left(W-\displaystyle \frac{\mu \phi }{{a}_{p}}\right){\boldsymbol{v}}\cdot {\boldsymbol{\nu }}\,{\rm{d}}{l}\\ & & -{\displaystyle \int }_{\partial {\rm{\Gamma }}}({\rm{\nabla }}{W}_{H}\cdot {\boldsymbol{\nu }}){v}_{n}\,{\rm{d}}{l}+{\displaystyle \int }_{\partial {\rm{\Gamma }}}{W}_{H}({\rm{\nabla }}{v}_{n}\cdot {\boldsymbol{\nu }})\,{\rm{d}}{l}-{\displaystyle \int }_{\partial {\rm{\Gamma }}}2{D}_{3}W({\rm{\nabla }}\phi \cdot {\boldsymbol{v}})({\rm{\nabla }}\phi \cdot {\boldsymbol{\nu }})\,{\rm{d}}{l}.\end{array}\end{eqnarray} \tag{ 23 }$$

In the literature dealing with the Cahn–Hilliard equation, related to our model, the second term in the fourth line is set to zero by requiring the natural high-order boundary condition ${\rm{\nabla }}\phi \cdot {\boldsymbol{\nu }}=0$ on ∂Γ [79]. As a result, the last integral over the edge also vanishes. Similarly, we can express the constraint of local area conservation as

$$\begin{eqnarray}&&{\int }_{{\rm{\Gamma }}}\sigma \ \mathrm{tr}\ {\boldsymbol{d}}\,{\rm{d}}{S}=-{\int }_{{\rm{\Gamma }}}{\rm{\nabla }}\sigma \cdot {\boldsymbol{v}}\,{\rm{d}}{S}-{\int }_{{\rm{\Gamma }}}\sigma {v}_{n}H\,{\rm{d}}{S}+{\int }_{\partial {\rm{\Gamma }}}\sigma {\boldsymbol{v}}\cdot {\boldsymbol{\nu }}\,{\rm{d}}{l}.\end{eqnarray} \tag{ 24 }$$

To obtain the governing equations of the system, we substitute the above relations in equation (19), minimize the Lagrangian with respect to $\{{\boldsymbol{w}},{\boldsymbol{v}},{v}_{n}\}$ and maximize it with respect to {σ, p}. Surface integrals in the two expressions above contribute to the Euler–Lagrange governing equations whereas boundary terms identify the Neumann boundary conditions. The resulting equations for the specific choice of W given in section 3.1 are discussed below.

Minimizing the Lagrangian with respect to the diffusive velocity, ${\delta }_{{\boldsymbol{w}}}{ \mathcal L }=0$, we obtain the diffusive flux of membrane proteins consistent with Fick's law

$$\begin{eqnarray}&&\phi {\boldsymbol{w}}=-\displaystyle \frac{\phi {\rm{\nabla }}\mu }{\xi }\ \ \mathrm{in}\ \ {\rm{\Gamma }}.\end{eqnarray} \tag{ 25 }$$

Substituting the above relation in equation (5) and using the inextensibility of the lipid membrane we obtain

$$\begin{eqnarray}&&\displaystyle \frac{\xi }{{a}_{p}}\dot{\phi }-{\rm{\nabla }}\cdot (\phi {\rm{\nabla }}\mu )=0,\end{eqnarray} \tag{ 26 }$$

which can be expanded recalling equation (12) and rearranged to yield the transport equation for proteins

$$\begin{eqnarray}&&\xi \dot{\phi }+{\rm{\nabla }}\cdot \left[{a}_{p}\kappa \bar{C}{\rm{\nabla }}H\phi -\left({\chi }^{\mathrm{eff}}\phi +\displaystyle \frac{{k}_{B}T{\phi }_{m}}{{\phi }_{m}-\phi }\right){\rm{\nabla }}\phi +{\rm{\Lambda }}\phi {\rm{\nabla }}({\rm{\Delta }}\phi )\right]=0,\end{eqnarray} \tag{ 27 }$$

where we have defined the effective interaction between proteins as ${\chi }^{\mathrm{eff}}=\chi +{a}_{p}\kappa {\bar{C}}^{2}$. For vanishing spontaneous curvature ($\bar{C}=0$), this governing equation ceases to depend explicitly on the curvature of the underlying surface, although it does depend on its geometry through the covariant derivative, and it reduces to a nonlinear Cahn–Hilliard equation [80] on the surface with a density-dependent diffusion coefficient

$$\begin{eqnarray}&&D(\phi )=\displaystyle \frac{1}{\xi }\left({\chi }^{\mathrm{eff}}\phi +\displaystyle \frac{{k}_{B}T{\phi }_{m}}{{\phi }_{m}-\phi }\right).\end{eqnarray} \tag{ 28 }$$

For χ = 0, Λ = 0 and ϕ ≪ ϕ_m, equation (27) reduces to the linear diffusion equation. For $\bar{C}\ne 0$, the curvature gradient introduces a bias in the diffusion with drift velocity

$$\begin{eqnarray}&&{{\boldsymbol{w}}}^{\mathrm{drift}}={a}_{p}\kappa \bar{C}{\rm{\nabla }}H/\xi ,\end{eqnarray} \tag{ 29 }$$

driving curved proteins along or against the curvature gradient depending on the sign of $\bar{C}$. At steady state, drift and diffusive transport must balance and yield a divergence-free flux.

Minimization of the Lagrangian respect to the in-plane velocity, ${\delta }_{{\boldsymbol{v}}}{ \mathcal L }=0$, yields

$$\begin{eqnarray}&&{\rm{\nabla }}\sigma -\displaystyle \frac{\phi {\rm{\nabla }}\mu }{{a}_{p}}+2\eta {\rm{\nabla }}\cdot {\boldsymbol{d}}={\bf{0}},\end{eqnarray} \tag{ 30 }$$

where the first term accounts for the tension required to impose lipid membrane inextensibility, the second term is a force density on the fluid membrane resulting from the relative motion of proteins, see equation (25), and the third term represents tangential dissipative forces due to membrane viscosity, which strongly depend on curvature and involve both tangential and normal velocities as discussed in detail elsewhere [67, 81].

In the common situation of an incompressible Stokes flow with a dilute diffusing species, the drag force density due to protein motion, the second term in equation (30), does not contribute to the hydrodynamics because ϕ∇μ can be expressed as a gradient and grouped with the Lagrange multiplier enforcing incompressibility in all the governing equations [59]. Here, however, this is not the case. Indeed, introducing equation (12) into (30) we obtain

$$\begin{eqnarray}&&{\rm{\nabla }}{\sigma }^{\mathrm{eff}}-\displaystyle \frac{{\rm{\Lambda }}}{{a}_{p}}{\rm{\nabla }}\cdot {\left({\rm{\nabla }}\phi \otimes {\rm{\nabla }}\phi \right)}^{\mathrm{dev}}+\displaystyle \frac{\xi }{{a}_{p}}\phi {{\boldsymbol{w}}}^{\mathrm{drift}}+2\eta {\rm{\nabla }}\cdot {\boldsymbol{d}}={\bf{0}},\end{eqnarray} \tag{ 31 }$$

where ${({\boldsymbol{a}}\otimes {\boldsymbol{a}})}^{\mathrm{dev}}={\boldsymbol{a}}\otimes {\boldsymbol{a}}-(| {\boldsymbol{a}}{| }^{2}/2){\boldsymbol{g}}$ is the deviatoric component of this rank-one tensor, and we identify the effective membrane tension (the hydrostatic part of the stress tensor) as

$$\begin{eqnarray}&&{\sigma }^{\mathrm{eff}}=\sigma +\displaystyle \frac{{k}_{B}T}{{a}_{p}}{\phi }_{m}\mathrm{log}({\phi }_{m}-\phi )-\displaystyle \frac{{\chi }^{\mathrm{eff}}}{2{a}_{p}}{\phi }^{2}+\displaystyle \frac{{\rm{\Lambda }}}{{a}_{p}}\phi {\rm{\Delta }}\phi .\end{eqnarray} \tag{ 32 }$$

Since the third term in equation (31) cannot be expressed as a gradient, the drift contribution to protein transport generates a tangential force density introducing an explicit coupling between hydrodynamics and protein transport. This protein-induced force just requires the presence of curved proteins and a curvature gradient and can drive flows out-of-equilibrium. Furthermore, as a result of the fundamental in-plane/out-of-plane coupling mediated by curvature [81], this force density can induce out-of-plane forces and shape changes. At steady state, since ${{\boldsymbol{w}}}^{\mathrm{drift}}$ is in general different from zero, equation (31) shows that we can expect a non-uniform effective membrane tension in the presence of gradients of curvature. The second term in equation (31) further contributes to a non-uniform and also non-hydrostatic membrane tension in equilibrium, in this case associated to gradients in ϕ.

Going back to the notion of effective tension σ^eff, one way to think about it is just as a Lagrange multiplier enforcing local inextensibility. However, equation (32) provides further insight about this tension. The first term, σ can be interpreted as the membrane tension of the lipids. The second term is the osmotic tension due to the presence of proteins. In the limit ϕ ≪ ϕ_m, this term becomes $-{k}_{B}T\phi /{a}_{p}$, recovering the Van't Hoff's equation. The third term accounts for the fact that attractive/repulsive proteins increase/decrease surface tension. The last term is a non-local tension that can become significant at phase boundaries. Thus, σ^eff is a measure of tension in the composite membrane-protein system. We refer to appendix A for a complementary and detailed derivation of the stress in the present theory.

Minimization with respect to the normal component of the velocity, ${\delta }_{{v}_{n}}{ \mathcal L }=0$, yields

$$\begin{eqnarray}&&\begin{array}{l}\kappa \left[{\rm{\Delta }}(H-\bar{C}\phi )+(H-\bar{C}\phi )\left(\displaystyle \frac{{H}^{2}}{2}-2K\right)-\displaystyle \frac{{H}^{2}\bar{C}\phi }{2}\right]\\ \quad -\,{\sigma }^{\mathrm{eff}}H-p+\displaystyle \frac{{\rm{\Lambda }}}{{a}_{p}}{\boldsymbol{k}}:{\left({\rm{\nabla }}\phi \otimes {\rm{\nabla }}\phi \right)}^{\mathrm{dev}}-2\eta {\boldsymbol{d}}:{\boldsymbol{k}}=0,\end{array}\end{eqnarray} \tag{ 33 }$$

where the first term corresponds to an out-of-plane force density due to membrane curvature elasticity, and hence coupled to the protein density, the second to fourth terms account for Laplace's law, and the last term is a normal viscous force density due to membrane shear, studied in detail in [67, 81].

In the chemo-mechanical problem studied here, we can impose Dirichlet boundary conditions in part of the boundary ∂Γ for the different fields. For instance, reasonable Dirichlet conditions for the mechanical part of the problem include fixing ${\boldsymbol{v}}\cdot {\boldsymbol{\tau }}={\hat{v}}_{\tau },{\boldsymbol{v}}\cdot {\boldsymbol{\nu }}={\hat{v}}_{\nu },{v}_{n}={\hat{v}}_{n}$ and ${\rm{\nabla }}{v}_{n}\cdot {\boldsymbol{\nu }}=\hat{\omega }$ on a part of ∂ Γ, where ${\hat{v}}_{\tau },{\hat{v}}_{\nu },{\hat{v}}_{n}$ and $\hat{\omega }$ are Dirichlet data. For the chemical problem, the Dirichlet boundary condition is prescribing $\phi {\boldsymbol{w}}\cdot {\boldsymbol{\nu }}={\hat{j}}_{b}$ on part of ∂Γ.

On parts of the boundary where Dirichlet boundary conditions are not specified, we obtain the following Neumann boundary conditions by extremizing the Lagrangian and collecting terms at the boundary:

$$\begin{eqnarray}\begin{array}{rcl}{\mu }_{b} & = & \mu ,\ \ \ \ \ \ M={W}_{H},\ \ \ \ \ \ {F}_{\nu }=W-\displaystyle \frac{\mu \phi }{{a}_{p}}+\sigma -{\kappa }_{\nu }M+2\eta \ {\boldsymbol{\nu }}\cdot {\boldsymbol{d}}\cdot {\boldsymbol{\nu }},\\ {F}_{n} & = & -{\rm{\nabla }}{W}_{H}\cdot {\boldsymbol{\nu }},\ \ \ \ \ \ {F}_{\tau }=-\tau M+2\eta \ {\boldsymbol{\tau }}\cdot {\boldsymbol{d}}\cdot {\boldsymbol{\nu }}.\end{array}\end{eqnarray} \tag{ 34 }$$

The first condition sets the chemical potential in the Neumann part of the boundary, whereas the next four equations set the applied torque and force per unit length. The in-plane force can be further recast in a form clearly showing the contributions of the effective surface tension and bending energy as

$$\begin{eqnarray}&&{F}_{\nu }=\displaystyle \frac{\kappa }{2}{H}^{2}+{\sigma }^{\mathrm{eff}}-{\kappa }_{\nu }M+2\eta \ {\boldsymbol{\nu }}\cdot {\boldsymbol{d}}\cdot {\boldsymbol{\nu }}.\end{eqnarray} \tag{ 35 }$$

For surfaces at equilibrium with flat boundaries, the forces normal to the boundary per unit length transmitted by the membrane is tangential and given by σ^eff.

The chemo-mechanical model described by the transport equation for proteins (27), membrane inextensibility (4), force balance in-plane (31) and out-of-plane (33) exhibits multiple intrinsic relaxation time-scales. To examine the competition of mechanical and chemical relaxation time-scales, we consider the simplest cases of shape disturbances whose relaxation is driven by bending elasticity and dragged by membrane viscosity, and of protein density disturbances entropically penalized and dragged by the friction between proteins and the lipid membrane. Such mechanical relaxation occurs during the characteristic time τ_m = ηℓ²/κ , where ℓis the characteristic size of the disturbance. The characteristic time for protein diffusion is ${\tau }_{p}=\xi {{\ell }}^{2}/({k}_{B}T)$. Their ratio is thus

$$\begin{eqnarray}&&\displaystyle \frac{{\tau }_{m}}{{\tau }_{p}}=\displaystyle \frac{{k}_{B}T\eta }{\kappa \xi }\approx \displaystyle \frac{\eta }{20\xi },\end{eqnarray} \tag{ 36 }$$

where we have used the common estimate for the bending stiffness $\kappa \approx 20\,{k}_{B}T$. In turn, ξ is related to the membrane surface viscosity η through the Saffman–Delbrück theory [82], which states that

$$\begin{eqnarray}&&\xi =\displaystyle \frac{4\pi \eta }{\mathrm{log}(2{L}_{\mathrm{sd}}/{{\ell }}_{a})-\gamma },\end{eqnarray} \tag{ 37 }$$

where L_sd ≈ 5 μm is the ratio of membrane and bulk viscosity, ℓ_a ≈ 1 nm is the effective radius of the protein and γ ≈ 0.577 is the Euler Mascheroni constant. Using this estimation we obtain, ξ ≈ 2πη, which results in ${\tau }_{m}/{\tau }_{p}\approx 1/(40\pi )$. Thus, mechanical relaxation occurs nearly two orders of magnitude faster than protein relaxation by diffusion. Although other phenomena accounted for in our general theory (such as protein self-interaction, drift by curvature gradients, or mechanical forces due to the presence of curved proteins on the membrane) can influence the dynamics of the system, this simple estimate establishes that in general protein transport can be expected to be the slow process.

We choose as the energy scale the bending rigidity of the membrane $\kappa =20\,{k}_{B}T$. As the length and time-scales, we choose ℓ₀ = 50 nm and ${\tau }_{p}={{\ell }}_{0}^{2}\xi /({k}_{B}T)$. Considering a membrane viscosity of $\eta =5\times {10}^{-9}\,{\rm{N}}\,{\rm{s}}\,{{\rm{m}}}^{-1}$ and the relation ξ = 2πη discussed in section 4.6, we obtain that $\xi \approx 3\times {10}^{-8}\,{\rm{N}}\,{\rm{s}}\,{{\rm{m}}}^{-1}$, the diffusion coefficient of proteins is $D={k}_{B}T/\xi \approx 1.3\times {10}^{-13}\,{{\rm{m}}}^{2}\,{{\rm{s}}}^{-1}$ and τ_p ≈ 0.02 s. In the absence of measurements, we choose ${\rm{\Lambda }}/{a}_{p}=1{k}_{B}T$ large enough so that, when phase separation occurs, domain boundaries have a finite thickness and simulations are devoid of numerical oscillations indicative of ill-conditioning, and small enough so that the dynamics of the problem are not significantly affected by this parameter. With these units, in our calculations we set the non-dimensional coefficients $\bar{\kappa }=1,\bar{\bar{C}}=2$ (corresponding to $1/\bar{C}\approx 25\,\mathrm{nm}$), ${\bar{a}}_{p}=0.04$ (corresponding to a_p ≈ 100 nm²), $\overline{{k}_{B}T}/{\bar{a}}_{p}=1.25$ , $\bar{\xi }/{\bar{a}}_{p}=1.25$ and $\bar{\eta }=1/(40\pi )$. We finally note that, to avoid numerical solutions with unreasonably thin necks, thinner than the bilayer thickness, we introduce a term that limits the minimum radius of a neck structure to about ∼ 7.5 nm by adding an energy contribution of the form

$$\begin{eqnarray}&&{{ \mathcal F }}_{\mathrm{neck}}={\int }_{{{\rm{\Gamma }}}_{s}}\displaystyle \frac{\gamma }{{\left(\rho -\epsilon \right)}^{2}}\,{\rm{d}}{S},\end{eqnarray} \tag{ 48 }$$

where Γ_s is the entire surface excluding a small region near the poles and γ = 0.1k_BT. We checked that this potential only affected the solutions close to the neck.

Curvature sensing is a phenomenon by which curved membrane proteins migrate to regions of the membrane with higher/preferred curvature. Hence, a necessary condition is the existence of a curvature gradient. We first considered a prolate vesicle as shown in figure 2(a). This vesicle is obtained by minimization of bending energy at constant area $S=4\pi {R}_{0}^{2}$, with R₀ = 500 nm, at fixed reduced volume v = 0.93. At t = 0, the vesicle is covered with a homogeneous area fraction of curved proteins ($\bar{\phi }=0.15$) with spontaneous curvature $\bar{C}=1/25\,{\mathrm{nm}}^{-1}$. We assume that the proteins are non-interacting and thus choose $\chi =-{a}_{p}\kappa {\bar{C}}^{2}$ so that χ^eff = 0. The initial homogeneous distribution of proteins is preferred entropically, but is not optimal from the point of view of bending energy, which favors protein migration towards the poles. The competition between these two free energy contributions leads to a non-uniform chemical potential of proteins and drives protein transport. Since ${\rm{\nabla }}\phi ={\bf{0}}$ at t = 0, protein transport is initially due exclusively to gradients in curvature with the diffusive velocity coinciding with the drift velocity ${{\boldsymbol{w}}}^{\mathrm{drift}}={a}_{p}\kappa \bar{C}{\rm{\nabla }}H/\xi$. Estimating the average gradient of mean curvature from the prolate shape, we estimate the time required for drift transport to induce a gradient in protein density as ${\tau }_{1}\approx {R}_{0}/| {{\boldsymbol{w}}}^{\mathrm{drift}}| \approx 0.3\,{\rm{s}}$. Subsequently, the dynamics are governed by a competition between drift and diffusive transport, driving proteins towards equilibrium over a time scale ${\tau }_{2}\approx {R}_{0}^{2}/D(\bar{\phi })\approx 1.36\,{\rm{s}}$. Thus, we estimate that the total time scale of relaxation is given by τ ≈ τ₁ + τ₂ ≈ 1.66 s. These estimates are consistent with the results shown in figure 2(a), where we show a few selected snapshots of protein distribution during the equilibration dynamics, along with the time-evolution of the changes in the total energy, ${\rm{\Delta }}{ \mathcal F }$, and of different components of it. The figure shows that equilibration takes place in a time commensurate to τ, and that protein migration towards the poles is driven by bending energy, which decreases during the dynamics, but opposed by ${{ \mathcal F }}_{p}+{{ \mathcal F }}_{{nl}}$. In this example, where the total number of proteins on the membrane is low, their area fraction in the protein-rich poles is far from the saturation area fraction ϕ_m = 1 and membrane shape does not change.

Zoom In Zoom Out Reset image size

To examine the ability of proteins to generate membrane shape, we revisited the previous example but increased the amount of protein by setting an initial homogeneous area fraction of $\bar{\phi }=0.35$, see figure 2(b). At early times, the dynamics parallel those of the previous example, with drift motion of proteins towards the poles, followed by balancing diffusive fluxes. However, now the amount of protein creates a sufficiently large spontaneous curvature $\bar{C}\phi$ to modify the shape of the vesicle, which develops a symmetry-breaking transition (ii). At this point, a positive feedback is established during which higher curvature attracts more proteins, which in turn locally increase curvature, and so on. We observe that during this process, the systems develops a cascade of rapid pearling events of duration τ_m, which create new curvature gradients and thus are followed by the partial equilibration of protein coverage over a time-scale of τ_p. Similar pearled tubular morphologies have been experimentally and computationally observed in bilayers with isotropic spontaneous curvature caused by anchored polymers [11, 86], in cells as a result of crowding of the grycocalyx [12] or by asymmetric lipid swelling due to changes in pH [87]. We note that if proteins induce anisotropic spontaneous curvature, for instance because they are elongated and adopt nematic order, experiments and molecular models suggest that one can expect tubular protein-rich protrusions of uniform radius [15, 28, 30] rather than pearled protrusions as we find here. We leave models capturing nematic ordering of curved proteins for future work. This pearling cascade and positive feedback loop between curvature and protein coverage continues until proteins almost reach their saturation density ϕ_m in the highly curved domain. In equilibrium, the system reaches a heterogeneous state where a protein-rich and highly curved pearled tube coexists with a depleted vesicle.

To further examine the role of the saturation area fraction ϕ_m in setting the equilibrium state, we repeated the previous simulation considering different values of ϕ_m. Figure 2(c) shows that the depth of the pearling cascade is indeed controlled by this parameter. For ϕ_m = 0.75, the saturation density and equilibrium is reached after the first pearl has formed. As the saturation area coverage is increased, the number of pearls, and the tube length and curvature progressively increase whereas for lower values of ϕ_m the system does not even pearl. Thus, in a model governed by the energies in equations (6) and (7), protein saturation controlled by ϕ_m limits the positive feedback loop between curvature and area coverage.

Curvature generation by membrane proteins involves recruitment of membrane area into protein-rich protrusions, and therefore should depend on membrane tension as shown experimentally [19]. To examine this mechanical coupling, we started from an equilibrium state showing a highly curved protein-rich tube and increased the pressure difference in steps, thus allowing for volume changes in the vesicle. Initially, p₀ = 13 Pa and σ^eff ≈ 0.003 mN m⁻¹. Figure 2(d) shows that as pressure, and thus tension, increase, membrane area is released from the protein-rich tube, which becomes shorter and more concentrated (p₁ = 55 Pa, p₂ = 65 Pa). Beyond p₃ = 250 Pa corresponding to σ^eff ≈ 0.064 mN m⁻¹, the entire protrusion is eliminated and the proteins uniformly spread over the membrane. This example thus shows the mechanically-induced dissolution of a protein-rich curved domain.

The transition between states of low curvature and homogeneous protein distribution and localized states has been classically analyzed assessing the linear stability of the uniform state [19, 42], summarized in appendix B. According to this analysis, a purely mechanical instability (Euler buckling) takes place when σ^eff < 0, while a purely chemical instability (phase separation) takes place when $D(\bar{\phi })\lt 0$, see equation (28). In addition to these standard instabilities, the system can also exhibit a chemo-mechanical instability involving shape and protein patterning, which in the ideal case of a planar infinite membrane can happen when (appendix B) [19]

$$\begin{eqnarray}&&\kappa | \bar{C}| -{\left[\displaystyle \frac{\kappa }{{a}_{p}}\mathop{\underbrace{\left({\chi }^{\mathrm{eff}}+\displaystyle \frac{{k}_{B}T{\phi }_{m}}{({\phi }_{m}-\bar{\phi })\bar{\phi }}\right)}}\limits_{{a}^{\mathrm{eff}}}\right]}^{\tfrac{1}{2}}\geqslant \sqrt{\displaystyle \frac{{\sigma }^{\mathrm{eff}}{\rm{\Lambda }}}{{a}_{p}}}.\end{eqnarray} \tag{ 49 }$$

This equation allows us to understand qualitatively the mechanically-induced disappearance of a mechano-chemical pattern as σ^eff increases in figure 2(d), as well as the emergence of such a pattern when $\bar{\phi }$ increases as in figures 2(a) and (b), where χ^eff = 0.

To further study the mechano-chemistry of membrane-protein interactions, we then considered a setup that mimics controlled in vitro experiments, where a curvature gradient is established by pulling a highly curved membrane tether out of a vesicle [16, 18, 44, 88]. We consider the same vesicle size and reduced volume as in the previous examples, and gradually increase the distance between the poles. As in experiments [89], our simulations show that beyond an extension, the system breaks symmetry and a thin tube elongates from one of the poles, see figure 3(a)-i. Starting from this configuration, we load the vesicle with a uniform distribution of proteins with area coverage $\bar{\phi }=0.15$. The protein dynamics are similar to the previous examples, with a progressive enrichment of the highly curved tube over a time period of about the diffusive time-scale ${\tau }_{2}\approx 4\pi {R}_{0}^{2}/D(\bar{\phi })\approx 17\,{\rm{s}}$, where now proteins need to migrate a longer distance on average to one of the two poles, figure 3(c). At equilibrium and for this low protein coverage, proteins have barely modified the shape of the vesicle-tube system, figure 3(a)-iii. However, their presence has a noticeable mechanical effect in the force required to keep the tether in place, which decreases by more than two-fold, figure 3(b). This kind of behavior has been experimentally observed for BAR proteins and dynamin [16, 18, 88]. For this protein coverage, however, the amount of protein is insufficient to stabilize the tube and following the release of the displacement constraint, the tube retracts and the protein-rich domain dissolves into the vesicle, figures 3(a)–(c). At higher protein coverage, the proteins drawn to the tube are able to modify visibly its shape by inducing slight pearling, they further reduce the tether force, and the larger protein amount is able to stabilize highly curved and protein-rich protrusions upon release of the constraint, figures 3(d) and (e), in agreement with in vitro experiments [18, 88]. Furthermore, the transition to a strongly pearled protrusion upon force removal (iii–iv) closely mimics the shape transformations in membrane protrusions bent by crowding of the glycocalyx upon disassembly of enclosed actin filaments [12].

Zoom In Zoom Out Reset image size

Buds constitute a prototypical membrane motif, and are involved in endo/exocytosis [2] or in tensional buffering of the plasma membrane through caveolae [4]. Although the formation of such buds requires the synergistic interaction of multiple proteins and lipids, they can be abstracted as curved protein-rich membrane domains [76, 90]. To understand the fundamental mechanism of bud formation, we consider a flat discoidal patch of membrane covered with a homogeneous distribution of proteins ($\bar{\phi }=0.1$) at t = 0. We assume that the membrane is flat at its edge and that proteins cannot flow in or out from the boundary. We apply radial tractions at the boundary corresponding to an isotropic membrane tension of σ^eff ≈ 3 × 10⁻³ mN m⁻¹, at the lower end of membrane tension in mammalian cells [91].

According to equation (49), the parameters $\bar{C}=2/25\,{\mathrm{nm}}^{-1}$ and a^eff/a_p ≈ 0.4 mN m⁻¹ should make the initially uniform and flat state unstable, leading to a chemo-mechanical pattern. In agreement with this prediction from the linear stability analysis, our nonlinear, yet axisymmetric, calculations show the spontaneous formation of a budded protein-rich domain as shown in figure 4(a), reminiscent of caveolae. As shown in the figure, this process leads to a significant reduction in the projected area of the membrane patch, and thus during bud formation chemical energy is released to perform mechanical work against the applied tension.

Zoom In Zoom Out Reset image size

A critical function of caveolae is the mechano-protection of cells subjected to stretching of the plasma membrane [4]. These budded domains provide a membrane reservoir, which upon tension increase, can be released to buffer membrane tension and avoid lysis [92]. To test the ability of our model to reproduce this phenomenology, we suddenly increased membrane tension to 0.5 mN m⁻¹ within 0.6 s. As a result and in agreement with equation (49), the budded domain rapidly disassembles, leaving a flat patch with uniformly distributed proteins, figure 4(b), consistent with the increased mobility of caveolar components following tension-induced disassembly [4].

The curvature model considered up to this point, based on equation (6), assumes that protein cooperativity increases the spontaneous curvature of the surface. However, an alternative model can be conceived in which curved proteins have a stiffness of their own, which results in an effective density-dependent stiffness of the protein shell on the membrane [40]. In this case, the curvature energy of the composite membrane-protein system can be written as

$$\begin{eqnarray}&&{{ \mathcal F }}_{b}={\int }_{{\rm{\Gamma }}}\displaystyle \frac{\kappa }{2}{H}^{2}\,{\rm{d}}{S}+{\int }_{{\rm{\Gamma }}}\displaystyle \frac{{\kappa }_{p}(\phi )}{2}{\left(H-{C}_{p}\right)}^{2}\,{\rm{d}}{S},\end{eqnarray} \tag{ 50 }$$

where C_p is the intrinsic curvature of proteins and κ_p(ϕ) a density-dependent stiffness of the protein coat. The resulting membrane shape is hence the result of the competition between elastic bending energies of proteins and of the lipid bilayer. This model has been used to study the response of membranes with stiff protein coats [93]. Assuming a linear dependence of κ_p on the protein area fraction, ${\kappa }_{p}(\phi )={\bar{\kappa }}_{p}\phi$, and a chemical energy given by equations (7) and (8) as before, the chemical potential of proteins now takes the form

$$\begin{eqnarray}&&\mu ={\mu }_{0}+{a}_{p}{\bar{\kappa }}_{p}{\left(H-{C}_{p}\right)}^{2}+{k}_{B}T\mathrm{log}\displaystyle \frac{\phi }{{\phi }_{m}-\phi }+\chi \phi -{\rm{\Lambda }}{\rm{\Delta }}\phi .\end{eqnarray} \tag{ 51 }$$

Invoking Onsager's principle with the same dissipation potentials as before, we obtain an alternative protein transport equation

$$\begin{eqnarray}&&\xi \dot{\phi }+{\rm{\nabla }}\cdot \left[{a}_{p}{\bar{\kappa }}_{p}({C}_{p}-H){\rm{\nabla }}H\phi -\left(\chi \phi +\displaystyle \frac{{k}_{B}T{\phi }_{m}}{{\phi }_{m}-\phi }\right){\rm{\nabla }}\phi +{\rm{\Lambda }}\phi {\rm{\nabla }}({\rm{\Delta }}\phi )\right]=0,\end{eqnarray} \tag{ 52 }$$

where the density-dependent diffusion coefficient

$$\begin{eqnarray}&&D(\phi )=\displaystyle \frac{1}{\xi }\left(\chi \phi +\displaystyle \frac{{k}_{B}T{\phi }_{m}}{{\phi }_{m}-\phi }\right)\end{eqnarray} \tag{ 53 }$$

has the same structure as the one in equation (28). In contrast, the drift velocity ${{\boldsymbol{w}}}^{{\rm{drift}}}={a}_{p}{\bar{\kappa }}_{p}({C}_{p}-H){\rm{\nabla }}H/\xi$ is now qualitatively different, in that its sign relative to ∇H can change in space and time depending on the sign of C_p − H, whereas in the previous model (${{\boldsymbol{w}}}^{{\rm{drift}}}={a}_{p}\kappa \bar{C}{\rm{\nabla }}H/\xi$) it just depended on the sign of the constant $\bar{C}$. Focusing on the case in which both $\bar{C}$ and C_p are positive for concreteness, in the model based on equation (6) the drift term always favors protein transport towards regions of higher curvature, whereas in the model based on equation (50) this will be the case only as long as membrane curvature is smaller than the preferred protein curvature. As a result, in the model presented in this section the positive feedback between curvature and protein coverage stops once membrane curvature reaches C_p. Recall that, as discussed in section 6.2, in the previous model this positive feedback was only stopped by the saturation of protein coverage as ϕ approached ϕ_m. We tested this idea computationally by examining the equilibrium shape predicted by the model based on equation (50) with ${\kappa }_{p}=40{k}_{B}T$ and C_p = 1/25 nm⁻¹, of a slightly deflated vesicle with the same reduced volume v = 0.93 and an average protein concentration $\bar{\phi }$. Figure 5 shows that the system equilibrates at a state with a protein-rich domain where H ≈ C_p/2 and where ϕ ≈ 0.5 is lower than ϕ_m = 0.75. The curvature of the protein-rich domain is controlled by the competition of membrane and protein elasticity, which in turn depends on protein coverage. This calculation shows that in this model ϕ_m does not select the curvature of the protein-rich domain. To further confirm this, we observe that increasing ϕ_m = 1 does not change the equilibrium state. In contrast, increasing C_p by 17.5% and 20% led to more curved protrusions with a larger number or pearls. Eventually, as protrusions become increasingly concentrated in protein at high values of C_p, protein density reaches saturation and hence ϕ_m starts playing a role.

Zoom In Zoom Out Reset image size

Up to this point, we have assumed that proteins interact at the mid-plane of the lipid membrane. However, this approximation clearly breaks down for membrane proteins with bulky partially disordered domains [10, 17] or anchoring long polymers [12]. In this situation, the interaction between proteins leading to bending can be overwhelmingly dominated by the entropic repulsion of these bulky partially disordered domains or polymers, which interact a few nanometers away from the lipid membrane. In the case of polymers attached to a membrane, in addition to their positional entropy, one must account for the changes in conformational entropy of the polymers themselves, which can transition from a mushroom regime to a brush regime as local density increases [12, 33, 94]. Here, we consider the case of proteins with a bulky off-membrane domain, which may be partially disordered but whose conformation does not change significantly with lateral packing. In this situation, we can ignore the entropy of conformational changes of the proteins and their main contribution is due to mixing entropy. Their ability to bend the membrane is related to the fact that curvature modifies the area fraction of proteins at the off-membrane surface where they interact, see figure 6(a).

Zoom In Zoom Out Reset image size

Thus, ignoring reconfigurations of the disordered domain/polymer blob [33], we assume that the proteins interact on a surface ${{\rm{\Gamma }}}^{+}$ at a distance d from the surface representing the lipid membrane, Γ. The free energy of the proteins can then be written as

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal F }}_{p} & = & {\displaystyle \int }_{{{\rm{\Gamma }}}^{+}}\displaystyle \frac{{k}_{B}T}{{a}_{p}}\left[{\phi }^{+}{\rm{log}}({\phi }^{+})+({\phi }_{m}-{\phi }^{+}){\rm{log}}({\phi }_{m}-{\phi }^{+})\right]\,{{\rm{d}}{S}}^{+}\\ & & +{\displaystyle \int }_{{{\rm{\Gamma }}}^{+}}\displaystyle \frac{\chi }{2{a}_{p}}{\left({\phi }^{+}\right)}^{2}\,{{\rm{d}}{S}}^{+}+{\displaystyle \int }_{{{\rm{\Gamma }}}^{+}}\displaystyle \frac{{\mu }_{0}}{{a}_{p}}{\phi }^{+}\,{{\rm{d}}{S}}^{+},\end{array}\end{eqnarray} \tag{ 54 }$$

where ${\phi }^{+}$ is the area fraction of proteins on ${{\rm{\Gamma }}}^{+},{a}_{p}$ is the area on this surface of each bulky protein domain, and ϕ_m is the saturation area fraction of bulky domains on ${{\rm{\Gamma }}}^{+}$.

We next refer this energy to the bilayer mid-surface. If the separation between Γ and ${{\rm{\Gamma }}}^{+}$ is small, dH ≪ 1, then the area element of ${{\rm{\Gamma }}}^{+}$ is related to that of Γ according to

$$\begin{eqnarray}&&{{\rm{d}}{S}}^{+}\approx (1+{d}{H})\,{\rm{d}}{S}.\end{eqnarray} \tag{ 55 }$$

Denoting by ${n}^{+}={\phi }^{+}/{a}_{p}$ the number density of proteins on ${{\rm{\Gamma }}}^{+}$, the above relation shows that we can express the number density on Γ using the relation ${n}^{+}\approx n/(1+{d}{H})\approx n(1-{d}{H})$. This relation clearly shows how curvature changes density, as illustrated in figure 6(a) where positive/negative curvature increases/decreases ${n}^{+}$ at fixed n. Even if the area fraction does not make strict sense on Γ, we can formally define it on the bilayer mid-plane as $\phi ={a}_{p}n$ and hence

$$\begin{eqnarray}&&{\phi }^{+}\approx \displaystyle \frac{\phi }{1+{dH}}\approx \phi (1-{dH}).\end{eqnarray} \tag{ 56 }$$

Denoting by ${\boldsymbol{w}}$ the diffusive velocity of proteins relative to the bilayer velocity at the membrane mid-plane, protein balance of mass is still given on Γ by equation (5).

Using equations (55) and (56), noting that for small dH we have that $\mathrm{log}(1+{dH})\approx {dH}$, further assuming that $\mathrm{log}[{\phi }_{m}-\phi (1-{dH})]\approx \mathrm{log}({\phi }_{m}-\phi )+\phi {dH}/({\phi }_{m}-\phi )$, which holds true provided that ϕ is not too close to ϕ_m, and neglecting terms proportional to (dH)², we can rewrite equation (54) as integral over the lipid surface Γ as

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal F }}_{p} & = & {\displaystyle \int }_{{\rm{\Gamma }}}\displaystyle \frac{{k}_{B}T}{{a}_{p}}[\phi {\rm{log}}\phi +({\phi }_{m}-\phi ){\rm{log}}({\phi }_{m}-\phi )+{\phi }_{m}{d}{H}{\rm{log}}({\phi }_{m}-\phi )]\,{\rm{d}}{S}\\ & & +{\displaystyle \int }_{{\rm{\Gamma }}}\displaystyle \frac{\chi }{2{a}_{p}}{\phi }^{2}(1-{d}{H})\,{\rm{d}}{S}+{\displaystyle \int }_{{\rm{\Gamma }}}\displaystyle \frac{{\mu }_{0}}{{a}_{p}}\phi \,{\rm{d}}{S}.\end{array}\end{eqnarray} \tag{ 57 }$$

Combining this chemical energy with a simple Helfrich energy for the bare bilayer, ${{ \mathcal F }}_{b}=\kappa /2{\int }_{{\rm{\Gamma }}}{H}^{2}\ {\rm{d}}{S}$, we obtain

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal F }}_{b}+{{ \mathcal F }}_{p} & = & {\displaystyle \int }_{{\rm{\Gamma }}}\displaystyle \frac{\kappa }{2}H[H-2C(\phi )]\,{\rm{d}}{S}+{\displaystyle \int }_{{\rm{\Gamma }}}\displaystyle \frac{{k}_{B}T}{{a}_{p}}[\phi {\rm{log}}\phi +({\phi }_{m}-\phi ){\rm{log}}({\phi }_{m}-\phi )]\,{\rm{d}}{S}\\ & & +{\displaystyle \int }_{{\rm{\Gamma }}}\displaystyle \frac{\chi }{2{a}_{p}}{\phi }^{2}\,{\rm{d}}{S}+{\displaystyle \int }_{{\rm{\Gamma }}}\displaystyle \frac{{\mu }_{0}}{{a}_{p}}\phi \,{\rm{d}}{S},\end{array}\end{eqnarray} \tag{ 58 }$$

where we have defined the protein-induced spontaneous curvature by crowding of off-membrane bulky domains/polymer blobs as

$$\begin{eqnarray}&&C(\phi )=\displaystyle \frac{d}{{a}_{p}\kappa }\left(\displaystyle \frac{\chi }{2}{\phi }^{2}-{k}_{B}T{\phi }_{m}\mathrm{log}({\phi }_{m}-\phi )\right).\end{eqnarray} \tag{ 59 }$$

The formal similarity between this expression for ${{ \mathcal F }}_{b}+{{ \mathcal F }}_{p}$ and that obtained in section 3.1 is remarkable, the only difference being that before, the density-dependent spontaneous curvature was simply $C(\phi )=\bar{C}\phi$ and now it is given by equation (59). Note that the term proportional to ${(\bar{C}\phi )}^{2}$ in the bending energy of the previous model can be dropped by re-defining χ as χ^eff, and thus it does not affect the structure of the free energy.

In vitro experiments examining membrane bending by protein crowding [10, 17] confined proteins to membrane domains. Not being able to diffuse freely, proteins became increasingly confined, leading to severe membrane remodeling. See figure 6(b). Here, we consider this situation, by allowing ϕ to differ from zero only over a subdomain Γ_p ⊂ Γ. Over the interface given by ∂ Γ_p, we thus have an initial jump in protein area fraction of $\bar{\phi }$, the initial average area fraction over the subdomain. Across the interface ∂Γ_p, forces and moments need to be continuous. Since the energy depends on curvature, jumps in the normal are not allowed but finite jumps in curvature are [69, 95]. Using the expression for the bending moment derived in equation (34) adapted to the free energy density in equation (58), continuity of bending moments across the interface leads to the condition $\kappa (H-C(\phi )){| }_{i}=\kappa H{| }_{o}$, where the subscripts indicate whether the quantity is evaluated on the inside or on the outside of the interface. We thus conclude that the jump in mean curvature across the interface needs to coincide with the protein-induced spontaneous curvature inside the protein-rich domain, $[\kern-2pt[ H]\kern-2pt] =C(\phi )$. To test these ideas, we considered various average protein area fractions, $\bar{\phi }=\{0.1,0.3,0.6,0.9\}$, within a domain of diameter 250nm in a membrane patch of diameter 2.5 μm. We assumed that d = 1 nm and χ/a_p = 6 mN m⁻¹ (net repulsive protein–protein interaction). As shown in figure 6(c) increasing the number of proteins leads to increasing curvature, going from very shallow caps to buds, which in all cases very precisely follow the predicted relation for the jump in mean curvature, figure 6(d).