On the role of mechanosensitive binding dynamics in the pattern formation of active surfaces

The actin cortex of an animal cell is a thin polymeric layer attached to the inner side of the plasma membrane. It plays a key role in shape regulation and pattern formation on the cellular and tissue scale and, in particular, generates the contractile ring during cell division. Experimental studies showed that the cortex is fluid-like but highly viscous on long time scales with a mechanics that is sensitively regulated by active and passive cross-linker molecules that tune active stress and shear viscosity. Here, we use an established minimal model of active surface dynamics of the cell cortex supplemented with the experimentally motivated feature of mechanosensitivity in cross-linker binding dynamics. Performing linear stability analysis and computer simulations, we show that cross-linker mechanosensitivity significantly enhances the versatility of pattern formation and enables self-organized formation of stable contractile rings.


I. INTRODUCTION
The actin cortex is a thin active polymeric layer that is formed right underneath the plasma membrane of animal cells [1][2][3] . It is tethered to the plasma membrane and serves as a mechanical shield as well as an active regulator of cell and tissue shape [1][2][3] . In particular, the actin cortex was shown to play an important role in cell division, cell migration and embryogenesis (e.g. during epithelial folding) [1][2][3] . Therefore, force generation and pattern formation in the actin cortex are an important prerequisite for the successful accomplishment of cellular function in our body. Previous theoretical and experimental studies suggest that corresponding biological pattern formation and shape dynamics of cortical surfaces could, at least in part, be driven through self-organized dynamic processes of cortical constituents, see e.g. [4][5][6][7][8][9][10][11][12][13] . The actin cortex is scaffolded by a network of actin polymers decorated by a multitude of associated proteins 1,2 . In particular, active motor proteins, such as myosin II, cross-link actin polymers and exert tensile forces on them 1,2 . In this way, motor proteins act as active force dipoles generating active contractile prestress in the cortical layer. Also (passive) actin cross-linker proteins such as α-actinin were shown to contribute to the generation of active cortical prestress 3, [14][15][16][17] . Furthermore, the degree of cross-linking of actin networks was shown to sensitively tune the network's mechanical resistance towards shear deformations 18,19 . In recent years, several experimental studies have provided direct or indirect evidence that actin cross-linkers exhibit mechanosensitive binding dynamics [20][21][22][23][24][25] . In particular, two scenarios have been suggested for cross-linker binding to actin 23,24 ; i) tension-induced increase of bond lifetimes -also called catch bond behavior and ii) tension-induced decrease of bond lifetimes -also called slip bond behavior. Therefore, in a cortex model that captures mechanosensitivity of stress-regulating cortical molecules, the unbinding rate of such molecules needs to depend on the local mechanical stress in the cortical surface. This relationship introduces an additional mechanochemical feedback into the system. Cortical constituents were observed to undergo dynamic turnover on time scales of seconds to tens of seconds [1][2][3] . In accordance, with that observation, cortical mechanics was shown to be time scale dependent and to be dominantly fluid-like on time scales larger than minutes 15,26 . Accordingly, self-organization and pattern formation in the actin cortex has been modeled as the dynamics of an active viscous fluid, see e.g. [4][5][6][7][8]10,11 . Experimental results suggest that the two-dimensional viscosity of the cortical fluid layer is on the order of 10 4 − 10 5 µm · Pa s on long time scales (estimated as the product of parameters K h and τ max from references 15,27 ). By comparison, viscosities of the cytoplasm have been estimated in the range of 10 −2 − 10 Pa s 28,29 . The ratio between cortical and cytoplasmic viscosity is therefore expected to be on the order of 10 3 µm or larger. Hence, the hydrodynamic length scale L h of a) Electronic mail: Corresponding author: sebastian.aland@math.tu-freiberg.de b) Electronic mail: Corresponding author: elisabeth.fischer-friedrich@tu-dresden.de the coupling between the cytoplasmic bulk and the cortex is expected to be orders of magnitudes larger than typical diameters of animal cells. In previous research, Mietke et al. suggested a minimal model of active spherical cortex surfaces 5 . There, the authors reported the ability of the active model cortex to spontaneously polarize and form patterns. However, spontaneous formation of contractile ring conformations were only reported for the case of small hydrodynamic length scale L h = η b /η, where η b is the surface bulk viscosity of the active surface and η is the viscosity of the enclosed passive bulk fluid. Here, we reconsider pattern formation of this minimal model of active spherical fluid surfaces for the case of large hydrodynamic length scale L h taking mechanosensitive binding of cortical molecules into consideration as a new feature.

