Collective membrane dynamics emerging from curvature-dependent spatial coupling

Membrane curvature has been recognized as an active participant of fundamental biological processes including vesicular transport and organelle biogenesis, but its effects on membrane remodeling are typically local. Here we show membrane curvature plays a critical role in propagating cortical waves and modulating mesoscale dynamics in living cells. We employ a membrane shape-dependent mechanochemical feedback model to account for the observed oscillatory travelling waves of Cdc42, F-BAR proteins and actin. We demonstrate that oscillatory membrane shape changes accompany and are required for such spatiotemporal patterns. In addition, modulating the curvature preference of the F-BAR proteins or membrane tension perturbs wave propagation. Our findings identify a distinct role of membrane curvature in mediating collective dynamics of cortical proteins and provide a molecular framework for integrating membrane mechanics and biochemical signaling in the context of subcellular pattern formation.


Introduction
Dynamic patterns as waves are widely observed in many biological systems. They reflect the pulse of underlying signal transduction [1]. Actin waves on plasma membrane, which underlie many important cellular functions [2], are especially intriguing. Many different models have been proposed for waves of actin, most of which are based on reactiondiffusion mechanism [3,4] or its variants with additional factors including cytoplasmic flow [5], myosin-dependent actin transport [6,7] and motor-independent treadmilling [8,9]. Models with actin transport or treadmilling are likely not directly applicable to the actin waves observed in Dictyostelium [10,11] or immune cells [12], where sequential cycles of dendritic actin assembly and disassembly are involved [13]. These actin-centric wave models differ significantly in terms of whether actin provides both positive and negative feedbacks [14], positive feedbacks alone with an unknown diffusing inhibitor [15], or negative feedbacks alone with an upstream activator [4].
Commonly in these models, the plasma membrane is treated as a passive and flat twodimensional manifold. However, it has been recognized that membrane shape change can feed back onto biochemical pathways in governing many cellular processes [16][17][18][19][20].
Interestingly, membrane curvature was theoretically predicted to induce traveling wave in active membrane when little was known about the molecular mechanisms of membrane curvature sensing or generation [21]. In this mechanical wave model, both types of curvature (concave and convex) are necessary and they act as positive and negative feedbacks, respectively. With the increasing interests in the involvement of actin waves in cellular processes involving membrane morphogenesis, such as cell migration, spreading and division, membrane geometry becomes an attractive parameter for modeling [22,23], as well as for experiments [24,25]. Membrane shape changes were also contemplated for actin wave-related phenomena such as cell edge protrusions [26,27] and dorsal ruffles, a wave-like propogation on the apical side of cell which involves membrane protrusions [28]. Because of a lack of knowledge of specific proteins for modulating membrane shape change, these models also differ in the specific requirements of mechanisms for curvature sensing or generation, where changes in membrane shape could be a passive result of forces applied by the actin filaments. More importantly, protrusion of the membrane by actin polymerization is an essential aspect of positive feedback and force generation mechanism. Therefore, it is not clear whether similar mechanisms could apply to actin waves propagating on the ventral surface of the cell, where membrane protrusion is presumably inhibited by the substrate. The role of membrane shape in actin waves at the basal surface, either passive or active, remains obscure due to a lack of direct experimental observation of plasma membrane deformation in these waves [25].
While all those models can generate morphologically similar wave patterns, they differ significantly in the mechanisms of wave propagation. In particular, whether actin plays an activating or antagonistic role, and whether membrane mechanics plays an active function in driving waves remain as fundamental open questions [14,29]. In this work, we combine theory and experiments to study the nature of collective membrane and protein waves in mast cells. Here we propose a ventral actin wave model in which the waves are driven by membrane shape undulation and fed back by curvature-dependent biochemical pathways, instead of driven by protrusive forces generated by actin polymerization. We show experimentally that membrane curvature is indispensable and that the sensitivity of cortical proteins to membrane shape renders the ventral cortex itself an excitable medium.
Consequently, membrane shape deformation potentiates long-range traveling waves that emerges from this cortical excitability.

