Spontaneous formation of chaotic protrusions in a polymerizing active gel layer

The actin cortex is a thin layer of actin filaments and myosin motors beneath the outer membrane of animal cells. It determines the cells’ mechanical properties and forms important morphological structures. Physical descriptions of the cortex as a contractile active gel suggest that these structures can result from dynamic instabilities. However, in these analyses the cortex is described as a two-dimensional layer. Here, we show that the dynamics of the cortex is qualitatively different when gel fluxes in the direction perpendicular to the membrane are taken into account. In particular, an isotropic cortex is then stable for arbitrarily large active stresses. If lateral contractility exceeds vertical contractility, the system can either from protrusions with an apparently chaotic dynamics or a periodic static pattern of protrusions.


Introduction
The cytoskeleton is a dynamic meshwork of filamentous protein assemblies that is of crucial importance for living cells [1][2][3]. The filaments, formed by the proteins actin or tubulin, interact with molecular motors like myosins or kinesins. By using a chemical fuel, these are able to generate mechanical stresses in the filament network making it an active gel [4]. Physical studies of active matter often use a hydrodynamic approach [5]. For the cytoskeleton, this approach has been used to study structures and dynamics that are either alien to conventional matter or behave differently from their passive counterparts. Examples include topological point defects [4,[6][7][8][9][10][11], spontaneous flow transitions [12][13][14], or buckling instabilities of contracting gels [15,16].
Physical tools and concepts have also been used to study structures formed by the cytoskeleton in living cells. In this context, the actin cortex-a thin actomyosin layer beneath the plasma membrane of animal cells-is of particular interest. It determines the mechanical properties of these cells [17,18] and plays an important role in cellular morphogenesis [19]. Flows in the actin cortex also play a role during development, for example, for defining the head and the tail of the nematode worm Caenorhabdtits elegans [20] or for the remodeling of cellcell junctions [21]. In addition, the cortex hosts a number of cytoskeletal structures [22]. A particularly striking example is the contractile actomyosin ring that cleaves animal cells into two daughter cells during division. The ring might be generated by a dynamic instability of the cortex [23,24] and its dynamics when connected to a membrane have been studied [25,26].
The cortex is formed by actin filaments that polymerize directly at the cell membrane [3]. Its thickness, which has been reported to be of a few hundred nanometers [27,28], is determined by the actin polymerization velocity and its disassembly rate [29]. Due to the large aspect ratio, in theoretical analysis, the cortex is typically considered as a two-dimensional surface [26,[30][31][32][33][34][35][36]. Such an approach has notably been used to measure physical properties of the cortex [37,38]. However, a formal reduction of the dynamics to the two lateral dimensions as in a thin-film approximation is not always possible. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
In this work, we explore the dynamics of a layer of a contractile active gel on a rigid substrate with a constant influx of material at the substrate. Contrary to previous work, we also account for gel flows in the dimension perpendicular to the substrate. We find that an isotropic layer is always stable against perturbations. In case active contractile stresses are stronger in the lateral direction than perpendicular to the substrate, we observe an instability leading to a local increase in gel density. This is similar to the dynamics when neglecting the gel thickness [36,39]. However, the states emerging in our system differ fundamentally from those reported in these works: we find periodic stationary states and states that are apparently chaotic. In all cases, they exhibit a characteristic length scale. When neglecting the gel thickness, such states can be obtained only after modifying the dynamics, for example, by introducing nonlinear gel assembly terms [40].

Dynamics of the actin cortex
We consider a three-dimensional viscous active gel, growing into the half space z0 at the surface (x, y, z=0), figure 1. We describe the active gel as a continuum and use a hydrodynamic approach to describe its dynamics [5]. The state of the gel is determined by the density ρ and the velocity field v=(v x , v y , v z ). For the sake of simplicity, we assume invariance along the y-direction and v y =0. Mass conservation then reads: where k d is the gel's disassembly rate. It is on the order of 0.05/s [17]. Treating the cortex as a viscous material implies that we restrict attention to the dynamics on time scales longer than 10 s. To determine the velocity field, we employ momentum conservation. Inertial effects can typically be neglected for the dynamics of the actin cortex. To estimate the Reynolds number of the gel flow, we use the formin-induced actin polymerization rate, which is the fastest process of actin growth and is about 400 nm s −1 [17], as the characteristic velocity. The characteristic length is et by the cortex thickness of 500 nm. For the kinematic viscosity, we take that of water, that is, 10 −6 m 2 s −1 , such that the Reynolds number is 2 × 10 −9 . Momentum conservation is thus expressed by force balance. In absence of external forces, we have s  =0 tot · , where σ tot is the total mechanical stress.
The total mechanical stress can be decomposed into several contributions [4]. As we describe the cortical actin gel as an active viscous fluid, it has a viscous part, a hydrostatic part, and an active part. Whereas the first two parts are well-known from the studies of conventional fluids, the latter results from the action of molecular motors and from actin-filament assembly. These processes transform the chemical energy released in the hydrolysis of adenosine triphosphate (ATP) into mechanical work and thus turn the actin cortex into an active material [5]. In particular, such materials can spontaneously contract in the absence of external compressive forces.
The different parts of the total stress tensor can in turn be decomposed into a diagonal and a traceless part. For the viscous part, we assume that the shear and the bulk viscosity, which respectively characterize the traceless and the diagonal parts of the dissipative stress due to velocity gradients, are equal. For the hydrostatic and the active stresses, the non-diagonal parts are related to gradients in the average orientation of the filaments. On average, the actin filaments point with their plus ends towards the membrane, such that the cortical actin gel is polarized along the z-direction. Furthermore, contraction of the actin cortex likely leads to filament alignment parallel to the membrane, such that there is nematic order of the gel along the x-direction. We will neglect gradients in the polarization along the x-direction as well as gradients in the nematic order along the z-direction, Figure 1. Illustration of an anisotropic active gel polymerizing on a surface. Assembly of the filaments (small yellow circles) at z=0 is mediated by nucleators (big blue circles) of density ρ 0 . This leads to a flux of gel k p δρ 0 perpendicular to the surface, where k p is the polymerization rate and δ the length added to a filament upon addition of a subunit. The gel disassembles in the bulk at a constant rate k d . Active stress is generated by motor proteins (green connected circles).
such that the hydrostatic and the active parts of the stress tensor are diagonal. Furthermore, we assume that the average orientation of the filaments is constant, such that the active and the hydrostatic stresses can be expressed as function of the actin gel density only. We discuss the limitations of our assumptions at the end of our article.
Explicitly, we write the force balance equations as Here, η denotes the viscosity and P x z , are the components of the non-viscous contribution to the total stress. As discussed in the previous paragraph, the latter has two contributions: an effective hydrostatic pressure and the stress generated by active processes in the gel. We write these contributions in a form similar to that chosen in [29]: with b>0. The first term accounts for the activity of the system, where we take > a 0 x z , reflecting the contractile nature of the active stress.
We chose the active term to be proportional to the third power of the density, because for a motor to generate active stresses it has to connect two filaments. If the motor binding and unbinding dynamics is fast, the density of motors bound to filaments is proportional to the filament density. Together this yields an active stress that is proportional to ρ 3 analogously to the law of mass action. However, other choices could have been made. For example, as in a generic expansion of the hydrostatic and active stress in terms of the gel density, we could have started with a linear term in ρ. The essential features are that Π x and Π z decrease with ρ for low enough densities and increase with ρ as r  ¥. As we have checked, however, our results are largely independent of the specific form chosen for P x z , . It suffices that the dependence of P x z , on ρ has a shape similar to that given in equation (4). The dynamic equations are completed with boundary conditions. As mentioned in the introduction, the polymerization of cortical actin is essentially restricted to the direct vicinity of the membrane [3]. To describe the assembly of the cortical actin gel at the membrane, we thus pose ρ(z=0)=ρ 0 and v z (z=0)=v p , where ρ 0 is the density of nucleators at the membrane and v p the polymerization speed. In addition, the tangential stress at the surface is balanced by friction of the gel along the membrane: where ξ is an effective friction constant. On a molecular level it results, for example, from the motion of the actin filaments close to the membrane, which can induce flows of lipids or proteins therein. Furthermore, we impose periodic boundary conditions in the x-direction with period L. In the z-direction, the system extends to infinity and r  0 for  ¥ z . If we ignore the z-direction in equations (1)-(3), we arrive at a dynamic system similar to others that have been used before in the analysis of the actin cortex or other thin active gel layers [36,39]. For high enough activity, these systems typically generate contracted stationary states with a single region of high gel density unless nonlinear assembly terms are considered [40]. However, there is no way to integrate the dynamic equations (1)-(3) over z to get a closed system for the averaged fields Here, h is a length characteristic of the average thickness of the gel layer. For example, integrating equation (2) yields ( ) cannot readily be expressed in terms of r x ( ) and v x ( ), such that further approximations are necessary. For example, a thin-film approximation was introduced in [37], where the gel was assumed to be incompressible. In that case, filament turnover was effectively accounted for by an ad hoc dynamics for the cortical thickness h.
In the following we will use a dimensionless form of the dynamic equations (1)-(4) and scale time by , concentration by ρ 0 , and stress by ηk d . Hence, the relevant parameters are ) . We will use the same notation for the original and the rescaled parameters. For future use, we introduce the anisotropy parameter =  a a x z .

Results
Before we analyze equations (1)-(4), let us briefly recall the properties of the stationary states in case the gel is invariant with respect to lateral translations [29]. Then, two qualitatively different stationary states exist: either the density decays exponentially or it jumps from a finite value ρ e to zero at z=h∼ℓ. For the 'exponential profile', the velocity v z converges to a finite value for  ¥ z . In contrast, v z =0 for z>h for the ʼstep profile'.
The latter exists only for large enough contractility, that is, for a z above a critical value  a 0 c , and when the density of nucleators ρ 0 exceeds a critical value ρ c . Otherwise, the profile decays exponentially. The values of ρ e and ρ c are given by the zeros of the function h P + k 2 d .

Numerical integration of the dynamic equations
We start our analysis of the general case by numerically solving the dynamic equations (1)-(4). In each time step, we first determine the velocity field through the force balance equations (2) and (3), where we use Fourier decomposition along x and finite-differences along z. We then update the density. The contribution of r ¶ v z z ( )is obtained by an up-wind finite-differences scheme in real space. Explicitly, r r r r , where the prime indicates the value at time t+Δt. All other quantities are taken at time t. To improve the stability of the scheme, we have added a diffusion term with diffusion constant D=10 −3 to the mass conservation equation (1). In all simulations, we use Δx=0.004, Δz=0.007, and Δt=0.0005. We checked that our results do not change for smaller values. As initial condition, we take r ) are independent and uniformly distributed random numbers between −0.05 and 0.05. Note, that through equations (2) and (3) the initial condition on the density ρ also fixes the initial velocity field (v x , v z ).
We first consider the limit of zero friction, i.e.
. If a z >a c , then there is a critical anisotropy ò c , such that the step profile is unstable for ò>ò c . In figure 2 and movie1 5 , we present the solution for a z = 7 > a c and ò = 1.18 > ò c . Locally, the gel contracts laterally, which leads to the formation of finger-like protrusions growing into the direction of increasing z. Continuously, neighboring protrusions merge and new ones appear. As can be seen in the kymograph, figure 2, the system does not settle in a periodic state. Instead the dynamics appears to be chaotic.
To further investigate whether the state is chaotic or not, we numerically solved the dynamic equations with slightly different initial conditions. Explicitly, we used r r  figure 2(c). To firmly establish the chaotic nature of the system, we would need to investigate sensitivity to the initial conditions for a dense set of initial conditions. In this work, we refrain from doing so and rather analyze a additional property showing consistency with chaotic behavior, namely the relaxation of the protrusion auto-correlation C(t), figure 2(d). Here, we have defined the auto-correlation function as s t s s t s = á -á ñ , where σ(x, t) is 1 inside a protrusion and 0 outside and averages are performed over time and space.
Our numerical results indicate that the chaotic state presented in figure 2 emerges directly at the critical anisotropy ò c . Such a direct transition to chaotic dynamics, which is impossible for systems with a finite number of degrees of freedom, is known to exist for other spatially extended dynamic systems [41,42]. Often this is linked to a group of continuous symmetry in addition to translational invariance, see for example, the Nikolaevskii equation [42]. Our system does not present an obvious additional continuous symmetry and further analysis is necessary to decipher the origin of the direct transition to spatiotemporal chaos.
Note that for ξ=0, ℓ is the only length scale present in our system. It fixes the average distance between neighboring protrusions and the thickness of the step profile. In a cell, we thus expect that the average distance between two protrusions is about 500 nm. When reducing the lateral extension L sufficiently, the system forms a single protrusion for ò>ò c , which is now stationary. In that case, new material that is introduced into the system at z=0 is drawn into the protrusion by the contractile lateral activity and cannot form a new protrusion.
In presence of friction between the gel and the membrane, that is, when x ¹ 0, the step profile is also unstable if a z >a c and ò>ò c . The state emerging after the instability can be chaotic as for ξ=0. However, if the friction coefficient exceeds a critical value, then a stationary state emerges at the instability. It consists of a periodic arrangement of protrusions, figure 3 and movie2 (see footnote 5). As for lower values of the friction coefficient, the average distance between protrusions is of the order of h.