II. CONSTITUTIVE EQUATIONS FOR AN ACTIVE GEL ON A CURVED SURFACE WITH MECHANOSENSITIVE MOLECULAR REGULATORS
In the following, we will write down the hydrodynamic equations that govern surface deformation and material flow for a closed active fluid surface in a regime of overdamped dynamics as previously introduced in 5 . The two-dimensional domain of the active surface is denoted as Γ. Let (s 1 , s 2 ) be local coordinates of the manifold Γ. Correspondingly, Γ can be locally parameterized as surface position vectors X(s 1 , s 2 ) ∈ R 3 . An induced basis of the tangent space is given by the vectors e i = ∂ i X, i = 1, 2. The surface normal n is n = e 1 × e 2 /|e 1 × e 2 |. The tangential velocity of a surface element in Γ can then be written as v = v i e i . The normal surface velocity will be denoted as v ⊥ = v n n. The curvature and metric tensor are defined as C ij = −n · ∂ i e j and g ij = e i · e j , respectively. The strain rate tensor of the active surface is v ij = 1 2 (∇ i v j + ∇ j v i ) + C ij v n and the in-plane stress tensor reads where η s and η b are the area shear viscosity and area bulk viscosity of the active fluid and ξf (c)g ij is an active stress contribution that is tuned by the contractility parameter ξ and a function f (c) that grows monotonically as a function of the local concentration c of molecular regulators. We have assumed that surface bending makes negligible in-plane stress contributions using the established membrane approximation of thin closed shells 30 . Furthermore, we did not take into consideration contributions of curvature to active in-plane surface stress. This choice was motivated by experimental findings showing that cell surface tension is largely independent of the degree of cell squishing in a cell confinement assay in spite of increasing curvature of the unsupported cell surface 31 . The force balance equations in tangential and normal direction are given by (2) External forces in the above Eq. (2) are assumed to stem from stresses in the inner bulk fluid occupying the space enclosed by the active surface Γ. Therefore, f ext = −n · σ| Γ , where σ = η ∇u + (∇u) T − p1 is the stress tensor, η is the viscosity and u is the velocity field in the bulk fluid. The Stokes Equation demands Furthermore, we impose bulk fluid incompressibility, i.e. ∇ · u = 0 and no slip boundary conditions at the interface with the active surface, i.e.
The surface concentration field c of the molecular regulator species is anticipated to undergo a dynamic evolution due to binding and unbinding events that give rise to an exchange with a pool of free molecules in the inner (cytoplasmic) bulk volume Ω. We anticipate that this inner molecular pool is well-stirred due to fast bulk diffusion. By imposing particle conservation, the concentration in the inner bulk is therefore given by c f ree = N tot − Γ c dΓ / Ω dΩ, where N tot is the overall number of molecules in the system. The time evolution of the concentration field c is modeled as where D is the diffusion constant in the active surface, k on is a binding rate of free molecules and k of f (t i i ) is an unbinding rate of surface-associated molecules. Furthermore, the term ∇ i (cv i ) accounts for concentration changes due to tangential advective flux on the surface and the term C k k v n c captures concentration changes by dilution/condensation through normal surface movement in the presence of curvature. Motivated by experimental observations, see e.g. [21][22][23][24] , we anticipate that the unbinding rate of molecular regulators of active stress depends on the local trace of the in-plane stress t i i . In particular, we make the choice of a Bell model dependence where σ 0 = 2ξf (c 0 ) is the trace of the stress in the steady state (c = c 0 ) and σ c is a characteristic stress that characterizes the nature of the load dependence of the molecular regulator. A decrease of the unbinding rate with stress, i.e. σ c > 0, is known as catch bond behavior while the increase of k of f in dependence of stress, i.e. σ c < 0, is known as slip bond behavior 32,33 .
In the case of the presence of more than one molecular regulator species, the above description of active hydrodynamics can be generalized in a straightforward manner using an active stress contribution ξf (c), where f (c) is a function of a vector of concentration fields with one vector component for each molecular regulator. Furthermore, time evolution equations in the form of Eq. (4) must be defined for the concentration fields of all molecular regulator species.

