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

The function of biological membranes is controlled by the interaction of the fluid lipid bilayer with various proteins, some of which induce or react to curvature. These proteins can preferentially bind or diffuse towards curved regions of the membrane, induce or stabilize membrane curvature and sequester membrane area into protein-rich curved domains. The resulting tight interplay between mechanics and chemistry is thought to control organelle morphogenesis and dynamics, including traffic, membrane mechanotransduction, or membrane area regulation and tension buffering. Despite all these processes are fundamentally dynamical, previous work has largely focused on equilibrium and a self-consistent theoretical treatment of the dynamics of curvature sensing and generation has been lacking. Here, we develop a general theoretical and computational framework based on a nonlinear Onsager’s formalism of irreversible thermodynamics for the dynamics of curved proteins and membranes. We develop variants of the model, one of which accounts for membrane curving by asymmetric crowding of bulky off-membrane protein domains. As illustrated by a selection of test cases, the resulting governing equations and numerical simulations provide a foundation to understand the dynamics of curvature sensing, curvature generation, and more generally membrane curvature mechano-chemistry.


Introduction
Animal life is characterized by a hierarchical compartmentalization into separate functional units delimited by interfaces. At cellular and sub-cellular scales, the fundamental structure supporting compartmentalization is the lipid membrane. To form specialized organelles, lipid membranes must adopt a variety of shapes including spherical vesicles, sheets, tubes, and complex assemblies involving several of these elements [1]. Furthermore, biomembranes need to dynamically reconfigure their shapes and topologies to accomplish vesicular transport [2], cell division [3], or to unfold membrane reservoirs under stress [4]. The plasma membrane is also a mechanical organizer of the cell [5], and supports a variety of mechanotransduction mechanisms [6]. To perform all these functions, lipid membranes rely on a combination of solid-like mechanical properties, e.g.bending and stretching elasticity, in-plane fluidity, and the interaction with a myriad of curved membrane proteins. Molecularly, membrane proteins impinge curvature on the membrane through various mechanisms, which include curved scaffolding domains, wedge-like insertions [7][8][9], or asymmetrical crowding of bulky disordered domains [10] or anchored polymers [11] as in the glycocalyx [12] interacting at a distance from the bilayer mid-plane. Here, we refer to all these objects as curved proteins. Besides enabling shape remodeling, membrane fluidity allows curved proteins to migrate to regions of favorable curvature during curvature sensing [13], to cluster and reshape the membrane during curvature generation [7], while membrane tension controls the (dis)assembly of protein-rich curved domains [4].

Setup, kinematics and balance laws
We model the lipid bilayer as a material surface Γ parametrized by r t , q a ( ), where (θ 1 , θ 2 ) are Lagrangian coordinates labeling material particles and t denotes time. The model proposed here does not explicitly describe each of the two monolayers [52,53]. Accounting for the bilayer architecture may be pertinent to some molecular curving mechanisms, such as shallow insertions into one of the monolayers, whereas the model developed here could apply to interactions that equally affect both monolayers such as scaffolding or full insertion of transmembrane proteins [54]. Using standard differential geometry [55], we use this parametrization to define the tangent vectors at each material point as g r q = ¶ ¶ a a , which form the natural basis of the tangent space, and the metric tensor with covariant components g g g = ab a b · . The components g βγ of the inverse of the metric tensor follow from the relations g g d = ab bg a g . The unit normal to the surface is n g g g 1 2 =( ) , where g g det = ab . The local curvature of the surface k is characterized by the second fundamental form, which measures changes in the normal and whose components in the natural basis are n g . The invariants of the second fundamental form are its trace k H k g tr = = ab ab (twice the mean curvature) and its determinant k K k g det det = = ab bg (Gaussian curvature). Throughout the text, ∇ denotes the covariant derivative or surface gradient and · the surface divergence. The dynamics of the surface are determined by its Lagrangian velocity, which can be decomposed into tangential and normal components as As a result of this flow, the metric tensor changes with time. Its material time derivative, a partial derivative in our Lagrangian setting, is called the rate-of-deformation tensor of the surface [53,56], and includes the usual term accounting for deformation resulting from tangential flows, and a term accounting for deformation resulting from shape changes, which involves the normal velocity and the curvature. The rate of change of local area follows as We adopt a mean field description of proteins in terms of a scalar field f(θ α , t) measuring their local area fraction. In doing so, we assume that proteins are isotropic, or in a regime in which entropic effects dominate over orientational order. We leave for future work a general continuum theory accounting for nematic order, pertinent for instance to elongated membrane proteins with BAR domains [15,28,30].