Linear stability analysis
To better understand the mechanism underlying the instability of the states homogenous in x, we performed a linear stability analysis of the step profile with ρ 0 =ρ e . In the dimensionless units introduced above, we then have ρ(x, z)=1 for 0z1 and zero otherwise. The stationary velocity profile is in this case =v z 1 z and v x = 0. Note that this state only exists if  a a z c and let a = r r P = x d d 1 x | and a = r r Let us remark that for obtaining the step profile, we necessarily have α z >0 [29]. We consider small perturbations around this profile and write them in the form for 0z1 and 0 otherwise. The growth rate s is given by with real wave numbers q x and q z . Note that the expression for s is independent of the value of the friction coefficient ξ. Its value only affects the range of values of q z that satisfy the boundary condition (5) at z=0. As there is no unique unstable mode but a continuous ensemble (see below), such a discrete selection does not affect the instability threshold. From equation (10) it follows that the step profile is unstable if α x is negative. Since α z >0 and recalling (4), this implies ax >a z as a necessary condition, i.e. an anisotropic active stress that is more contractile in the xdirection. We can give a hand-waving picture of the instability: in this case, a local increase of the gel density leads to a local increase of the contractility in x-direction that is stronger than that in the z-direction. In turn, this generates flows toward the region of increased density, leading to a further increase.
The most unstable modes q=(q x , q z ) are (q x ,0) for any q x and ¥ q , z ( ) for any finite q z . In both cases s=−α/2. In particular, if there is no friction between the gel and the membrane, that is, if ξ=0, then at q z = 0 modes with any wavenumber q x are unstable with the same growth exponent. This is contrary to the usual case with one or a few dominant unstable wavelengths. If the dominant wavelength is finite, then a pattern with a corresponding characteristic length scale emerges, whereas coarsening into a pattern determined by the system size is typically observed if the dominant wavenumber is zero. Since in our system modes at all scales become simultaneously instable and with grow initially with the same rate, a chaotic state can appear directly at the onset of instability of a stationary state. This is the scenario suggested by our numeric solution for the parameter values corresponding to figure 2. The case of finite friction is essentially the same as fro vanishing friction, because all modes (q x ,0) are still allowed and are unstable with the same growth rate as in the case ξ=0.
In figure 4, we present the stability diagram in the plane of the activity parameter a z and the anisotropy ò. . Furthermore,  c saturates for  ¥ a z . For P x z , given by equation (4), we find Saturation is independent of the precise from of P x z , , providing that it can be expanded in powers of ρ: whatever the exponents in (4), dimensional analysis leads to ò=O(1).
The exponential profiles are unconditionally stable. The crucial difference compared to the step profiles is that, here, the steady-state velocity v z (z) reaches a non-zero velocity ¥ v when  ¥ z . Hence, any perturbation dr will be transported at a finite velocity away from the substrate at z=0. As the gel degrades with time at rate k d , the perturbation will reach in finite time a region where the density ρ is exponentially small. The perturbation then satisfies the advection-degradation equation dr dr dr , which implies that it eventually disappears.
Since there is not a single dominant mode that grows fastest, the linear stability analysis cannot yield information about the state after the instability. Even to study whether it is stationary or not requires a nonlinear analysis. This is challenging as the standard techniques of weakly nonlinear analysis cannot be applied in our case. Instead, we can give a hand-waving reasoning for distinguishing between stationary and chaotic states. Let us consider a stationary periodic pattern r x z , per ( )to which we add a perturbation δρ. By expanding the massbalance equation close to the surface, we see that the fate of δρ is determined by the competition of two transport processes. On one hand, potentially unstable terms arise due to v x per and dv x , describing the attraction of neighboring protrusions thanks to contractile effects. On the other hand, stabilizing terms emerge due to the transport of material towards increasing z at velocity v p . Even though v x per and dv x cannot be explicitly determined, the boundary condition implies that they both are of order h x º v v p approach 1 ℓ ( ) . The transition between the chaotic and the periodic pattern occurs when v p and v approach are of the same order, that is, for h x~ℓ. Indeed, For a z <a c with a c denoting the critical activity for the transition from an exponential decaying density to one jumping to zero at z=h [29], the steady state is linearly stable. For a z >a c the steady state is linearly stable as long as the anisotropy ò is below a critical value. Parameters are b=4.9 and ρ 0 =ρ e for a z >a c . For a z <a c , the system evolves into an exponential profile independently of the value of ρ 0 .