III. PATTERN FORMATION OF AN ACTIVE SPHERICAL SURFACE WITH MECHANOSENSITIVE MOLECULAR REGULATORS
In the following, we will discuss self-organized pattern formation and surface deformation for the active hydrodynamics as described by Eqn. (2) and (4) emerging from a spatially uniform steady state of a spherical shape with radius R 0 . Therefore, initially, the uniform surface concentration is c 0 = 3N tot k on /(4πR 2 0 ) × 1/(3k on + R 0 k 0 of f ) and velocities vanish on the surface and in the bulk. To identify the stability of this steady state, we perturb the concentration field and the radius of the active surface. As spherical harmonics Y m (θ, ϕ) ( = 0, 1, . . . , ∞, m = − , . . . , ) are an orthonormal basis of the function space on the sphere, we can write δc = ,m δc m Y m (θ, ϕ) and δR = ,m δR m Y m (θ, ϕ) where θ and ϕ denote the polar and azimuthal angle of the sphere, respectively. Correspondingly, we expand the emerging tangential velocity field in the form δv = ,m (δv are the tangential vector spherical harmonic functions. Eqn. (2), (3) and (4) are then expanded in the perturbations δc, δR and δv to linear order and the growth rates λ of each harmonic mode are determined as previously described by Mietke et al. 5 (see Supplementary Section 1). In the linear regime, mode amplitudes δc m , δR m and δv m change in time proportional to exp(λ t), i.e. if λ > 0, the perturbation mode grows and if λ < 0 the perturbation mode shrinks over time. (Note that the growth rate is independent of mode parameter m 5 .) For the limit case of vanishing bulk viscosity η (i.e. diverging hydrodynamic length scale L h = η b /η), an analytical expression for the growth rate is given in Supplementary Section 2. Furthermore, in this limit case, we find that the induced change of the trace of the stress tensor δt i i is to first order for harmonic perturbations with mode degree . Here, λ is the corresponding growth rate. Accordingly, the trace of the stress tensor does not change through the dynamics for modes with = 1 (up to linear order). This behavior results from the lack of area shear contributions for = 1 perturbations for vanishing bulk viscosity (see Supplementary Section 3). We conclude that the unbinding rate k of f (t i i ) stays constant for these modes (in space and time) and the associated mechanosensitivity of the molecular regulator does not affect the binding dynamics (up to linear order). Correspondingly, the stress-amplified binding as mediated by a catch-bond mechanosensitivity of a molecular regulator can be used to specifically amplify modes with > 1.

A. One mechanosensitive regulator species
Using linear stability analysis, we first considered the stability of harmonic modes with ≥ 1 for one molecular regulator species in the limit of large hydrodynamic length scale L h with the specific choice f (c) = 2c 2 /(c 2 0 + c 2 ). For both, catch and slip bond mechanosensitivity, we find that the steady state becomes unstable if the contractility parameter ξ is sufficiently large (see Fig. 1a,b). For slip bond mechanosensitivity (σ c < 0), the mode with = 1 is predicted to be always dominant, i.e. its growth rate is always the largest (see Fig. 1a). This behavior is analogous to the situation of no mechanosensitivity 5 . For catch bond mechanosensitivity (σ c < 0), we find that the versatility of predicted patterns is greatly enriched in that parameter regimes with different dominant modes exist. In particular, parameters can be chosen such that arbitrary -modes can become dominant (see Fig. 1b,c). A mode is denoted as dominant if the real part of the growth rate is maximal at its corresponding degree . To examine the behaviour of the system away from the linear stability region, we solved the Eqn. (2),(3) and (4) numerically in an axisymmetric scheme using a finite-element framework and the Arbitrary-Lagrangian-Eulerian (ALE) method (see Supplementary Section 5) [34][35][36] . We find that numerical simulations agree with predictions of the linear stability analysis with regards to stability of the uniform state (see green and red dots in Fig. 1a,b). Furthermore, numerical simulations reproduce the amplification of different harmonic modes from random initial perturbations of the steady state as predicted by the linear stability analysis for catch bond molecular regulators (see Fig. 1b,d,e,f). B. Two molecular regulator species with opposite mechanosensitivity Some actin cross-linkers such as α-actinin-4 and myosin have been shown to exhibit a catch-bond behavior 21,23,37 . For other actin cross-linkers such as fascin, previous data indicate a slip bond behavior 38 . Therefore, it is a plausible scenario that molecular regulators of cortical activity may exhibit both kinds of mechanosensitivity. In fact, admitting for at least two molecular regulator species with equal but opposite mechanosensitivity, the richness of pattern forming behaviors can be further extended; in particular, it is possible to find parameter regimes, where growth rates exhibit a positive real-valued part and a non-vanishing imaginary part indicating oscillatory behavior.
In the following, we will consider a concentration field c 1 describing a catch bond molecular regulator (σ 1 > 0) in combination with a second concentration field c 2 describing a slip bond molecular regulator (σ 2 < 0). Making the specific choice of f (c 1 , c 2 ) = 2c 1 c 2 /( √ c 0,1 c 0,2 (c 1 + c 2 )) and k 0 of f,1 = 10k 0 of f,2 , we present a phase diagram as predicted by linear stability analysis in Fig. 2a (vanishing bulk viscosity) and 2b (small bulk viscosity ηR 0 /η b = 0.1). For > 1, oscillatory dynamics of mode perturbations can be found (see Fig. 2d,e,f,g). However, the real part of the growth rate of the = 1 mode is always greater than the growth rate of an oscillatory higher order mode. In accordance with this, we see stationary solutions with one concentration maximum at the pole emerging at longer times (see Fig.  2d,e). It is noteworthy that for flat active surfaces, there are no pattern-forming eigenmodes that avoid shear strains in the linear regime. Correspondingly, persistent oscillatory concentration patterns are predicted by the linear stability analysis (see Supplementary Section 4). Making the choice that where σ 0 = 2ξf (c 0 ) is the trace of the stress in steady state (see Supplementary Section 4). The emergence of oscillations in this parameter regime can be understood as follows: perturbative local enrichments of molecular regulators will give rise to local contractile advection because of locally increased f (c). In the presence of shear deformation, the induced flow will give rise to a local increase of the trace of the stress tensor t i i . In turn, the catch bonding property of the first molecular regulator will at first locally increase c 1 and local contractility. After some delay (k 0 of f,1 > k 0 of f,2 !), the locally increased k of f,2 of the slip bonding regulator kicks in reducing the local c 2 and correspondingly also local contractility. Therefore, concentration bumps are first amplified (due to the fast catch bonding regulator) and then shrink with a delay (due to the slow slip bonding regulator) giving rise to oscillatory concentration changes. It is noteworthy that this effect relies on the presence of a non-vanishing shear viscosity and shear deformation. Correspondingly, oscillatory eigenvalues do not occur for η s = 0. Furthermore, it is noteworthy that the here reported mechanism of oscillatory pattern formation relies on the well-known mechanism of a fast activator (catch-bonding regulator) and a slow inhibitor (slip-bonding regulator) as was previously in depth characterized for reaction-diffusion systems 39,40 .

IV. RING CONSTRICTION ON ACTIVE SPHERICAL SURFACES DEPENDS ON THE RATIO BETWEEN SHEAR AND BULK VISCOSITY
During the division of animal cells, an initially round cell constricts through a contractile ring that forms as part of the actomyosin cortex of the cell 41,42 . This contractile ring is known to be enriched in molecular contractility regulators (myosin motors) 41,42 . It has been proposed that mitotic cell constriction might be driven by self-organized pattern formation of the active actomyosin cortex at the cell periphery, see e.g. [4][5][6][7][8][9][10][11][12][13] . Furthermore, myosin motors have been previously suggested to exhibit catch bond mechanosensitivity in their actin binding 21,25 . Therefore, we asked to which extent active spherical surfaces in our model can form a ring constriction by one molecular regulator with catch bond mechanosensitivity. Overall, in the unstable parameter regime, small ringshaped concentration perturbations (proportional to a negative = 2 mode, see e.g. Fig. 1e) lead to a surface constriction that is enhanced by the contractility parameter ξ, see 34 . Here, we show that constriction increases also in dependence of the ratio η s /η b between area shear viscosity and area bulk viscosity (see Fig. 3a-d). Intuitively, this can be understood as follows; in a ring-shaped concentration peak, there is a continuous flux from the poles to the ring 34 , leading to a unidirectional compression of material in the direction orthogonal to the ring. If shear viscosity is high, shear deformations (i.e. aspect ratio changes of surface elements) tend to be avoided by the dynamics of the system. Correspondingly, unidirectional compression of the material orthogonal to the ring tends to trigger compression along the ring. Therefore, the ring gets smaller in diameter, i.e. it constricts.
For growing values of η s /η b , harmonic modes with increasing value of become dominant. In particular, for the parameter choice in Fig. 3, mode dominance transitions from the = 2 mode to the = 3 mode once η s /η b > 0.66. Hence, we see that large symmetric ring constrictions are not a stationary state but start to evolve to an asymmetric distribution aftert = 0.3 or greater. We had seen that in certain parameter regimes ring-shaped concentration distributions may emerge from random concentration perturbations of catch bonding cross-linkers (see Fig. 1b,e). However, we observed in simulations that, in these cases, the emerging ring is not stable in the nonlinear regime and slips to one of the poles such that, eventually, the system emerges to a steady state conformation with two concentration maxima at the poles (see Fig. 3e,f and Movie 1). We conclude, that with the current model of one mechanosensitive catch bond cross-linker in the limit of low cytoplasmic viscosity (large hydrodynamic length scale), we could not find parameters for which we obtained stable constricted rings emerging from random perturbations in a self-organized manner.

V. CONCENTRATION-DEPENDENT SHEAR VISCOSITIES CAN ENHANCE THE ROBUSTNESS OF RING FORMATION
In the search for mechanisms that would stabilize the 'ring' pattern on an active surface, we further investigated modifications of the model that exclusively affected the nonlinear range of concentration and shape perturbations. In particular, we decided to investigate the case of concentration dependent shear viscosities. The concentration dependence of shear viscosity can be motivated by two arguments: 1) It has been observed experimentally that the complex shear modulus increases upon increased cross-linking, see e.g. 18,19 . (In particular, cross-linking was associated to a disproportionate increase of shear stiffness over bulk stiffness as can be inferred from the observation of a crosslinking induced decrease of the Poisson ratio 43 .) 2) The second argument is provided by the experimental finding that active stress leads to shear-stiffening of biopolymeric networks as was suggested by in vivo and in vitro studies 15,[44][45][46][47] . This stress-stiffening might result from the shear-induced nematic alignment of cortical filaments 48,49 or from other stress-induced structural changes 50 . (As active stress in the active surface increases monotonically as a function of concentration of molecular regulators, the shear viscosity increase in dependence of active stress can be translated into a viscosity increase in dependence of molecular concentration c.) We make the specific choice of η s = η 0 s (p + 1)c 2 /(pc 2 0 + c 2 ), where p is a parameter that tunes the concentration dependence of shear viscosity. Here, p = 0 corresponds to the case of no concentration dependence and p = ∞ corresponds to the strong concentration dependence η s = η 0 s c 2 /c 2 0 (see Fig. 4a). We find that increasing values of p lead to increasing constriction of ring distributions in our simulations (see Fig. 4b,c). This is expected from our finding that high ratios of area shear and bulk viscosity η s /η b promote ring constriction, as this viscosity ratio grows with parameter p inside the ring (see Fig. 3a and Fig. 4a,c). While for low values of p the contractile ring slips toward the poles at later times (see Fig. 4e, top panel and Movie 2), we find that for a sufficiently strong rise of η s with concentration c, this ring slippage does not take place anymore (see Fig. 4d,e). In fact, we find that initially off-centered ring conformations can be centred over time by the dynamics of the system (Fig. 4e, Fig. S1 and Movie 3 and 4). This behavior can be understood as follows; in the ring, regulator concentrations are high which increases viscosity locally and slows ring dynamics; when the ring is in an off-center position, the neighboring larger segment of the sphere is drained to a lesser extent of molecular regulators through advective flux into the ring. Therefore, in the larger segment, the regulator concentration starts to become higher than in the smaller segment of the sphere on the other side of the ring (see e.g. time pointt = 0.6 in Fig. 4e, mid  panel). Therefore, the ring is 'fed' more from the side of the larger segment. In turn, the ring moves towards the side of the larger segment over time. If the ring motion is not slow enough, this phenomenon can lead to an overshoot, that means a ring movement beyond the center (see e.g. time pointt = 0.9 in Fig. 4e, mid panel). A damped left-right oscillation of the ring follows until the ring stays in the center (see Fig. 4e, mid panel, and Movie 3). For sufficiently large values of p, viscosities in the ring are so high, that no overshoot happens and the ring centers more quickly (see Fig. 4e, bottom panel, Fig. S1 and Movie 4). In spite of the described effect of ring centering, we observed that the emergence of a ring conformation depends on the nature of the random perturbation. Choosing ten different random seeds, we see ring emergence in seven out of ten cases. In the remaining three cases, concentration aggregates start to emerge at the pole or close to the pole. In these cases, two concentration maxima at the pole are forming as steady state conformation. Therefore, we conclude that emergence of a ring conformation from a uniform state requires perturbations that favor ring-shaped concentration patterns over polar concentration peaks.

