Spontaneous formation of chaotic protrusions in a polymerizing active gel layer

TL;DR

This paper investigates the dynamics of the actin cortex in cells, revealing that considering vertical fluxes leads to new stability behaviors and the formation of chaotic or periodic protrusions.

Contribution

It introduces a three-dimensional model of the actin cortex, showing how vertical fluxes influence stability and pattern formation, which was not addressed in previous two-dimensional models.

Findings

Isotropic cortex remains stable under large active stresses when vertical fluxes are included.

Lateral contractility exceeding vertical contractility causes protrusions with chaotic or periodic patterns.

Abstract

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.

Equations

$${\partial }_t\rho +{\partial }_x(\rho v_x)+{\partial }_z(\rho v_z)$$

$$\eta [2{\partial }_{xx}{v}_x+{\partial }_{xz}{v}_z+{\partial }_{zz}{v}_x]$$

$$\eta [2{\partial }_{zz}{v}_z+{\partial }_{xz}{v}_x+{\partial }_{xx}{v}_z]$$

$${\Pi }_{x,z}(\rho )=-a_{x,z}{\rho }^3+b{\rho }^4$$

$${\eta {\partial }_z{v}_x|}_{z=0}$$

$$s=-\frac{{\rho }_e}{2}\left(\beta +(\alpha -\beta )\frac{q_x^2}{q_x^2+q_z^2}\right)$$

Full text

The cytoskeleton is a dynamic meshwork of filamentous protein assemblies that is of crucial importance for living cells Bray (2001); Howard (2001); Alberts et al. (2008). 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 Kruse et al. (2004). Physical studies of active matter often use a hydrodynamic approach Marchetti et al. (2013). 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 Kruse et al. (2004); Sanchez et al. (2012); Giomi et al. (2013); Thampi et al. (2013); Pismen (2013); Keber et al. (2014); Pearce et al. (2018), spontaneous flow transitions Voituriez et al. (2007); Fürthauer et al. (2012); Duclos et al. (2018), or buckling instabilities of contracting gels Hannezo et al. (2014); Ideses et al. (2018).

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 animal cells Fritzsche et al. (2016); Chugh et al. (2017) and plays an important role in cellular morphogenesis Salbreux et al. (2012). In addition, the cortex hosts a number of cytoskeletal structures Blanchoin et al. (2014). 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 Zumdieck et al. (2005); Salbreux et al. (2009) and its dynamics when connected to a membrane have been studied Sedzinski et al. (2011); Turlier et al. (2014).

The cortex is formed by actin filaments that polymerize directly at the cell membrane Alberts et al. (2008). Its thickness, which has been reported to be of a few hundred nanometers Clark et al. (2013); Clausen et al. (2017), is determined by the actin polymerization velocity and its disassembly rate Joanny et al. (2013). Due to the large aspect ratio, in theoretical analysis, the cortex is typically considered as a two-dimensional surface Mayer et al. (2010); Naganathan et al. (2014); Berthoumieux et al. (2014); Mietke et al. (2019). However, a formal reduction of the dynamics to the two lateral dimensions as in a thin-film approximation is not always possible.

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 Bois et al. (2011); Mietke et al. (2019). 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 non-linear gel assembly terms Hannezo et al. (2015).

We consider a three-dimensional viscous active gel, growing into the half space $z\geq 0$ at the surface $(x,y,z=0)$, Fig. 1. For the sake of simplicity, we assume invariance along the direction $y$. The state of the gel is determined by the density $\rho$ and the velocity field $\mathbf{v}=(v_{x},v_{z})$. Mass conservation and force balance yield:

[TABLE]

In these equations, $k_{d}$ is the gel’s disassembly rate, $\eta$ denotes the viscosity and $\Pi_{x,z}$ are the components of the non-viscous contribution to the total stress, which we assume to have only diagonal components. It has two contributions: an effective hydrostatic pressure and the stress generated by active processes in the gel. In general, both contributions could be anisotropic due to local filament alignment. For simplicity, we consider an isotropic hydrostatic pressure and an active stress with fixed anisotropy. We discuss the limitations of this assumption at the end of our manuscript. Similar to Ref. Joanny et al. (2013), we choose

[TABLE]

with $b>0$. We take $a_{x,z}>0$ reflecting the contractile nature of the active stress. Consequently, $\Pi_{x}$ and $\Pi_{z}$ decrease with $\rho$ for low enough densities and increase with $\rho$ as $\rho\to\infty$. We checked that our results do not depend on the specific form chosen for $\Pi_{x,z}$. It suffices that the dependence of $\Pi_{x,z}$ on $\rho$ has a shape similar to that given in Eq. (4).

The dynamic equations are completed with the boundary condition $\rho(z=0)=\rho_{0}$ and $v_{z}(z=0)=v_{p}$, where $\rho_{0}$ is the density of nucleators at the boundary and $v_{p}$ the polymerization speed. In addition, the tangential stress at the surface is balanced by friction of the gel along the membrane:

[TABLE]