Propagation Speed of Waves does not Scale with Actin Polymerization Rate
A number of actin wave models argue that protrusive force driven by actin polymerization on the membrane is the activating factor for travelling wave formation [14,26,29]. If the spatial propagation is driven by actin polymerization, wave propagation speed is expected to be proportional to the actin protrusion velocity, which depends on the rate of actin polymerization [2]. Our previous findings introduced a number of explicit components of the actin wave in mast cells [30]. Upon antigen stimulation, not only actin, but also a few membrane-interacting proteins, including N-WASP, Cdc42 and F-BAR domain-containing protein FBP17, exhibit coordinated traveling wave behavior within the ventral cortex. These proteins allow us to investigate the effect of perturbing actin turnover in cortical waves using reporters other than actin itself. With low doses of Latrunculin A (LatA), which prevents actin polymerization by sequestering G-actin, we observed that in mast cells, the oscillation period increased, but FBP17 waves persisted with the same propagation speed ( Figure 1A). Lack of the dependence of FBP17 wave speed on actin polymerization rate suggests that the waves are not driven by actin polymerization. In addition, lifetime of FBP17 puncta was prolonged (red arrow in Figure   1A). Together those data suggest that the actin likely plays an antagonistic role in the cortical waves of mast cells.

Mechanochemical Model Recapitulates Wave Propagation
We decided to identify mechanisms that could positively activate wave propagation on the ventral membrane. Considering the facts that F-BAR proteins are members of BAR superfamily proteins that can sense and generate membrane curvature [31][32][33], and the theoretical concept that membrane shape can feed back to cortical protein recruitment in waves [18,19,21,28,34], we reason that these cortical proteins may mediate membrane shape changes that are important to rhythmic propagation. We set out to determine theoretically whether a model involving membrane shape propagation and curvaturedependent biochemical feedbacks can recapitulate our experimental observations quantitatively.
The mathematical essence of the model is formulated into a set of partial differential equations that govern membrane shape changes, cortical reaction and molecular diffusions (Equations [2.1] - [2.8] and Section 2 of Materials and Methods). Briefly, the model assumes a transient activation signal that increases the local amount of GTP-bound Cdc42 at the cortex ( Figure 1B). Active GTP-Cdc42 at the membrane is assumed to recruit effector proteins N-WASP and F-BAR domain-containing proteins (e.g., FBP17, CIP4 or Toca-1) [30,35,36], which we refer to a generic entity as F-BAR. Waves of active Cdc42 were synchronized with the FBP17 and N-WASP waves [30]. A similar activation signal may be evoked by the intrinsic cellular state via a self-assembly process that does not require external stimulation, as in the case of spontaneous travelling waves [30]. Oscillations in the model arise due to the non-linear mechano-chemical pathways downstream of Cdc42 activation.
Cortical Cdc42/N-WASP/F-BAR (CWF) complex could elicit two responses in the model. First, when the F-BAR level exceeded a threshold value, F-BAR would generate membrane shape changes. We postulate that the resulting local curvature induces further F-BAR recruitment to the cortex due to their curvature-sensing ability [37], which serves as the basis of positive feedback between activation of CWF complex and membrane shape changes. Second, active N-WASP recruits Arp2/3 complex and promotes actin polymerization [38][39][40]. This is based on the observation that waves of active Cdc42 preceded actin waves [30]. Increased actin polymerization decreases local membrane deformability and antagonizes F-BAR cortical recruitment [41][42][43]. This results in a delayed negative feedback of actin polymerization to cortical recruitment of CWF complex, forming the basis of the temporal oscillations. The model shows that oscillations arise with parameters that are within the range of measured values under physiological conditions (Tables S1 and S2). In agreement with experiments, it successfully recapitulates the local oscillation of cortical F-BAR and actin with phase shift by total internal reflection fluorescence (TIRF) microscopy ( Figure 1C). According to the model, travelling wave appearance of cortical proteins can develop as the out-of-plane membrane undulation conveys the membrane shape deformation in the activation patch to the neighboring cortical area, where such deformation invokes mechanochemical feedback to locally recruit proteins from the cytoplasm ( Figure 1D).
As such, wave propagates well beyond the original cortical activation site, like a ripple moving in the pond. The model predicts two types of propagation depending on the rate of local cortical activity ( Figure S1A). For the "ripple wave", the wavefront displays saltatory instead of smooth movement. That is, both the speed and the cortical protein level at the wavefront oscillate, and there is a phase shift between them ( Figure S1A, Right). The saltatory behavior of the wavefront stems from the fine balance between rate of cortical protein recruitment and rate of the attended membrane shape change. Away from the balance arises the "continuous wave", in which the wavefront travels at a roughly constant speed and cortical protein intensity ( Figure S1A, Middle). The buildup of cortical protein, which drives nucleation, is in synchrony with wave propagation. If the cortical protein level at the wavefront is below threshold, the wave is predicted to decay at the edge of the cortical patch and eventually die out ( Figure S1A, Left). Those predictions match with experimental observations ( Figure S1B).
The delayed actin-dependent negative feedback supports a travelling wave regime in a wide range of parameter space ( Figure S2). While wave propagation has contributions from both lateral diffusions and the local mechanochemical reactions through recruitment from the cytosol, our calculation shows that ~90% of the increase in the local protein concentration at the wave front stems from the recruitment reactions from cytoplasm, while the protein lateral diffusion contributes to ~10% ( Figure S3). This suggests that the wave propagation reflects more of local protein assembly from cytosol than protein diffusional drift along the cortex. This prediction is consistent with the lack of laternal movement of the individual FBP17 puncta in the wave ( Figure 1E). The model also recapitulates the dependence of oscillation period and the relative insensitivity of propagation speed on actin polymerization rate ( Figure 1F), in consistent with LatA experiments ( Figure 1A).

