Pulsatile contractions and pattern formation in excitable actomyosin cortex

The actin cortex is an active adaptive material, embedded with complex regulatory networks that can sense, generate, and transmit mechanical forces. The cortex exhibits a wide range of dynamic behaviours, from generating pulsatory contractions and travelling waves to forming organised structures. Despite the progress in characterising the biochemical and mechanical components of the actin cortex, the emergent dynamics of this mechanochemical system is poorly understood. Here we develop a reaction-diffusion model for the RhoA signalling network, the upstream regulator for actomyosin assembly and contractility, coupled to an active actomyosin gel, to investigate how the interplay between chemical signalling and mechanical forces regulates stresses and patterns in the cortex. We demonstrate that mechanochemical feedback in the cortex acts to destabilise homogeneous states and robustly generate pulsatile contractions. By tuning active stress in the system, we show that the cortex can generate propagating contraction pulses, form network structures, or exhibit topological turbulence.

Introduction previous theoretical models that are either purely biochemical [17,[27][28][29] or based on the mechanics of active gels [18,33,35,36,38]. Using this mechanochemical model, we ask how pulsatory flows arise from the coupling between a fast diffusing activator (RhoA) and a slow diffusing inhibitor (actomyosin), how mechanochemical feedbacks stabilise contractile instabilities into localised patterns, and how active stresses propagate local contractile signals or drive turbulent dynamics. Our model builds upon recent experimental observations of an activator-inhibitor relationship between RhoA and actomyosin [17], by introducing myosindriven contraction to provide a mechanical feedback. By tuning just two biologically relevant parameters, the basal rate of RhoA production and the magnitude of active contractile stress, our model is able to capture a wide range of dynamic phases observed in the actomyosin networks, from travelling waves and pulsatile contractions to stable contractile networks and turbulent dynamics. We find that mechanochemical feedback acts to destabilise stationary states to robustly generate pulsed contractions. At high contractile activity and low basal rates of RhoA production, stable patterns of actomyosin emerge. Furthermore, the mechanochemical system encodes memory of transient perturbations that allows local signals to be translated into propagating contraction waves or stable patterns.

Mechanochemical feedback generates robust pulsatile contractions
To elucidate the role of mechanochemical feedback in the generation of dynamic behaviours and instabilities in the actin cortex, we first study an ordinary differential equation model for the coupling between RhoA, actomyosin, and mechanical strain ( Fig 1A). Here, for simplicity, we neglect spatial variations in chemical concentrations and study the coupled dynamics of RhoA and actomyosin in a locally homogeneous region of the cortex. Recent experiments on Xenopus oocytes [27] and C. elegans embryos [17] suggest that actomyosin pulsation in the cortex is regulated by the excitable dynamics of RhoA GTPase-the upstream regulator of actomyosin assembly and force production. Local autocatalytic activation of RhoA drives rapid initiation of RhoA pulses, followed by F-actin assembly and myosin recruitment. As actomyosin concentrations increase, F-actin dependent accumulation of the RhoA GTPaseactivating proteins RGA-3/4 terminate the pulse through a delayed negative feedback [17].
Here, for simplicity, we represent F-actin and myosin as a single species-actomyosinwhose properties combine active force production (by myosin) and inhibition of RhoA (by Factin). This representation is justified since the time course of appearance and disappearance of F-actin and myosin intensities measured locally during pulses of RhoA activity in the C. elegans cortex are remarkably similar, and the cortical lifetimes of F-actin and myosin II measured by single molecule imaging are also remarkably similar [17]. Based on experimental observations [17], we make the following assumptions in our model. First, RhoA is activated at a constant basal rate, S. Second, active RhoA promotes further production of active RhoA via an autocatalytic feedback, and also promotes the production of actomyosin. Third, actomyosin (F-actin) promotes local inactivation of RhoA by recruiting GAPs. With these assumptions, the rate of production of RhoA-GTP can be written in terms of the RhoA-GTP concentration r and actomyosin concentration m: where S is the applied stimulus (or basal rate of RhoA production) representing the activity of Rho-GEF, a is the rate of autocatalytic production of RhoA, and g is a negative feedback parameter arising from F-actin-driven accumulation of the GAP RGA-3/4 that inactivates RhoA. The constants r a and r g represent the threshold concentrations of RhoA above which the rate of RhoA production is independent of RhoA concentration. The rate of actomyosin production is given by: where S m is the basal rate of actomyosin production, k a is the actomyosin assembly rate and k d is the disassembly rate. While the model parameters S, S m , k a and k d are directly available from single-molecule data [17], the rest are calibrated by fitting our model to the experimental data for the concentration of RhoA and actomyosin during one contraction pulse ( Fig 1B). The overall architecture of the biochemical circuit, namely autocatalytic production of RhoA and delayed feedback inhibition, is similar to recently proposed models for cortical RhoA dynamics [27,29].
To introduce mechanical feedback into this model, we describe a locally homogeneous region of the cortex as an active viscoelastic material with strain u, with contractile stresses , and contractile phase (S = 0.075 s −1 ). (F) Phase diagram of the system, for varying contractile stress, σ a /E, and applied stimulus, S, with η L /E = 5s. (G) Phase diagram of the system, for varying rate of autocatalytic production of RhoA, a, negative feedback parameter, g. See Table 1 where E is the compressional elastic modulus, η L is the local viscosity of the cytoskeletal element, σ a (> 0) is the maximum active stress arising from actomyosin-driven contractions, and m 0 is the concentration of actomyosin at half-maximum stress. The viscoelastic time-scale, τ = η L /E and relative contractility σ a /E are estimated from available data on actin cortex of C. elegans [39]. Mechanical stress feeds back to the dynamics of both RhoA and actomyosin through conservation of mass; if the size of the system doubles then the concentration of chemical species must halve. Written explicitly, @ u (c(1 + u)) = 0, which implies that @ u c = −c/(1 + u), where c is the chemical concentration. While the contraction of a material in one or two dimensions should lead to an expansion in other dimensions, such deformation effects are not relevant for regulating the concentrations of membrane-bound RhoA and cortical actomyosin. It is likely that local contraction or dilation can drive variations in the thickness of cortical actomyosin, but how magnitude of cortical stress depends on thickness remains poorly understood and likely depends on details that vary across cells [40]. The governing equations for the coupled dynamics of RhoA and actomyosin are given by: In the above equations, actomyosin provides additional feedback to both RhoA and itself through changes in mechanical strain. An increase in local actomyosin concentration induces a local contraction, which in turn increases the actomyosin concentrations. By contrast, when actomyosin concentration decreases, there is local strain relaxation leading to a decrease in actomyosin concentration. Using the parameters calibrated from experiments, we simulated these dynamics numerically [41], for different values of the activity parameters S and σ a , to understand the role of mechanochemical feedback in actin cortex dynamics. The model yields four distinct dynamic behaviours-quiescent, excitable, pulsatile, and contractile, depending on the magnitude of the applied stimulus, S (S1 Fig). For low S, the system displays an excitable behaviour-a single pulse of RhoA, followed by a pulse of actomyosin and contractile strain buildup, before reaching a steady-state ( Fig 1E). This excitable behaviour can be visualised in phase space as a trajectory about the intersecting nullclines ( Fig 1C). For this excitable pulse, we observed a loop about the high (r, m) fixed point before settling to a steady state at the lower fixed point.
As the applied stimulus is increased, other behaviours emerge (Fig 1E and 1F). For a small increase in S, we observe sustained pulsatile contractions, when the RhoA nullcline shifts up (Fig 1D), resulting in a single fixed point at high concentrations and the system is trapped in a limit cycle. At even higher values of S, the limit cycle is unstable and the pulse amplitude decays until the system settles in a contracted state with high actomyosin concentration and strain. Finally, for very low applied stimulus, we observe a quiescent mode; both RhoA and actomyosin steadily increase to a fixed set-point.
Previous work has shown that changes in concentration due to mechanical contraction can trigger oscillations in systems that otherwise remain stationary [42]. While pulsatile contractions are observed in the absence of contractile stress (σ a = 0), we find that including mechanical feedback helps to sustain oscillations over a wider range of the parameter space ( Fig 1F). As actomyosin concentration drops after a pulse, the mechanical strain relaxes, further reducing both actomyosin and RhoA concentrations away from a fixed point, allowing another pulse to occur. These pulsatile states occur above a threshold value for the autocatalytic positive feedback parameter a and for moderate values of the negative feedback parameter g (Fig 1G). When the negative feedback is higher than the positive feedback we observe quiescent behaviour, where any changes in RhoA concentration are quickly slowed down by the negative feedback and brought to its equilibrium value. When a is high compared to g we observed several distinct regions of dynamic behaviours. When both a and g are low, the system settles to a contractile steady-state. When both a and g are high the system is excitable, and when a is high and g is moderate the system becomes pulsatile (Fig 1G).

Spatial propagation of pulsatile flows and contractile pattern formation
To investigate the role of active mechanical stresses and RhoA signalling on spatial patterns and flows within the actin cortex, we develop a continuum active gel model of the actin cortex, coupled to the excitable RhoA signalling network and a viscous cytosol ( Fig 1A). Using the C. elegans embryo as a model system, we develop a one-dimensional description of the actin cortex with periodic boundary conditions, assuming azimuthal symmetry around the long axis of the cell. The actin cortex is modelled as a porous Maxwell viscoelastic material, behaving elastically at short times and remodelling over longer time scales, with flows in the cortex advecting RhoA and actomyosin. The constitutive equation defining the time-evolution of stress field σ (x, t) is given by: where the spatial coordinate x denotes distance along the surface of the cell, τ is the viscoelastic relaxation time scale, η is the actomyosin network viscosity and v(x, t) is the actomyosin flow velocity. The second term on the right hand side of the above equation is the active stress term with σ a the active contractile stress, and m 0 is the actomyosin concentration at half-maximum active stress. Local balance of viscoelastic forces with cytosolic drag, pressure and actomyosingenerated active contractile forces can be written as: where Γ is the frictional drag coefficient between the actomyosin network and the cytosol, v cyto is velocity of the cytosolic fluid, P is the cytosolic pressure, and ϕ is the cytosol volume fraction. Conservation of mass implies, ϕv cyto = −(1 − ϕ)v. Darcy's law applied to the poroviscoelastic continuum yields a relationship between the pressure gradient and relative velocity between the cytosol and the cytoskeletal network [43,44] G�ðv cyto À vÞ ¼ À @ x P : ð8Þ By taking a spatial derivative in Eq (6), and using Eqs (7) and (8) alongwith mass conservation we obtain the following equation for the cytoskeletal velocity field: where γ = Γ(1 + ϕ)/ϕ is the effective drag coefficient. Mechanical feedback is then introduced through myosin-induced F-actin flows which advect both RhoA and myosin [18,45], leading to their local accumulation due to convergent flows or depletion via divergent flows. These advective flows compete with reaction and diffusion of RhoA and actomyosin: where D r and D m are the diffusion coefficients for RhoA and actomyosin, taken from Nishikawa et al [18]. The default parameter for active stress and viscoelastic time scale are taken from Saha et al [39]. The model equations are then numerically integrated [46] in a periodic box of length L = 10λ, where l ¼ ffi ffi ffi ffi ffi ffi ffi Z=g p is the hydrodynamic length scale. As shown in Fig 2A-2C, the numerical solutions predict a wide diversity of dynamic states as the active contractility s 0 a ¼ s a =g and RhoA stimulus S are varied-from stationary patterns to propagating waves and pulsatile flows (S2 Fig). At low s 0 a , diffusion dominates over advection, giving rise to spatially uniform concentration profiles that exhibit excitable, oscillatory, or contractile dynamics (Fig 2C, top row), as observed in the cellular-scale model (Fig 1). In the absence of mechanical feedback, reaction-diffusion alone cannot generate spatial patterns because the activator, RhoA, diffuses much faster than the inhibitor, actomyosin (D r � D m ). In this regime, the homogeneous state becomes unstable and Turing patterns emerge, only when D r � D m (S3 Fig).
As s 0 a is increased, contractile instabilities develop due to local accumulation of actomyosin, allowing finite wavelength patterns to emerge. At low S and high s 0 a , we observe stable localized peaks of actomyosin and RhoA (Fig 2C, bottom left). Initially, autocatalytic positive feedback of RhoA creates small RhoA concentration peaks (Fig 3A, 1st panel), which in turn produce actomyosin (Fig 3A, 2nd panel). As actomyosin begins to accumulate, large inward flows are generated, further increasing both RhoA and actomyosin concentrations (Fig 3A, 3rd panel). Finally, actomyosin-induced inhibition of RhoA results in RhoA localization on either side of the actomyosin peak ( Fig 3A, 4th panel), in contrast to Turing patterns where activators and inhibitors overlap. With no RhoA advection, RhoA diffuses away from the peak, producing actomyosin behind it and generates a travelling wave (S4A Fig). With no actomyosin advection, actomyosin concentrations remain too low to prevent RhoA from diffusing and create a uniform steady state (S4B Fig). When several contractile actomyosin foci exist, they attract each other and merge into a single peak. This phase is reminiscent of equatorial RhoA zones during cell division [47,48], and medial [49] and junctional [50] RhoA domains in polarized epithelia.
At higher S and moderate s 0 a , we observe propagating waves and pulsatile flows (Fig 2C,  bottom). A higher level of excitation in RhoA generates a localized actomyosin peak with high contractility. Away from the actomyosin peak, RhoA concentrations are higher and are advected towards actomyosin (Fig 3B, 1st panel). Advected RhoA produces actomyosin as it moves, such that the newly assembled actomyosin generates flows away from the centre, reducing the actomyosin concentration at the centre (Fig 3A, 2nd panel). Once the initial contraction dissipates, the two remaining actomyosin peaks merge, completing a cycle (Fig 3A,  panels 3-4).
In the propagating waves state, we observe the periodic formation of RhoA pulses, which travel in waves before annihilating as two waves meet, as observed in the starfish oocyte [16,27]. As s 0 a is increased, actomyosin pulses generate large contractile flows that advect neighbouring pulses, creating chaotic, aperiodic motion. The pulses are highly localised, with small regions of high actomyosin concentration that form before dispersing, followed by a new pulse elsewhere, as in the C. elegans embryo [17]. Mechanical feedback through advection is necessary for the waves to form (S4D Fig). Without advection of RhoA however, we may still observe waves, since the actomyosin generated by RhoA forms clusters that travel with RhoA (S4C Fig). However, these waves display much less chaotic motion when compared to the system with advection ( Fig 2C, bottom right).
Linear stability analysis of the continuum model reveals the role of active stress in destabilising the homogeneous state of the system (Fig 2B, S5 Fig). At low S, three fixed points exist (Fig 1B), with the lower fixed point being stable for high s 0 a . As S is increased, the RhoA

PLOS COMPUTATIONAL BIOLOGY
nullcline shifts up until only one fixed point remains and the system enters the pulsatile regime ( Fig 1C). Increasing s 0 a leads to contractile instabilities that manifest as propagating waves and pulsatile flows (Fig 2C). While active stress is required for wave propagation, higher s 0 a leads to lower wave speeds (Fig 2B).

Response to local bursts of contractile activity
While our model can capture many of the dynamic states observed in the actin cortex, cells must have the ability to actively switch between flowing and contractile states during  To understand how RhoA signals can propagate through space and induce state transitions, we locally turned on RhoA activity and examined the output response. Starting from rest with no applied stimulus (S = 0), we applied an increased stimulus at the central region (S = 0.03 s −1 ) (Fig 4A). By changing active stress in the system, we were able to regulate both the ability for the input RhoA to propagate in space, and for the system to remember the spatial location of the signal. We observed three distinct phases (Fig 4A): (i) propagation of a bistable front, where the memory of the signal location is lost and a global increase in RhoA concentration is observed, (ii) a soliton phase with transient spatial memory, and (iii) a high memory phase where a stable RhoA pattern is maintained.
At low active stress (s 0 a ¼ 12:5 μm 2 /s), RhoA is excited into a pulsatile state within the activation region, while spreading laterally through diffusion (Fig 4A, left). A front of highly concentrated RhoA travels away from the source, increasing the total RhoA and actomyosin in the system (Fig 4B). This is reminiscent of a bistable front propagation in classical excitable systems [50,51], where the system switches to the high concentration stationary state ( Fig 2C).
As active stress is increased, we find that mechanical feedback is able to tune the properties of the biochemical system away from classical excitable systems. At moderate active stress (s 0 a ¼ 25 μm 2 /s), RhoA pulses spread out as two solitons before annihilating as they meet ( Fig  4A, middle), as seen in surface contraction waves [19,20]. In contrast to classical excitable systems, actomyosin generated behind the RhoA wavefront increases its own concentration through contractile flows. This region of highly concentrated actomyosin behind the wave acts as a barrier that inhibits RhoA, preventing the system from switching to the high fixed point, and instead creates a soliton. Such a barrier can be seen at s 0 a ¼ 12 μ 2 /s, although it is too weak to prevent diffusion, with the band of reduced RhoA concentration behind the wavefront (Fig 4A, left).
At high active stress (s 0 a ¼ 50 μm 2 /s), the contractile forces within the activated region are strong enough to maintain a spatially localized state that persists after the activation, remaining stationary in a fixed location at the centre of activation (Fig 4A, right). This suggests a potential mechanism for cells to direct the locations of contractility. Without the guidance of externally induced RhoA activation, the system in this parameter regime is incapable of spontaneously forming a stable pattern.
To quantify the input-output relationship of the system, we measured the correlation between the input signal and the output RhoA concentration at long times (Fig 4C). At low s 0 a , RhoA spreads outwards, leading to a loss of correlation between input and the output signal, akin to a memoryless system. For the soliton case, the shape of the input signal is remembered, resulting in a negative correlation as the waves travel away from the source. Finally, in the spatially patterned phase, a high-memory state emerges where RhoA remains localized at the centre of activation, with a strong positive correlation between the input and the output. These results suggest that mechanical stresses can play an important role in biochemical signal propagation, and in retaining the spatial memory of activity. For low active stress, signals propagate the fastest with global changes in contractility and no memory of the spatial location of the signal. As active stress is increased, we observe contraction waves propagating away from the source, displaying transient memory. At higher stresses, a high memory state develops, where transient local RhoA activations create localized contractile states.

Network formation, pulsatile flows and topological turbulence
We now proceed to analyze our model in two spatial dimensions to investigate how contractile patterns form and propagate on the surface of a cell. The constitutive equation defining the time-evolution of the stress field σ(x, t) is given by where v is the velocity vector, η s and η b are the shear and bulk viscosities, and we have assumed that the bulk and shear components have the same relaxation timescale, τ. Using local force balance we find that the equation governing the temporal evolution of the velocity field v is given by We assume that the bulk viscosity is higher than the shear viscosity, such that Z b ¼ 3 4 Z and

PLOS COMPUTATIONAL BIOLOGY
where η is the experimentally measured viscosity [18]. The reaction-diffusion equations for RhoA and actomyosin concentration fields are given by _ r þ r � ðrvÞ ¼ R r ðr; mÞ þ D r r 2 r; ð15Þ In the two-dimensional model, as we vary active stress (σ a ) and the applied stimulus of RhoA (S) we observe similar dynamic phases as in the one-dimensional model, though with more varied spatial organisation. At low applied stimulus and high active stress, we observe the formation of stable contractile networks (Fig 5A and S1 Video), analogous to the pattern formation phase observed in the one-dimensional model. Here, an initial excitable pulse forms many small cluster of highly concentrated actomyosin (Fig 5A, left panel). Over time, the clusters stretch and merge with other nearby clusters (Fig 5A, middle panel), eventually creating a space-spanning network of actomyosin, which gradually condenses into a stable configuration of fully connected edges and vertices (Fig 5A, right panel). These condensates of actomyosin are stabilized by contractile flows, similar to the localized contraction peaks observed in the one-dimensional model (Fig 3A). This sequence of actomyosin patterning is reminiscent of apical microridge formation in some epithelial cells [52].
As the applied stimulus S is increased we observe pulsatile contractions and propagating waves of RhoA and actomyosin. At moderate active stress we observe propagating waves ( Fig  5B and S2 Video), in which waves of actomyosin propagate with constant velocity until two waves collide and annihilate, with new waves periodically produced. Finally, at high applied stimulus and high active stress, we observe the pulsatile contraction phase (Fig 5C and S3  Video). Pulses of actomyosin contract into highly concentrated foci before dispersing and new pulses are formed. These dynamics are distinct from propagating waves as the pulses of actomyosin move more erratically as clusters merge, and show much more spatial variation in concentration. The corresponding velocity and strain fields are shown in S6 Fig. The propagating waves of RhoA show features of topological turbulence, as recently reported in both the Xenopus eggs and embryo [27] and the starfish oocyte [16]. In this regime, RhoA and actomyosin exist with a limit cycle (Fig 1D), allowing us to extract the phase of oscillation using the relation � ¼ tan À 1 rÀ � r mÀ � m À � , where � r and � m are the mean RhoA and actomyosin concentrations (Fig 6A). We find that numerous topological defects exist in the phase field, indicating locations where the phase field is discontinuous, as reported in experiments [16,27].
These may be parameterised by the winding number: the total amount that the phase changes by following an anti-clockwise contour around the defect. From the phase velocity field v ϕ = rϕ (Fig 6B), the integer winding number is calculated as 1 where C is a closed curve around the defect containing no others. In a +1 defect, the phase field spirals anticlockwise about the defect, while in a -1 defect it spirals in a clockwise direction (Fig 6D). In our model, we find approximately the same number of +1 and -1 defects. In addition, oppositely charged defects attract one another, and annihilate when close enough, conserving the overall charge of the system (S4 Video).
The phase velocity also describes an effective velocity field due to a potential flow with point vortices of magnitude ω = r × v ϕ , with +1 and -1 defects corresponding to the center of vorticies with clockwise or anticlockwise rotations. We find a linear scaling between the effective kinetic energy � E ¼ 1 2 hv 2 � i and the effective enstrophy � O ¼ 1 2 ho 2 i, a measure of vorticity, for different values of active stress, in agreement with experimental data on RhoA in starfish oocytes [16]. Furthermore, when we quantify the phase velocity statistics, we find a power law tail in the probability density function which decays like v −3 , likely due to interacting vortices, consistent with experimental results [16] and theoretical predictions for vortex interactions [53]. Thus our model can quantitatively capture the key statistics and scaling laws of defect mediated turbulence in the RhoA phase field.

Discussion
Activator-inhibitor coupling between RhoA and actomyosin has been shown to generate excitable dynamics and oscillations in the concentration of RhoA, F-actin and myosin-II [17,18,27,29]. Purely biochemical feedback models with reaction-diffusion dynamics are thus sufficient to reproduce pulsatile dynamics of cortical RhoA and actomyosin concentrations as observed in experiments [17,18,27]. However, in contrast to Turing pattern-forming systems [26], the activator RhoA is fast-diffusing while the inhibitor actomyosin is slow-diffusing, thus the system should be expected to form spatially uniform states. Stochastic noise is therefore required to trigger the formation of waves in these excitable chemical systems without additional feedback mechanisms [17,27,29].
Mechanical forces generated by myosin-II in F-actin networks play an essential role in regulating the emergent dynamics of the cortical actomyosin networks. Mechanochemical feedback motifs in which mechanical strain locally upregulates actomyosin concentration through a mechanochemical coupling can lead to the formation of stress patterns [32,54]. Alternatively, the contractile stresses generated by actomyosin drives advective flows of RhoA and actomyosin up the concentration gradients, and thus providing an additional mechanical feedback through active stress. Active fluid models involving the diffusion and advection of active stress producing components, are cable of displaying both pulsatile and stable contractions without explicit coupling to signalling networks [33,35]. In particular, it has been suggested that spontaneous pulsatory patterns may emerge when the fast-diffusing chemical species upregulates and the slow-diffusing species down-regulates the active stress, or the up-regulator of active stress turns over faster compared to the down-regulator [35]. These conclusions, however, stand in contrast to the known properties of cortical RhoA and actomyosin, where the both the fast-diffusing RhoA and slow-diffusing actomyosin up-regulate the active stress.
Recent work suggests that feedback between chemical signaling and mechanical forces are essential to regulate pulsatile actomyosin dynamics in the cellular actin cortex [18]. Without coupling to chemical signaling, a positive feedback between contractile stresses generated by actomyosin and advective flows can lead to a density instability in the cortex, leading to the formation of high density myosin regions. By regulating myosin assembly and disassembly by an independent oscillatory signal, the contractile instability can be controlled, generating a pulsatile pattern of myosin-driven contractions. While this model neglects feedback between excitable RhoA and actomyosin, it illustrates how the coupling between signaling and mechanics can lead to spontaneous formation and propagation of RhoA waves without noise. At the cellular scale, dilution effects due to contraction could be capable of inducing oscillations in otherwise stationary systems [42]. However, the role of myosin-induced forces on wave propagation may be system dependent. In the Xenopus embryo, the myosin inhibitor blebbistatin had no significant effect on the observed pulses [27], but in the starfish oocyte, applying blebbistatin prevented surface contraction waves from propagating [20].
In this work, we couple actomyosin mechanics and biochemistry with RhoA signalling in a unified theoretical framework, where mechanical stresses and chemical signalling feedback to each other, regulating the emergent dynamics of the cortex. With this mechanochemical model, we are able to reproduce waves, patterns and pulsatile dynamics observed in the actin cortex by changing only two control parameters, namely the amount of active stress, and the basal rate of RhoA production. We demonstrate that mechanical feedback helps to robustly sustain pulsatile contractions, and provides a new mechanism for spatial patterning of RhoA and myosin, distinct from classical Turing patterns [26]. Having an excitable cortex with tunable dynamics and mechanochemical feedback can be beneficial in a number of morphogenetic contexts. Sustained pulsed contractions can act as a mechanical ratchet to sequentially reduce cell area and drive tissue bending [11-14, 55, 56]. Mechanochemical feedback may be important in enabling robust control of this morphogenetic ratchet. Cells may also display excitable behaviours which are governed by mechanics, for example, with active contractions triggered by large mechanical strains on the cell, helping to prevent tissue rupture [57,58].
At low RhoA activity and high contractile activity, myosin forms stationary patterns, with RhoA localised on either side of the myosin peak. In two spatial dimensions, we observe the formation of many small clusters of high myosin concentrations, which stretch towards neighbouring clusters and condense into a stable contractile domains of actomyosin and RhoA. This is reminiscent of stable patterns observed in vivo on the apical surfaces of epithelial cells, including microridges [52], stable junctional and medial RhoA domains [49,50]. As contractile activity is increased, the waves of myosin become highly peaked, creating pulsing contractile networks, reminiscent of pulsatile contractions in the C. elegans cortex [5]. At moderate activity, we observe periodic waves of myosin which annihilate before new wavefronts are generated in regions with low actomyosin concentration. These waves show features of topological turbulence, reminiscent of topological turbulence in RhoA concentration waves observed in Xenopus and starfish eggs [16,27]. Besides generating patterns and flows, the level of active contractile stress regulates the system's response to local RhoA signals, enabling phase transitions and memory entrainment as the activity is increased. Together, these results highlight the importance of considering both mechanical forces and chemical reactions when modelling the actin cortex, and they reveal a variety of ways in which cells can tune the dynamic coupling between RhoA activity, force production, and advective transport to control morphogenetic behaviours.

Linear stability analysis
To understand how mechanical stress and biochemical stimulus work together to generate pulsatile behaviour and propagating waves, we perform linear stability analysis about the homogenous steady-state (or fixed point) to determine the stability of the state, and the onset of oscillatory instabilities. In particular, we use linear stability analysis to compute the frequency, wave length and wave speed of the wave-like states that emerge from perturbations about quiescent homogeneous state. A perturbation of the formỹe ikx about the steady state, (r, m, 0), results in the following Jacobian: By finding the eigenvalue with the highest real part over all k, λ(k � ), we obtain the fastest growing mode with wavelength 1/2πk � , frequency Im(λ(k)), and a wave speed of Im(λ(k))/2πk � .

Numerical method
The system of ordinary differential equations (Eqs 1-5) are numerically integrated using the SciPy package, which implements a Runge-Kutta 4(5) method. The 1D and 2D systems of partial differential equations are solved numerically using the FiPy python package, a finite volume PDE solver. The equations discretized in a periodic box with edge lengths L = 10λ, where l ¼ ffi ffi ffi ffi ffi ffi ffi Z=g p is the hydrodynamic length scale, using 150 grid points in each dimension. The system is then integrated forward in time using an implicity method, outputting in timesteps of dt = 1s.

Model parameters
Tables 1 and 2 list the parameters used in our simulations. The RhoA-actomyosin chemical reaction model and parameters are obtained by fitting data from Michaux et al [17]. Parameters for dimensionless active stress and viscoelastic time scale are taken from Saha et al [39], and myosin and RhoA diffusion coefficients are taken from Nishikawa et al [18].