where $\xi$ is the friction constant. Furthermore, we impose periodic boundary conditions in the $x$-direction with period $L$. In the $z$-direction, the system extends to infinity and $\rho\to 0$ for $z\to\infty$.

If we ignore the $z$-direction in Eqs. (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 Bois et al. (2011); Mietke et al. (2019). For high enough activity, these systems typically generate contracted stationary states with a single region of high gel density unless non-linear assembly terms are considered Hannezo et al. (2015). However, there is no way to integrate the dynamic equations (1)-(3) over $z$ to get a closed system for the averaged fields $\bar{\rho}(x)$ and $\bar{v}(x)$ in one dimension: a formal integration over $z$ creates term such as $\int v_{x}\rho\,\mathrm{d}z$, $\int\Pi(\rho)\,\mathrm{d}z$ and boundary terms such as $\partial_{z}v_{x}|_{z=0}$ or $\partial_{x}v_{x}|_{z=0}$. These terms cannot be readily expressed in terms of $\bar{\rho}(x)$ and $\bar{v}(x)$ for a compressible system.

In the following we will use a dimensionless form of the dynamic equations (1)-(4) and scale time by $k_{d}^{-1}$, length by $\ell=v_{p}k_{d}^{-1}$, concentration by $\rho_{0}$, and stress by $\eta k_{d}$. Hence, the relevant parameters are $a_{x,z}(k_{d}\eta)^{-1}\rho_{0}^{3}$, $b(k_{d}\eta)^{-1}\rho_{0}^{4}$, and $\xi v_{p}(k_{d}\eta)^{-1}$. We will use the same notation for the original and the rescaled parameters. For future use, we introduce the anisotropy parameter $\epsilon=a_{x}/a_{z}$.

In case the gel is invariant with respect to lateral translations, two qualitatively different stationary states exist: either the density decays exponentially or it jumps from a finite value $\rho_{e}$ to zero at $z=h\sim\ell$ Joanny et al. (2013). For the ”exponential profile”, the velocity $v_{z}$ converges to a finite value for $z\to\infty$. 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_{c}\geq 0$, and when the density of nucleators $\rho_{0}$ exceeds a critical value $\rho_{c}$. Otherwise, the profile decays exponentially.

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 $\partial_{z}(\rho v_{z})$ is obtained by an up-wind finite-differences scheme in real space. 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 $\Delta x=0.004$, $\Delta z=0.007$, and $\Delta t=0.0005$. We checked that our results do not change for smaller values. As initial condition, we take $\rho(x_{i},z_{j})=\rho_{0}(z_{j})(1+\tilde{\eta}(x_{i},z_{j}))$, where $\tilde{\eta}(x_{i},z_{j})$ are independent and uniformly distributed random numbers between $-0.05$ and $0.05$.

We first consider the limit of zero friction, i.e., $\partial_{z}v_{x}(z=0)=0$. If $a_{z}>a_{c}$, then there is a critical anisotropy $\epsilon_{c}$, such that the step profile is unstable for $\epsilon>\epsilon_{c}$. In Figure 2 and Movie 1 SM , we present the solution for $a_{z}=7>a_{c}$ and $\epsilon=1.18>\epsilon_{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, Fig. 2, the system does not settle in a periodic state. Instead the dynamics appears to be chaotic. The chaotic nature of our system is highlighted by its sensitivity to initial conditions, Fig. 2c, and the relaxation of the protrusion auto-correlation $C(t)$, Fig. 2d, with $C(t)=\langle(\sigma(x,\tau)-\langle\sigma\rangle)(\sigma(x,\tau+t)-\langle\sigma\rangle)\rangle$, where $\sigma(x,t)$ is 1 inside a protrusion and 0 outside. Furthermore, $\langle\ldots\rangle$ denotes a space and time average. Our numerical results indicate that the chaotic state presented in Fig. 2 emerges directly at the critical anisotropy $\epsilon_{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 Aranson and Kramer (2002); Tribelsky and Tsuboi (1996).

Note that for $\xi=0$, $\ell$ is the only length scale present in our system. It fixes the order of the density of protrusions and the thickness of the step profile. When reducing the lateral extension $L$ sufficiently, the system forms a single protrusion for $\epsilon>\epsilon_{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. Also in the case of large enough friction $\xi$, a stationary state emerges at the instability. It consists of a periodic arrangement of protrusions, Fig. 3 and Movie 2 SM .

We perform a linear stability analysis of the step profile with $\rho_{0}=\rho_{e}$. Then $\rho(x,z)=\rho_{0}$ for $0\leq z\leq h$ and zero otherwise. Note that this state only exists if $a_{z}\geq a_{c}$. Let $\alpha=\frac{d\Pi_{x}}{d\rho}|_{\rho_{e}}$ and $\beta=\frac{d\Pi_{z}}{d\rho}|_{\rho_{e}}$. Let us recall that for obtaining the step profile, we necessarily have $\beta>0$ Joanny et al. (2013). Then, the growth rate $s$ of a perturbation $\tilde{\rho}\exp(iq_{x}x+iq_{z}z)$ for $0\leq z\leq h$ is given by

[TABLE]

with real wave numbers $q_{x}$ and $q_{z}$. Note that the expression for $s$ is independent of the value of the friction coefficient $\xi$. Its value only affects the range of values of $q_{z}$ that satisfy the boundary condition at $z=0$.

From Equation (6), the step profile is unstable if $\alpha$ is negative. Together with $\beta>0$ this implies $a_{x}>a_{z}$, i.e., an anisotropic active stress that is more contractile in the $x$-direction. In this case, a local increase of the gel density leads to a local increase of the contractility in $x$-direction. In turn, this generates flows toward the region of increased density, leading to a further increase. The most unstable modes $\mathbf{q}=(q_{x},q_{z})$ are $(q_{x},0)$ for any $q_{x}$ and $(\infty,q_{z})$ for any finite $q_{z}$. Contrary to the usual case with one or a few dominant unstable wavelengths, where a pattern with a corresponding characteristic length scale emerges, here an infinite number of modes becomes simultaneously unstable. In this case, the steady state can lose stability directly in favor of a chaotic state as our numeric solution suggests for the parameter values corresponding to Fig. 2.

In Figure 4, we present the stability diagram in the plane of the activity parameter $a_{z}$ and the anisotropy $\epsilon$. Even though the precise form of the stability boundary $\epsilon_{c}(a_{z})$ depends on the detailed expression of the functions $\Pi_{x,z}$, we always have $\epsilon_{c}-1\propto\sqrt{|a_{z}-a_{c}|}$ for $a_{z}\gtrsim a_{c}$. Indeed, $\frac{d\Pi_{z}}{d\rho}|_{\rho_{e}^{c}}=0$, where $\rho_{e}^{c}\equiv\rho_{e}$ for $a_{z}=a_{c}$, so that for $a_{z}=a_{c}+\delta a$ we get $\rho_{e}-\rho_{e}^{c}\propto\sqrt{\delta a}$. Furthermore, $\epsilon_{c}$ saturates for $a_{z}\to\infty$ independently of the precise from of $\Pi_{x,z}$. For $\Pi_{x,z}$ given by Eq. (4), we find $\epsilon=4b\rho_{e}/a_{z}=O(1)$ for $a_{z}\to\infty$.

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_{\infty}$ when $z\to\infty$. Hence, any perturbation $\delta\rho$ 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 a region where the density $\rho$ is exponentially small in finite time. The perturbation then satisfies the advection-degradation equation $\partial_{t}\delta\rho+v_{\infty}\partial_{z}\delta\rho=-k_{d}\delta\rho$, 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 non-linear analysis. This is challenging as the standard techniques of weakly non-linear 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 $\rho^{\mathrm{per}}(x,z)$ to which we add a perturbation $\delta\rho$. By expanding the mass-balance equation close to the surface, we see that the fate of $\delta\rho$ is determined by the competition of two transport processes. On one hand, potentially unstable terms arise due to $v_{x}^{\mathrm{per}}$ and $\delta v_{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}^{\mathrm{per}}$ and $\delta v_{x}$ cannot be explicitly determined, the boundary condition implies that they both are of order $v_{\mathrm{approach}}\equiv v_{p}\eta(\ell\xi)^{-1}$. The transition between the chaotic and the periodic pattern occurs when $v_{p}$ and $v_{\mathrm{approach}}$ are of the same order, that is, for $\eta/\xi\sim\ell$. Indeed, $0.2=v_{\mathrm{approach}}<v_{p}$ in Fig. 3, whereas $\infty=v_{\mathrm{approach}}>v_{p}$ in Fig. 2.

In this letter, 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. 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. In a full description of the cytoskeleton, the filament alignment is a dynamic field and the anisotropy should emerge spontaneously. Filament alignment induced by contraction is a natural feature of the dynamics of active polar or nematic gels Salbreux et al. (2009). In addition, the presence of a surface should also lead to a preferential alignment of the filaments at the boundary. We do thus not expect the dynamics to change qualitatively in presence of a dynamic orientational order field and postpone its analysis 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 Kruse and Jülicher (2003); Bois et al. (2011); Gowrishankar et al. (2012); Kumar et al. (2014); Banerjee et al. (2017). 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 Dombrowski et al. (2004); Sanchez et al. (2012); Neef and Kruse (2014); Ramaswamy and Jülicher (2016); Doostmohammadi et al. (2016). However, these system were incompressible, did not exhibit turnover, and displayed dynamic nematic or polar order. In the two-component system studied in Ref. Doostmohammadi et al. (2016) or in presence of non-linear gel assembly terms Hannezo et al. (2015), 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.

In addition to the activity and the friction, we can also use $\rho_{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 Fritzsche et al. (2017). Second, on the surface of suspended Jurkat T cells, short finger-like protrusions emerge continuously Fritzsche et al. (2017). Finally, many cells present “blebs” at their surface, which are membrane bulges in regions with a weakened cortex Charras (2008). In our system, the cortical density is reduced in the regions between the protrusions, which offers a possible mechanism for the appearance of blebs. In upcoming work, we will analyse 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.