Membrane Shape Oscillations Accompany Traveling Wave Propagation
We next attempted to further validate the model experimentally by testing several predictions. One defining feature is that membrane shape changes are involved in wave propagation. Figure 2A shows membrane shape deformation accompanied the traveling waves over time and space in the model. Cortical protein recruitment was in synchrony with membrane shape change, and marks the wave front propagation ( Figure 2B). Note that, the local F-BAR level closely followed that of membrane curvature, but not the membrane height ( Figure 2B). To obtain information on the local membrane shape and its relation to the local rhythm of cortical protein recruitments, we simultaneously imaged FBP17 by TIRF microscopy and the distance between cell membrane and the substrate by surface reflection interference contrast (SRIC) microscopy [24]. We found that local membrane shape oscillated and the deformation strongly correlated with the location of the wave marked by cortical FBP17 density ( Figure 2C). In addition, both theory and experiment showed a longer duration of membrane height oscillation compared to that of F-BAR oscillations ( Figure 2D), due to the delay between changes in membrane curvature and membrane height ( Figure 2B). Membrane oscillation had a phase lag relative to F-BAR ( Figures 2D, S4A and S4B). Importantly, F-BAR was in the same pace with the rate of membrane height changes (dh/dt) at any fixed location ( Figures 2D and   S4C), consistent with the reciprocal interaction of F-BAR recruitment and local membrane curvature underlying the wave propagation. Moreover, actin was anticorrelated with the rate of membrane height changes (dh/dt) ( Figure S4A and S4C), indicating its antagonistic effects on membrane.