Discussion
In this work, we have shown that the dynamics of an active gel polymerizing at a surface is fundamentally different from the dynamics of this system, when fluxes perpendicular to the surface are neglected. In particular, the polymerization velocity v p is not a helpful parameter in a two-dimensional description, such that it cannot give information about the spatial structure of patterns. Hence, results obtained within the thin-layer approximation have to be taken with caution. Let us comment on our assumption of fixed anisotropy of the active gel. In the cytoskeleton, the differences in the active stress in the direction lateral and perpendicular to the surface result from the alignment of the filaments along the membrane. As mentioned already in the part leading to equations (2)-(4), actin filaments will preferentially be oriented with their plus ends towards the membrane, because filament nucleation and elongation is restricted essentially to the immediate vicinity of the membrane. Furthermore, contraction of the actin cortex in the direction perpendicular to the membrane will likely align the filaments in the direction parallel to the membrane. Indeed, filament alignment induced by contraction is a natural feature of the dynamics of active polar or nematic gels [24].
In a full description of the cytoskeleton, the filament alignment should thus be a dynamic field and the anisotropy should emerge spontaneously. The polarization of the actin cortex perpendicular to the membrane and the nematic alignment of the filaments parallel to it can be affected by lateral contractions of the network. However, in the presence of a nearby surface as assumed in our work will probably act as an additional 'force' and lead to a preferential alignment of the filaments with the boundary. In the setting of our work and as long as the state of the system is close to the steady state obtained in [29], we do not expect the dynamics to change qualitatively in presence of a dynamic orientational order field, which could lead to the emergence of nondiagonal terms in the active and the hydrostatic stress. A detailed analysis of the dynamics of the orientational field in particular in the presence of a deformable membrane is postponed to future work.
Compared to previous analyses of thin layers of active gels, we found a direct transition to a chaotic state and a transition controlled by friction. Periodic patterns also emerge in multi-component active systems [31,35,39,43,44]. However, there the period results from balancing friction and diffusion, which is different from our system, where it is a consequence of filament turnover. Spatiotemporal chaos in absence of inertial effects has been reported in active systems [6,[45][46][47][48]. However, these system were incompressible, did not exhibit turnover, and displayed dynamic nematic or polar order. In the two-component system studied in [48] or in presence of nonlinear gel assembly terms [40], a transition between a static periodic and a chaotic state was obtained by changing the friction coefficient. However, in these cases the friction permeated the whole system, whereas in our case, it is restricted to a boundary. Together, these features show that adding a dimension can have effects similar to the introduction of additional physical properties. Note that we restricted our analysis to situations, where the solutions are invariant with respect to translations along the y-direction. In a full threedimensional analysis, we still expect the emergence of chaotic and stationary structures consisting of finger-like protrusions. However, in particular, when considering orientational order as a dynamic field, then the patterns might contain different structures, for example, ridges.
In addition to the activity and the friction, we can also use ρ 0 as a control parameter. In a living cell, this parameter is set by the density of active actin-filament nucleators. In this way the cell has a readily accessible controller to switch between different morphologies. Although a direct observation of this mechanism inside living cells may be tricky, some known phenomena could be manifestations of this instability. First, a 'ramified actin network' consisting of actin bundles emerging from the cortex has been observed in spreading Jurkat T cells [49]. Second, on the surface of suspended Jurkat T cells, short finger-like protrusions emerge continuously [49]. Also, actin cables emanating from the cortex have been observed in frog epithelial cells [50].
Finally, many cells present 'blebs' at their surface, which are membrane bulges in regions with a weakened cortex [51]. In our system, the cortical density is reduced in the regions between the protrusions. Compared to the steady state the thickness of the cortex is roughly halved between protrusions. This has two effects: first it reduces the coupling of the plasma membrane to the actin cortex, second, the actin cortex is predisposed to rupture more easily in regions of a thinned cortex. Both phenomena favor a detachment of the plasma membrane from the cortex and thus the formation of a bleb by hydrostatic pressure [51]. In such a situation one would expect coexistence between protrusions and blebs. Note, that the time needed for formation of a bleb is about 30 s and its lifetime roughly 120 s, which is compatible with the life time of the thinned cortical regions in our analysis. A quantitative analysis is required, though, to check whether the appearance of protrusions and the accompanying thinning of the cortex between these protrusions really offers a mechanism for the appearance of blebs. In upcoming work, we will analyze the dynamics of the actin cortex in case the substrate is deformable to investigate the connection between the cellular structures just mentioned and the cortical instabilities reported in this work.