Balance of mass
Denoting by ρ l (θ α , t) the lipid areal density, balance of mass of lipids requires d t tr 0 l l r r ¶ ¶ + = . Note that in this equation t l r ¶ ¶ 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 Since proteins are a diffusive species, balance of mass of proteins can be expressed as where w is the diffusive velocity relative to the Lagrangian coordinates and we have ignored protein sorption. Again, because of our Lagrangian parametrization, t f ¶ ¶ coincides with the material time-derivative ḟ .

Energetics, dissipation and power input
To describe the dynamics of protein-membrane interactions, we adopt a nonlinear Onsager's formalism of dissipative dynamics [57][58][59], according to which the time evolution of the system follows a minimization principle where energy release and dissipation compete. This nonlinear variational principle is useful to formulate complex models in a very transparent way, and has been invoked in a variety of contexts including low Reynolds number interfacial hydrodynamics [60][61][62], materials modeling [63][64][65], soft matter physics in general Figure 1. The state of the system, given by the shape of the bilayer mid-surface and by the protein area fraction, can evolve as a result of the membrane velocity (V ), tangential (v) and normal (v n ), and of the diffusive velocity of proteins relative to the membrane (w). [66], and membranes in particular [53,67,68]. In the present context, the state variables r and f determine the elastic energy (associated to bending) and the chemical energy (entropy and self-interactions) of the system, whereas the variables characterizing changes in the state of the system, here V and w, determine energy dissipation through shear viscous forces in the membrane and friction as proteins move relative to the membrane. Onsager's formalism provides a transparent method to derive the governing equations even in the presence of strong nonlinearity, here caused mechanically by the very large deformations of the membrane and chemically by molecular crowding in protein-rich domains.