Curvature Preference of F-BAR is Critical for Wave Propagation
As membrane shape did change detectably along the wave propagation path, we further determined whether membrane shape changes were necessary for wave propagation. In the model, wave propagation requires a threshold level of curvature-dependent cortical protein recruitment, below which it can not sustain the residual membrane shape deformation at the edge of the original cortical patch of activation ( Figure 3A). We first examined the effect of removing F-BAR proteins on the formation of traveling waves.
We reduced the level of endogenous curvature-generating proteins by shRNA. Knocking down FBP17 or its homolog CIP4 alone did not significantly reduce the percentage of the cells with traveling waves, but double knocking down (DKD) FBP17 and CIP4 abolished the formation of traveling waves ( Figure 3B). This indicates some redundancies among the Toca family of F-BAR proteins and demonstrates their essential role as a collective entity of the Cdc42 interacting F-BAR proteins in the traveling waves.
To determine whether the F-BAR domain of FBP17 was necessary or sufficient for the FBP17 wave, we generated multiple versions of truncated FBP17. Deletion of any of the three structural domains (F-BAR, HR1, or SH3) abolished the wave-like appearance of cortical FBP17 in wild-type cells ( Figure S5), indicating that the F-BAR domain was necessary, but not sufficient for generating the traveling waves. ΔSH3 and ΔHR1 mutants each formed some stable puncta on cell membrane, but ΔF-BAR did not ( Figure S5).
Because the SH3 domain and HR1 domain interact with actin machinery, the formation of stable FBP17 puncta in ΔSH3 and ΔHR1 mutants suggests the requirement of a link to actin for dissociation of FBP17 from cortex, thus supporting the actin-mediated negative feedback proposed by our model.
Because F-BAR proteins could function in regulating actin dynamics in addition to their membrane-remodeling ability, the indispensable role of F-BAR proteins in wave propagation did not necessarily prove a requirement for membrane shape changes. To further explore whether the wave propagation requires changes in membrane shape, we generated a series of domain-swapping mutant proteins which differed from FBP17 only in their curvature preferences and tested whether they could functionally replace FBP17 in the waves. Two types of mutants were designed based on biophysical properties of known curvature-generating proteins, which covered a range of curvature preferences from 0.1 nm -1 to 0.01 nm -1 ( Figure 3C). First, we replaced the F-BAR domain of FBP17 with other lipid-binding domains of different shapes: F-BAR domain of FCHo1, BAR domain of Endophilin2, ENTH domain of Epsin1 ( Figures 3C and 3D). These domains represented three major classes of curvature-generating domains (F-BAR, BAR and ENTH) [70][71][72]. Phospholipid-binding but curvature-insensitive PH domain of phospholipase C [73] or membrane-targeting sequence of tyrosine-protein kinase Lyn (Lyn10) were introduced as controls. Among these, the mutant with the constitutively membrane-targeting Lyn10 did not form waves in wild-type (WT) cells, consistent with a requirement of F-BAR dissociation from the membrane for waves ( Figure 3D, top). The rest of the domain-swapped mutants still localized to the waves in WT cells, indicating that these mutants were properly folded functional proteins. However, in DKD cells, only other F-BAR (F-BAR of FCHo1) could rescue wave formation ( Figure 3D, bottom). The fact that the BAR (from Endophilin2) or ENTH (from Epsin1) domain-swapped mutant could not rescue wave formation suggests that a specific range of curvature (similar to or lower than that of the F-BAR domains) is required for wave formation.
For the second type of mutants, we introduced single point mutations to perturb the curvature preference of the F-BAR domain of FBP17 ( Figures 3C and 3E). Mutants K56E, K33E and F276D had paralyzed tubule-generating ability in living cells [31,51].
Of particular interests are mutants K166A and K66E, which generated narrower (higher curvature) and wider tubules (lower curvature) compared with wild-type F-BAR domain, respectively, in in vitro liposome tubulation assays [51]. Both mutants K166A and K66E were localized in waves in WT cells ( Figure 3E, top), confirming that they were functional. However, only K66E point mutant of FBP17 (of shallower curvature preference) could rescue wave formation in DKD cells ( Figures 3E, bottom). Together with the domain-swapping data ( Figure 3D), these results indicate that the shape of the lipid-binding domain of FBP17 is important for wave formation, regardless of the origin of the curvature preferences being the intrinsic protein shapes or the conformation of protein oligomers. We conclude that membrane shape-mediated feedback is actively involved in the cortical wave propagation.

Wave Propagation Requires Optimal Plasma Membrane Deformability
The requirement of membrane shape-mediated feedback confirmed the mechanochemical nature of the wave propagation. An independent test for this was on the mechanosensitivity of the waves. As membrane deformability critically determined rhythmic propagation of cortical protein recruitment in our model, the wave propagation was predicted to diminish with increasing or decreasing membrane tension ( Figure 4A).
Hence, we investigated the effect of membrane mechanics on the waves. We conducted osmolarity shock experiments by cyclically perfusing cells with hypotonic or hypertonic buffers; these treatments increased or decreased membrane tension respectively, which subsequently affected membrane shape deformability. The kymographs show that osmolarity changes reversibly inhibited rhythmic wave propagation ( Figure 4B). Osmotic treatment could lead to additional effects including cell volume and cytosolic concentration changes ( Figure S6). To circumvent these potential side effects, we used the surfactant deoxycholate (DC) to specifically reduce the plasma membrane tension without changing cell volume [74]. We observed that FBP17, actin, and membrane shape waves disappeared ( Figure S7) in a DC dose-dependent manner ( Figure 4C). Such effects were reversible on the time scale of seconds upon DC washout, eliminating possibilities of permanent membrane compositional changes. Waves experiencing moderate DC or hyper-osmotic treatment displayed no significant changes in propagation speed, but prolonged oscillation periods ( Figure 4D), indicating that the shape-mediated mechanochemical feedback quantitatively determines cortex excitability.