VI. DISCUSSION
In this study, we revisited a minimal model of active cortical surfaces of animal cells as introduced by Mietke et al. 5 . Model cortices with spherical reference geometry were described as closed active viscous surfaces whose mechanics is determined by a shear and bulk surface viscosity η s and η b as well as an active stress contribution adjusted by the local concentration c of molecular regulators and an activity tuning parameter ξ (see Eq. (1)). Furthermore, an internal bulk fluid of significantly lower viscosity η was mechanically coupled to the active surface. The corresponding hydrodynamic length scale L h = η b /η was chosen to be large as suggested by measurements on cells 15,27 . The model by Mietke et al. was extended by the new feature of mechanosensitive molecular regulators. In particular, we discussed the case of a molecular unbinding rate k of f with an exponential dependence on the trace t i i of the local surface stress tensor (see Eq. (1)). Corresponding exponential growth of the unbinding rate was mimicking a slip bond behavior, while an exponential decline with t i i was mimicking a catch bond behavior of the molecular regulator. Therefore, the lifetime of catch bond regulators at the active surface is increased for higher values of t i i , while the lifetime of slip bond cross-linkers at the active surface is reduced at increased t i i . Using linear stability analysis and finite element numerical simulations, we analysed the pattern formation of our model in different parameters regimes. For the case of one molecular species with slip-bond mechanosensitivity, we find that the harmonic mode with = 1 always becomes unstable first analogous to the case of stress-independent k of f (see Fig. 1a) 5 . By contrast, for the case of one molecular regulator with catch-bond mechanosensitivity, the spectrum of emerging patterns significantly widens; for arbitrary (spherical) harmonic modes, we found parameter regions in which associated concentration perturbations are dominantly amplified (Fig. 1b). This effect originates from the growth rates of higher order modes ( > 1) being increased through the stress-dependence of k of f in the catch bond case. On the other hand, the growth rate of the = 1 mode remains unchanged because the stress trace t i i stays uniform across the sphere for = 1 perturbations (in the linear regime). In fact, only higher order modes generate area shear deformation in the linear regime, which allows for local variation of t i i and a corresponding variation of molecular lifetime 1/k of f (t i i ) at the cortex (see Eq. (6) and Supplementary Section 3). Admitting for two molecular regulators with mixed mechanosensitivity, we found that oscillatory higher order modes ( > 1) can be generated as transient patterns on the sphere (Fig. 3). For a flat two-dimensional active surface, stable oscillatory solutions are predicted by linear stability analysis for adequate parameter regimes (see Supplementary Section 4). Motivated by the myosin-rich contractile ring that constricts cells during cell division, we asked to which extent active spherical surfaces in our model can be constricted by one molecular regulator with catch bond mechanosensitivity. We observed that amplification of = 2 mode perturbations in ring conformation can give rise to self-organized central accumulation of molecular regulators and central constriction of the cell. We found that this constriction increases with the ratio η s /η b between area shear and bulk viscosity. However, we further observed that corresponding ring conformations of the cortex are not stable with regards to small asymmetric concentration perturbations in the nonlinear regime; over time, ring conformations transition into other patterns, such as two symmetric concentration maxima at the poles, via asymmetric ring slippage (see Fig. 3e,f). To investigate if 'ring' conformations of the active surface could be stabilized in the nonlinear regime, we further tested a modification of the model that exclusively affects the nonlinear range of concentration and shape perturbations. In particular, we investigated the case of area shear viscosities that depend on the local concentration of molecular regulators as motivated by experimental observations 18 . We found that the 'ring' conformation of the active cortex can be stabilised by this effect and even ring centering can be achieved over time (Fig. 4d,e). Numerous experimental findings have highlighted the key role of actin cross-linkers as key regulators of actin cortex mechanics in cells. In addition, experimental evidence is increasing that actin cross-linker binding dynamics is by itself mechanosensitive [20][21][22][23][24] . Here, we investigated the influence of mechanosensitive molecular regulators on pattern formation and deformation of active spherical surfaces. Our analysis shows that mechanosensitivity of molecular regulators significantly enriches the pattern spectrum of mechanochemical self-organization in the regime of large hydrodynamic length scales. Thus molecular mechanosensitivity of cortical regulators may serve cells as a mechanism to tune cortical pattern formation adapted to cellular function.