Energetics
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 where κ is the bending modulus and 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] In this equation, k T B is the thermal energy, a p is the area on the membrane of a single protein so that f/a p is the number density, the term involving f 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 μ 0 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 Λbeing a material parameter. The length scale of this interaction is of the order c L | |, 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 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 k f f   · [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  with respect to f as where D W i denotes the partial derivative of W with respect to its ith argument. Because this variation is taken at fixed shape, 2 with n 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, a W W a C H C k T log , 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.

Dissipation
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, d tr 0 = , the dissipation potential accounting for in-plane shear can be written as where η is the in-plane shear viscosity of the membrane and the rate-of-deformation tensor d is given by equation (2) [67]. Protein transport is characterized by the collective protein velocity 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 w x -, where ξ is the molecular drag coefficient. By superposition in a dilute approximation, a collection of proteins characterized by local number density f/a p experiences a drag force per unit area given by w a p xf - . The associated dissipation potential can then be written as 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 mm [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:

Power input
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 where μ b is a fixed chemical potential for proteins at the boundary, e.g.maintained by a protein reservoir, t is a unit tangent vector along ∂Γ so that n F , n t =´t, F ν and F n are traction components at the boundary, M is a bending moment per unit length, and ṅ represents the material time derivative of the surface normal. We can express the last integral in terms of our dynamical variables v and v n using the relation · is the geodesic torsion of the boundary curve and kn n k = n · its normal curvature [77,78]. Each of these terms can be defined on Neumann parts of the boundary ∂Γ.

Onsager's variational principle
Onsager's recipe to derive the dissipative dynamics is to minimize the Rayleghian functional defined as where  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 Since Γ is a material surface, the rate of change of energy can be written as where Ẇ is the material time-derivative of the energy density and is given by 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 on ∂Γ [79]. As a result, the last integral over the edge also vanishes. Similarly, we can express the constraint of local area conservation as To obtain the governing equations of the system, we substitute the above relations in equation (19), minimize the Lagrangian with respect to w 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.

Transport of proteins
Minimizing the Lagrangian with respect to the diffusive velocity, Substituting the above relation in equation (5) and using the inextensibility of the lipid membrane we obtain which can be expanded recalling equation (12) and rearranged to yield the transport equation for proteins where we have defined the effective interaction between proteins as a C p eff 2 c c k = +¯. For vanishing spontaneous curvature (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 For χ=0, Λ=0 and f = f m , equation (27) reduces to the linear diffusion equation. For C 0 1 , the curvature gradient introduces a bias in the diffusion with drift velocity driving curved proteins along or against the curvature gradient depending on the sign of C . At steady state, drift and diffusive transport must balance and yield a divergence-free flux.

In-plane force balance
Minimization of the Lagrangian respect to the in-plane velocity, 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 f∇μ 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 where a a a a a g 2 is the deviatoric component of this rank-one tensor, and we identify the effective membrane tension (the hydrostatic part of the stress tensor) as 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-ofplane coupling mediated by curvature [81], this force density can induce out-of-plane forces and shape changes. At steady state, since w drift is in general different from zero, equation (31) shows that we can expect a nonuniform 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 f. 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 f = f m , this term becomes k T a B p f -, 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.

Force balance normal to the surface
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].

Boundary conditions
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 , , n t nˆˆa nd ŵ are Dirichlet data. For the chemical problem, the Dirichlet boundary condition is prescribing w 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: 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 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 .

Time-scales of shape and protein density relaxation
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 =ηℓ 2 /κ , where ℓis the characteristic size of the disturbance. The characteristic time for protein diffusion is k T where we have used the common estimate for the bending stiffness k T 20 B k » . In turn, ξ is related to the membrane surface viscosity η through the Saffman-Delbrück theory [82], which states that 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 1 40 ). 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 selfinteraction, 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.

Axisymmetric formulation and numerical approximation
Under the assumption of axisymmetry, pertinent to many structures resulting from protein-membrane interactions, the shape of the membrane can be parametrized in terms of the distance to the z axis and the z coordinate of its generating curve (ρ(u, t), z(u, t)), where u is a Lagrangian coordinate labeling material particles in the interval [0, 1] and t is time. Protein area fraction does not depend on the azimuthal angle and thus can be expressed as f(u, t). For closed surfaces, smoothness of the surface at the poles is guaranteed by the conditions  [74]. Noting that the element of area is S a u d 2 d p r = , this expression allows us to compute the bending energy ]in equation (6). As a reference surface Ḡ with local radius r and speed ā is mapped to the current surface Γ with radius ρ and speed a, the local areal stretch at each point is J a a r r = ( ) (¯¯). Thus, membrane inextensibility can be expressed as J=1, or a a r r =¯¯. As shown in [74], the membrane dissipation potential in equation (13) is the surface Laplacian of f. To perform numerical calculations, we use a Galerkin finite element approach based on a B-Spline approximation of the different fields. We numerically represent the state variables as . This approximation provides C 2 continuity, enough for our formulation, which requires at least C 1 continuity for square integrable curvatures and protein Laplacians.
To move forward in time, we adopt a staggered approach in which we first evolve the protein density field at fixed membrane shape, and then update shape at fixed protein distribution. To obtain the concentration of proteins f n+1 at time t n 1 + , we assume a given shape of membrane {ρ n , z n } at time t t t n n n 1 = -D + , use a backward Euler approximation to discretize equation (40) in time, multiply this equation with a test function ψ (u), integrate over the surface and integrate by parts. For simplicity of our exposition, we assume no flux of proteins through the boundary to obtain We note that we have further integrated by parts the term involving H¢ to lower the smoothness requirements of the theory. Replacing equation ( Since the matrix K IJ depends on the unknown, the system of equations in equation (44) is nonlinear and we solve it using Newton's method.
To solve the mechanical problem, we fix protein area fraction to f n+1 and write down an incremental Lagrangian accounting for the rate of change of the free energy, for membrane dissipation, and for local area and volume constraints ; , The Lagrange multiplier n 1 s + is also discretized in space using B-splines. However, rather than cubic, we use quadratic B-spline basis functions for this field to obtain a stable formulation [84,85]. To move forward in time, we obtain the mechanical unknown at time t n+1 by numerically solving the algebraic optimization problem s described above, this Lagrangian method will in general lead to significant distortions of the numerical grid. For robustness and accuracy of the numerical methods, we periodically perform mesh reparametrizations of the generating curve.

Selection of parameters
We choose as the energy scale the bending rigidity of the membrane k T 20 B k = . As the length and time-scales, we choose ℓ 0 =50 nm and k T Considering a membrane viscosity of 5 10 N s m 9 1 h =´-and the relation ξ=2πη discussed in section 4.6, we obtain that 3 10 N s m 8 1 x »´--, the diffusion coefficient of proteins is D k T 1.3 10 m s 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 where Γ s is the entire surface excluding a small region near the poles and γ=0.1k B T. We checked that this potential only affected the solutions close to the neck.

Curvature sensing and generation starting from a prolate vesicle
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 c k = -¯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 0 f  = at t=0, protein transport is initially due exclusively to gradients in curvature with the diffusive velocity coinciding with the drift velocity w a C H 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 Subsequently, the dynamics are governed by a competition between drift and diffusive transport, driving proteins towards equilibrium over a time scale . Thus, we estimate that the total time scale of relaxation is given by τ≈τ 1 +τ 2 ≈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,  D , 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 . 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 f m =1 and membrane shape does not change.
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 0.35 f = , 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 Cf 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 f 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 f m in setting the equilibrium state, we repeated the previous simulation considering different values of f m . Figure 2(c) shows that the depth of the pearling cascade is indeed controlled by this parameter. For f 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 f m the system does not even pearl. Thus, in a model governed by the energies in equations (6) and (7), protein saturation controlled by f 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 0 =13 Pa and σ eff ≈0.003 mN m −1 . 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 1 =55 Pa, p 2 =65 Pa). Beyond p 3 =250 Pa corresponding to σ eff ≈0.064 mN m −1 , 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 0 f < (¯) , 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 This equation allows us to understand qualitatively the mechanically-induced disappearance of a mechanochemical pattern as σ eff increases in figure 2(d), as well as the emergence of such a pattern when f increases as in figures 2(a) and (b), where χ eff =0.