Discussion
Taken together, our results present a biological scenario where membrane shapemediated mechano-chemical feedback is critical for spatiotemporal pattern formation in the dynamical systems. Our data agree with the sequential actin assembly model proposed for actin waves in Dictyostelium [10] and Xenopus [75], which is in contrast to collective pattern formation of actin that requires bulk actin flow or transport of actin filaments [76]. Waves are not driven by the protrusive forces of actin on the membrane, as suggested by previous models [8,28,77]. The antagonistic role of actin on cortical wave is similar to what was proposed in other systems [4,12,75,78], but the activating signal in those models were assumed to be purely chemical. Here both our model and experiments show that the plasma membrane acts as an active medium for cortical traveling waves. Membrane shape changes establish a spatial cue for cortical protein recruitment, which reciprocally governs the cortical protein dynamics, rendering the plasma membrane an excitable system. Changes in membrane shape and consequently wave propagation depend sensitively on both the dynamic protein membrane interaction and membrane deformability. Softer membrane confers a faster F-BAR-mediated membrane shape deformation that reciprocally speeds up F-BAR recruitment to the cortex, which ends up deforming the membrane in a highly localized manner. Because membrane tension is low, smoothing of this deformation is slow. This introduces a further delay in the actin polymerization-mediated negative feedback on F-BAR, prolonging the oscillation period. However, if the membrane tension is reduced even further, it will not be able to smooth out the membrane deformation on a relevant timescale and consequently the cortical rhythm is disrupted. In this latter case, F-BAR puncta remains non-oscillatory and non-propagating, which phenocopies the nonpropagating puncta induced by the BAR domain mutants. Conversely, stiffening membrane (as in hypotonic treatment) inhibits membrane shape changes and F-BARprotein recruitment. When the F-BAR level is lower than the threshold value, cortical oscillation stops, and wave decays.
Travelling wave has been proposed to correlate information processing at integrated level in giant single cells or multi-cellular systems [79][80][81]. Recently, travelling waves of actin and GTPase including Rho [75], Rac [82] and Cdc42 [83] have emerged as hallmarks and likely regulators for cell cycle regulation, cell migration and chemotaxis.
Mechanochemical sensitivity of F-BAR proteins also appears to be important for cell polarity [84]. Such curvature-mediated spatial coupling, manifested as travelling waves, theoretically allows cell to rapidly scan its surface for mechanical and chemical cues.
This provides an efficient mechanism of cortical signal transmission that can be widely applicable in diverse cellular contexts. New Haven, CT) as previously described [30]. All plasmids were sequenced to vindicate their integrity.