FIG. 2.
Self-organized pattern formation of an active spherical surface with catch and slip bond molecular regulators of opposite mechanosensitivity. a,b) Phase diagrams of the linear stability of an active spherical surface with catch and slip bond molecular regulators in dependence of dimensionless unbinding rate k 0 of f R 2 0 /D and dimensionless active contractility ξR 2 0 /(η b D) for two different values of internal viscosity (a: ηR0/η b = 0, b: ηR0/η b = 0.1). Blue region: the uniform steady state is stable, grey regions: the uniform steady state is unstable and the eigenmode with degree = 1 is dominant. In the light grey area, also the eigenmode with = 2 is unstable but not dominant. The striped region further indicates a parameter region where the imaginary part of the = 2 mode growth rate does not vanish. Therefore, in the dashed parameter regime, oscillations of concentration perturbations proportional to the = 2 mode are predicted. c) Real (red) and imaginary (blue) part of dimensionless growth rates of different pattern modes in dependence of mode number according to the linear stability analysis. (For one of the eigenvalues, the real part was below −300 for all and is not shown.) Parameters are equivalent to parameters of the simulation presented in panel d and e. d,e) Concentration patterns emerging from a small = 2 mode perturbation (amplitude −5 · 10 −3 ). Note that the depicted concentration range grows over time. For increasing times, the pattern switches between an aggregation at the center (ring conformation) to aggregations at both poles. Once numerical errors build up, there is a spontaneous left-right symmetry breaking and the system is driven towards the dominant = 1 mode. f,g) Concentration perturbations at the surface center as predicted by linear stability analysis (blue line) and as obtained from simulations (orange dots). When the system starts to evolve towards the = 1 mode due to spontaneous symmetry breaking, there is an emerging discrepancy between linear stability analysis and simulations (t > 0.05). Parameters