Sensing on a tube and shape stabilization
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 , 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 proteinrich 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].

Bud formation and tension-induced dissolution
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 ( 0.1 f = ) 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 f x has the same structure as the one in equation (28). In contrast, the drift velocity w  Figure 5 shows that the system equilibrates at a state with a protein-rich domain where H≈C p /2 and where f≈0.5 is lower than f 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 f m does not select the curvature of the proteinrich domain. To further confirm this, we observe that increasing f 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 f m starts playing a role.

Membrane bending by protein crowding
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).
Thus, ignoring reconfigurations of the disordered domain/polymer blob [33], we assume that the proteins interact on a surface G + at a distance d from the surface representing the lipid membrane, Γ. The free energy of where f + is the area fraction of proteins on a , p G + is the area on this surface of each bulky protein domain, and f m is the saturation area fraction of bulky domains on G + . We next refer this energy to the bilayer mid-surface. If the separation between Γ and G + is small, dH = 1, then the area element of G + is related to that of Γ according to Denoting by n a p f = + + the number density of proteins on G + , the above relation shows that we can express the number density on Γ using the relation n n dH n dH ). 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 midplane as a n p f = and hence dH dH 1 1 .
Denoting by 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) ) , which holds true provided that f is not too close to f m , and neglecting terms proportional to (dH) 2 , we can rewrite equation (54) as integral over the lipid surface Γ as Corresponding jumps in mean curvature as a function of the density dependent spontaneous curvature C(f).
The formal similarity between this expression for b p   + and that obtained in section 3.1 is remarkable, the only difference being that before, the density-dependent spontaneous curvature was simply C C f f = ( )¯and now it is given by equation (59). Note that the term proportional to C 2 f (¯) 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 f 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 f , 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 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, H C f =   ( ). To test these ideas, we considered various average protein area fractions, 0.1, 0.3, 0.6, 0.9 f = { } , 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 −1 (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).

Conclusions
We have presented a nonlinear and self-consistent continuum model for the dynamics of membranes interacting with curved proteins. Our theory describes a biologically important instance of chemo-mechanical self-organization leading to surface shape dynamics, which coexists in cells with alternative shape patterning mechanisms [96]. By combining elementary ingredients into a nonlinear Onsager's formalism, we have systematically derived fully nonlinear governing equations exhibiting a tight interplay between geometry, protein transport, and mechanics. Previous simpler models appear as specialized limits of our theory. Our numerical simulations have demonstrated the ability of the model to describe curvature sensing, generation, stabilization, and tension-induced disassembly of protein-rich curved domains. We have developed three versions of the model. We have shown that a common model where spontaneous curvature is proportional to protein density develops a positive feedback between curvature and protein density, only stabilized by protein saturation. An alternative model accounting for the bending elasticity of the protein coat does not exhibit this feature. Finally, a variant of the model where bending is induced by crowding of bulky off-membrane protein domains is formally equivalent to the first model, albeit with a nonlinear relation between protein-induced spontaneous curvature and protein density. The work presented here can be the background for further studies, e.g.accounting for orientational order to model anisotropic proteins such as BAR domains [97], for the interaction of multiple curvature-active species, or developing computational methods to treat the governing equations in 3D [53]. We end by discussing the applicability of the model presented here. While our approach can efficiently treat large numbers of proteins during slow diffusive processes (over minutes and spanning microns), whether continuum models can be quantitative in small systems involving tens of proteins over 10s of nanometers remains to be systematically tested. Furthermore, it is unclear at this moment to what degree the chemical specificity required to understand the biophysics of membrane-protein interactions can be captured by the parameters of our model or variations of it. To address these issues, comparison with coarse-grained molecular dynamics may be particularly interesting, as these models can more easily connect with truly specific atomistic models while accessing larger systems for longer times. 1 , Taking the variation of the dissipation potential, a direct calculation allows us to identify g . Similarly, we can identify the component due to inextensibility as g g s s =  for r x y i j 0 = + and q n n , Qualitatively this condition is similar to that in equation (B.13). However, there is one subtle difference since now a eff and σ eff are independent of the spontaneous curvature of the proteins, which was not the case before.