Section 2. Model equations and parameters
The model describes rhythmic propagation emerging from the interplay between membrane-bound chemical reactions and membrane shape changes. For cortical protein recruitment from and turnover into the cytoplasm, the cytoplasm is treated as an unlimited reservoir. Once recruited to the ventral cortex, these proteins diffuse laterally; notably, cortical diffusion is much slower than cytoplasmic diffusion [56,86]. The model assumes the membrane as an elastic sheet. As the model concerns the cortical dynamics occurring on the ventral membrane, membrane shape is governed by membrane tension, bending modulus, membrane-substrate adhesion, and mechanical action of membranebound proteins, including F-BAR-mediated membrane shape deformation and actinmediated enhancement of membrane tension. The resulting membrane shape promotes F-BAR recruitment, which is depicted with Michaelis-Menten-like kinetics.
We determine the dynamical evolution of the system by integrating these equations from an initial condition: a flat membrane patch (40 µm across) and no cortical proteins. In the simulation, this membrane patch represents the ventral side of the plasma membrane, beyond which the model imposes a pinning boundary condition that ramps up a large resistance force, restoring the membrane height back to baseline, reflecting the hindrance effect of the cell edge, as the traveling wave in our system is observed on the cell ventral side. With an input of a localized transient pulse of GTP-Cdc42, the model output at each time step is membrane shape and local protein densities.
In the following, we describe the formulation of model equations. We focus on the simplest mechanism and only include the most essential components suggested by experiments. These components are: Cdc42, N-WASP, F-BAR, actin, Arp2/3, and membrane shape. Below, we will translate this qualitative picture into a set of coupled partial differential equations (PDEs). All the chemical reactions are modeled by Michaelis-Menten type kinetics, in which the local membrane shape influences the on rate of F-BAR. Conversely, the local F-BAR level sculptures the shape of the local membrane, and the local actin level modulates the membrane tension. In this way, the membrane mechanics is coupled with the local chemical reactions. These PDEs thus depicts the spatio-temporal dynamics of the mechanochemical coupling of the membrane shape and the membrane-bound chemical reactions. The model parameters are listed in Tables S1 and S2.  A 0 Threshold actin level to modulated membrane tension 0-1 0.2 [42,43] Cdc42 dynamics: N-WASP dynamics: F-BAR dynamics: Cortical recruitment of F-BAR by Cdc42 and N-WASP Here, the term k on Fh ⋅ e Ω h+ /Ω 0 ( ) − 1 ( ) depicts the membrane curvature-dependence of F-BAR which F-BAR cortical recruitment is significantly increased. We would like to emphasize that the F-BAR in the model refers to a collective functional entity of multiple F-BAR domain proteins, instead of referring to any individual F-BAR protein. This is consistent with our double knockdown experiment results in Figure 3B. Consequently, the characteristic membrane curvature Ω 0 is within the range of the values measured/inferred from different F-BAR proteins. Because F-BAR is more flexible than BAR-domain proteins, it could bind to the membranes with different curvatures that correspond to vesicle sizes ranging from 60 nm to microns in diameter [31,49,50]. An in vitro experiment demonstrated this curvature sensing of FCHo (a F-BAR protein), which corresponded to the size of lipid vesicle ~ 200-400 nm in diameter [49,50]. In the model, we thus chose the nominal value of Ω 0 to be ~ 1/200 nm -1 . In addition, our phase diagram study suggests that the model essence persists over a broad range of this threshold curvature Ω 0 ( Figure 3A).

Arp2/3 dynamics:
Note that N-WASP is known to activate Arp2/3 complex as a dimer [40]. Therefore, the Hill coefficient is 2.
Actin dynamics: We note that Arp2/3-mediated branched actin polymerization is autocatalytic. It involves nucleation phase with latency followed by rapid actin polymerization and, hence, is highly non-linear [38]. The combination of equations (2.2), (2.4) and (2.5) phenomenologically describes the nonlinear dynamics of actin polymerization with an effective time delay that lags behind the N-WASP activation.
Membrane shape dynamics: Equation ( drag from outside, whose value is measured to be ~ 2×10 9 Pa . s/m [48]. In addition, the functional derivative δ F δh at the membrane equals to the derivative of the F-BAR concentration field in the normal direction of the membrane surface. The F-BAR concentration in this direction has a gradient low in cytoplasm and high at the plasma membrane, so it is negative. In the simplest approximation, we assume that this gradient at the membrane is a constant -g F with g F >0. Combined together,  [88]. We therefore use the hill coefficient = 3 in our model. The qualitative feature of our model remains the same for any hill coefficient > 2 (data not shown). It also has been measured that cortical actin polymerization could increase the measured plasma membrane tension by 2-10 folds [42,43,68,69], which sets the range of the a A in equation (2.8).

Computational scheme
In the simulation, a free boundary condition was imposed. The initial condition was: all chemical components were at the equilibrium; membrane was at the resting state. At time zero, an activation signal was imposed to the system. We    Figure 2A). The membrane height was normalized relative to the highest membrane height. Right: Experimental measurement of the spatio-temporal changes in F-BAR cortical level, membrane height (by SRIC), and the rate of membrane height changes (dh/dt) at a fixed location during traveling wave propagation.