I. LINEAR STABILITY ANALYSIS
A derivation of the linear stability analysis of the active fluid model on a sphere has been previously presented by Mietke et al. 1 . For our study, we restrict this previously presented analysis to the experimentally motivated parameter regime for a cellular cortex. In turn, we will anticipate i) a negligible external bulk viscosity (η = 0):η models the viscosity of the aqueous medium around the cell while cortical and cytoplasmic viscosities were estimated to be orders of magnitude higher 2-5 .
ii) a negligible cellular surface tension in the absence of molecular regulators (γ = 0, see 1 ): cellular surface tension was shown to be strongly dependent on the activity of molecular motors 2,6,7 .
iii) a negligible contribution of curvature to the active in-plane stress in the surface (ξ = 0, see 1 ): cell surface tension was shown to be largely independent of the degree of cell confinement induced by cell squishing with the cantilever of an atomic force microscope in spite of growing curvature of the free-standing cell surface 6 .
iv) a negligible bending stiffness (κ = 0): This approximation corresponds to the well-established membrane approximation for thin closed shells 8 .
In the following, we discuss the linear stability of the homogeneous stationary state in which the active surface is given by a sphere of radius R 0 , the concentration of molecular regulators is c 0 on the sphere and there are no flows, i.e. v = 0. Following Mietke et al. 1 , we express the shape, concentration and flow perturbations as where Y m (θ, ϕ), Ψ ( m) (θ, ϕ) and Φ ( m) (θ, ϕ) are scalar and vector spherical harmonics, respectively. Including mechanosensitive unbinding of molecular regulators, one finds from the linear stability analysis of dynamic Eqn. (2)- (5), main text, the following system of linear equations in the perturbation amplitudes δc m , δv m and δR m η s Here, Eqn. (2) and (3) are derived from the balance of tangential forces in the surface, while Eq. (4) emerges from the balance of forces normal to the surface 1 . Furthermore, the linearized time evolution of the concentration field gives where σ 0 = t i i,0 = 2ξf (c 0 ) is the trace of the stress in the steady state (c = c 0 ) and the perturbation of the trace of the in-plane stress δt i i is given by Anticipating that perturbations either grow or shrink exponentially in time with a growth rate λ m , we can replace time derivatives δṘ m and δċ m in the above Eqn.

II. THE GROWTH RATE IN THE LIMIT OF VANISHING BULK VISCOSITY
For one molecular regulator, one finds for the growth rate λ of a harmonic perturbation with mode number > 0 (independent of the value of m) in the regime of vanishing bulk viscosity where σ 0 = 2ξf (c 0 ), k 0 of f and c 0 are the trace of the stress, the unbinding rate and the concentration in the steady state, respectively. For = 0, we find

III. THE STRAIN RATE TENSOR IN THE LIMIT OF VANISHING BULK VISCOSITY
For vanishing bulk viscosity (limit of diverging hydrodynamic length scale), the strain rate tensor v i j reads in spherical coordinates (θ, ϕ) to first order where we have omitted the argument (θ, ϕ) of the spherical harmonic function Y m (θ, ϕ) for brevity. Correspondingly, the strain rate tensor is proportional to the identity matrix for = 1, i.e. that the shear contribution to the strain rate vanishes for = 1 to first order.

IV. OPPOSITE MECHANOSENSITIVITY OF TWO MOLECULAR REGULATORS IN FLAT TWO-DIMENSIONAL SPACE
Consider a flat, thin active surface with two molecular regulators of surface contractility described by concentration fields c 1 and c 2 . Let c 1,f ree and c 2,f ree denote the concentration fields of freely diffusing unbound motor and passive cross-linker proteins. Correspondingly, we have the following reaction dynamics As there is no curvature, the in-plane mechanical stress t is where c = (c 1 , c 2 ). Consider a perturbation in the form of a Fourier mode with wave vector k, i.e. δc 1/2 = δc k 1/2 exp(ik · r), δv = δv k exp(ik · r) with k > 0 and r = (x, y). By the force balance equation ∇ · t = 0, we have Therefore, the velocity contribution orthogonal to k has to vanish. Further and, correspondingly, we have The dynamic equation for the perturbation is where σ 0 = 2ξf (c 0 ) is the trace of the stress in steady state. In the following, we will omit the argument σ 0 of k of f,1/2 and k of f,1/2 for brevity. Using Eqn. 12 and 13 in Eqn. (14) and (15), we find finally Let's anticipate that D = D 1 = D 2 , ∂ c1 f (c 0 ) = ∂ c2 f (c 0 ) and c 0,1 = c 0,2 . In this case, we find where χ = c 0,1 ξ∂ c1 f (c 0 )/(η s + η b ). The corresponding eigenvalues of the time derivative operator are Let's consider the particular case of a catch bonding molecular regulator described by c 1 , i.e. k of f,1 < 0 and slip bonding molecular regulator described by c 2 with k of f,2 = −k of f,1 . In this case, where ∆k of f,1 = (k of f,1 −k of f,2 ). For the case ∆k of f,1 > 0, i.e. if k of f,1 > k of f,2 , we obtain an imaginary contribution in the eigenvalue if and therefore, oscillating solutions are predicted to emerge.

V. NUMERICAL METHOD
To examine the behaviour of the system away from the linear stability region, we solve the Eqn. (2)-(4) in an axisymmetric scheme. A detailed description of the numerical method can be found in 9 . Here, we summarize briefly the main features of the method. The method is based on a finite-element framework, implemented in the AMDiS toolbox 10,11 . Rewriting the equations in axisymmetric form reduces the effective computational domain to two dimensions, see 12 . The viscous fluid surface is embedded into two bulk fluids (intracellular and extracellular) to simplify the solution procedure. Intracellular bulk fluid viscosity is set to R 0 η/η b = 0.1. We also included an extracellular bulk Stokes fluid with viscosity R 0η /η b = 0.001 in the numerical model. The finite viscosity of this external bulk fluid ensured stability of the numerical model, but the value was chosen so small that it did not influence velocity fields on the active surface. By comparison,η was set to zero in the theoretical calculations. The two-dimensional numerical domain is discretized by using a surface grid for the cell surface and a bulk grid for the surrounding fluids. Both grids match at the surface by construction. Grid size is set to 0.08 at the surface and coarser afar. The time stepping is constructed by operator splitting, i.e. hydrodynamics and surface concentration equations are solved subsequently in every time step. The incompressible Stokes equations in the surrounding fluids are discretized by standard Taylor-Hood elements (polynomial basis functions of degree two and one for velocity and pressure, respectively). The hydrodynamic equations for both fluid phases are assembled in a unified way following 13 , leading to a continuous velocity field, but allowing for a discontinuous pressure field across the cell surface. The surface viscous stress is treated explicitly and enters the bulk Stokes equations as a boundary condition. The explicit treatment imposes time step restrictions which are handled by underrelaxation (relaxation factor ω = 10 −3 ). The concentration equation is solved on the surface grid and discretized with piece-wise linear basis functions. A conservative weak form of the concentration equations is employed to ensure exact mass conservation under strong tangential surface flows. To accommodate cellular shape changes and lateral translation of the cell, the discretization is performed on a moving finite-element grid using the Arbitrary-Lagrangian-Eulerian (ALE) method 14 . Normal movement of the surface grid is imposed conforming with the flow, while tangential grid movement is chosen such that grid distortion is limited. The grid of the surrounding fluids follows surface motion by harmonically extending surface grid movement into the fluid phases. A preceding time step study was conducted to ensure time stepping errors are small at the chosen time step size of ∆t ∈ [10 −7 , 10 −5 ] depending on occurring stress magnitude. The evolution of patterns was found robust with respect to time stepping errors. Yet, we find that time step size still somewhat influences the time scale at which the instability develops and changes in patterns occur. Further details about the numerical method can be found in 9 . In our study, simulations were used to verify predictions about the system pattern forming behavior according to the linear stability analysis (Fig. 1, main text). To this end, we were running simulations for different model parameters with initial (dimensionless) concentrations corresponding to the steady state concentration 1 perturbed by an axisymmetric concentration profile with random values drawn from the interval −5 · 10 −3 , 5 · 10 −3 ( Fig. 1, main text). Simulations were carried out for a total simulation time of 1. To decide about whether the system was evolving as stable or unstable (see red and green dots in Fig. 1a,b, main text), we calculated the difference ∆c between the maximum and the minimum concentration on the sphere as a function of time. If at timet = 1 the slope of ∆c was positive, the simulation was classified as unstable, otherwise stable.