Introduction

Actin networks are remarkably dynamic and versatile and organize in a variety of architectures to accomplish crucial cellular functions [1]. For instance, isotropic thin actin gels form the cell cortex, which largely determines cell shape [2] and motility in confined non-adherent environments [3, 4]. Polar structures at the edge of adherent cells, either forming filaments as in filopodia or sheets as in lamellipodia [5], enable cells to probe their environment and crawl on substrates. Nematic actin bundles conform a variety of contractile structures [6], including sub-cellular rings during cytokinesis [7] and cortical repair [8], supra-cellular rings during wound healing [9] or development [10], bundle networks during cellularization [11], or stress fibers in adherent cells [1214]. Nematic bundles consist of highly aligned and densely packed actin filaments of mixed polarity connected by a diversity of crosslinkers. In vivo and in vitro observations show the key role of actin nucleators and regulators, of myosin activity and of crosslinkers in the assembly and maintenance of actin bundles [1, 2, 5, 6, 13, 1521].

Various studies have emphasized the morphological, dynamical, molecular and functional specificities of different types of actin bundles such as dorsal, transverse, and ventral stress fibers or contractile rings [13, 14, 2224]. Here, we ask the question of whether, despite this diversity, the ubiquity of actin bundles in different contexts can be explained by the intrinsic ability of the active actomyosin gel to self-organize into patterns of dense nematic structures. Suggestive of such active self-organization, stress fibers often form dynamic highly organized patterns, e.g. involving families of fibers along orthogonal directions [12, 14, 19, 22, 2527]. Other kinds of actin bundles also develop patterns of remarkable regularity. For instance, parallel arrangements of 2 μm-spaced actin bundles serve as templates for extracellular matrix deposition during butterfly wing morphogenesis and determine their iridescence [28]. Similarly, the morphogenesis of the striated tracheal cuticle in Drosophila is pre-patterned by a parallel arrangement of actin bundles spaced by ~1μm, spanning from sub-cellular to supra-cellular and organ scales [21, 29]. Muscle-like actin bundles form regular parallel patterns spanning organs as in C. elegans [19] or the entire organism of hydra [30]. Furthermore, dense nematic bundles have been shown to assemble de novo from the sparse isotropic cortex in a process controlled by myosin activity [20], and to form a mechanically integrated network with the cortex [31].

Theoretical models of stress fibers often assume from the outset the existence of dense actin bundles, possibly embedded in an isotropic network [32], although previous work has studied the emergence of nematic rings as a result of localized activity gradients [33]. In the field of active nematics, many theoretical studies have examined the well-known hydrodynamical bend/splay instability arising in uniform active systems with pre-existing long-range orientational order [3436]. However, the question of how patterns of orientational order arise from an isotropic and uniform system as a result of activity has received very little attention, and has been based so far on models that do not capture fundamental physical features of actomyosin gels as discussed below [37, 38]. Here, to understand the mechanisms underlying the self-assembly of dense actin bundles from a low-density isotropic cortex, we develop a theory for the self-organization of dense nematic structures in the actomyosin cytoskeleton based on a nematic active gel theory accounting for density variations and compressibility. In our theory, nematic patterning is driven by activity rather than by a more conventional crowding mechanism. Linear stability analysis and fully nonlinear simulations show that, when coupled to nematodynamics, the well-known patterning mechanism based on self-reinforcing flows [29, 3943] leads to a rich diversity of patterns combining density and nematic order. The geometry and dynamics of the emergent patterns are very similar to those observed in diverse cellular contexts. Finally, we test key assumptions of our phenomenological theory leading to such self-organization using discrete network simulations.

Theoretical Model

The actomyosin cytoskeleton can be understood as a compressible, active and viscous fluid gel with orientational order undergoing turnover [44, 45]. Symmetry-breaking and pattern formation in actomyosin gels is often mediated by the emergence of an advective instability leading to compressible self-reinforcing flows. According to this mechanism, fluctuations in the density of cytoskeletal active units generate gradients in active stress, driving flows, which in turn advect the active units reinforcing the initial fluctuation [29, 3941]. This kind of advective instability has been invoked to explain cell polarization and amoeboid motility [4, 41, 46], the formation of the cytokinetic ring [42], the formation of periodic dense actin structures during tracheal morphogenesis in Drosophila [29], or the emergence of self-sustained dynamical states in actomyosin gels extracted from cells [47, 48], and has been recently reproduced to some degree in confined reconstituted gels from purified proteins [49]. In all of these examples, the actomyosin gel develops sustained compressible flows converging towards regions of high density. Furthermore, the compressive strain rate induced by these active flows has been shown to drive nematic order [7, 33]. Finally, observations on adherent cells show that active contractility is required for actin bundle formation [19, 20, 22]. Therefore, a theoretical model to understand the self-organization of dense nematic structures in actomyosin gels should consider a compressible and density-dependent fluid capturing the advective instability mentioned above. Furthermore, this model should acknowledge the active nature of the assembly of dense nematic structures, and permit extended isotropic phases commonly observed in the actin cortex, possibly coexisting with dense nematic phases [20, 31], rather than thermodynamically enforce high nematic order everywhere except at topological defects [36, 5052].

Previous models for dry and dilute aligning active matter develop density patterns [37, 53, 54] but fail to capture the hydrodynamic interactions of actomyosin gels, whereas models for active nematic fluids either ignore density [38, 5561], or account for the density of active particles suspended in an incompressible flow [34, 35, 62, 63], and therefore cannot describe the advective instability and self-reinforcing flows of actomyosin gels. Previous models describing the emergence of nematic patterns from uniform and isotropic states either ignore hydrodynamics [37], or the density and flow compressibility characteristic of actomyosin gels [38], and therefore result in very different instability mechanisms to those presented here. To model a thin layer of actomyosin gel, we summarize next a minimal theory for 2D density-dependent compressible active nematic fluids. In [64], we provide a systematic derivation of this theory based on a variational formalism of irreversible thermodynamics. This model can be understood as a density-dependent and compressible version of the active nematic theory presented in [57], or as a nematic generalization of the isotropic theory used in [29]. As elaborated in [64], it is possible to develop an alternative compressible and density-dependent active nematic theory based on the Beris-Edwards formalism [38, 51, 55, 56, 58, 59, 63].

In our model, the local state of the system is described by the areal density of cytoskeletal material ρ(t, x) and by the network architecture given by the symmetric and traceless nematic tensor q(t, x), see Fig. 1a, both of which depend on time t and position x. A more detailed model could consider separate densities for actin, myosin, and possibly other structural or regulatory proteins. Likewise, in principle, the orientational order of actin and myosin filaments could be described by separate nematic tensors. The nematic tensor can be expressed as qij = S (ninjδij/2), where n is the average molecular alignment, the degree of local alignment about n and δij is the identity. We denote by v(t, x) the velocity field of the gel. The rate-of-deformation tensor measures the local rate of distortion of the fluid, whereas measures its local spin, where is the nabla differential operator. The deviatoric part of the rate-of-deformation tensor is defined by . The rate of change of q relative to a frame that translates and locally rotates with the flow generated by υ is given by the Jaumann derivative .

Because we consider a bi-periodic domain Ω, we ignore other boundary conditions. Balance of cytoskeletal mass for a compressible fluid undergoing turnover takes the conventional form [29]

where the second term models advection of cytoskeletal material by flow, the third term models diffusion with D being an effective diffusivity, and the last term models cytoskeletal turnover, where ρ0 is the steady-state areal density and kd is the depolymerization rate. Note that for a uniform, steady-state and quiescent gel, the first three terms vanish and ρ(t, x) = ρ0.

Key model ingredients

(a) The local state is defined by areal density ρ and by orientational order quantified by the nematic parameter S and by the nematic direction n. (b) Isotropic active tension λ when the network is isotropic (S = 0), and (c) anisotropic tension when S ≠ 0, controlled by κ. Active tension is positive (contractile) in all directions whenever |κ| < 1, but its deviatoric part is contractile when κ > 0 and extensile when κ < 0. Orientational order is driven by (d) active forces conjugate to nematic order and characterized by parameter λ⊙ and by (e) passive flow-induced alignment in the presence of deviatoric rate-of-deformation with coupling parameter β.

Force balance in the gel can be expressed as

where γ > 0 models a viscous drag with the surroundings (e.g. the plasma membrane) and σ = σanem + σdiss + σact is the stress tensor in the gel, which in 2D has units of tension and which includes a contribution coming from the nematic free-energy, a dissipative contribution, and an active contribution.

The nematic stress σnem follows from a standard derivation adapted here to a density-dependent material. It derives from the free-energy ö = Ωρf(q, ∇q)dS, where

is the classical Landau expansion of free-energy density per unit mass, with a and b > 0 susceptibility parameters, and L > 0 the Frank constant. When a < 0, the first term favors equilibrium nematic ordering, e.g. due to crowding in very dense gels of elongated filaments. Otherwise, the susceptibility terms entropically favor isotropic states with small S. In the actin cytoskeleton, this can model the random orientation of filaments as a result of the entropic fluctuations of filaments and their nucleators. The last term penalizes sharp gradients in the nematic field, which can result from the bending rigidity of actin filaments. A lengthy but standard calculation leads to the so-called molecular field [64]

and to the explicit form of the nematic non-symmetric stress tensor

The dissipative part of the stress

includes a viscous stress controlled by the gel viscosity parameter η [33] and a stress resulting from changes in the nematic field controlled by the coupling parameter η < 0 [7]. The term involving η can be understood as a stress in the gel arising from the drag induced by filaments as they reorient relative to the underlying hydrodynamic flow. Finally, we assume that the active stress resulting from the mechanical transduction of chemical power has an isotropic component and an anisotropic component oriented along the nematic tensor following

Activity parameter λ controls the isotropic tension and is contractile for λ > 0, as assumed here. The additional activity parameter λaniso, or equivalently κ = λaniso/λ, controls the deviatoric (traceless) component of active tension. This component is anisotropic and can be positive or negative parallel or perpendicular to the nematic direction depending on the sign of λ. When |λ| < 1, then the total active tension remains positive in all directions, with a larger magnitude parallel to the nematic direction when λ > 0 (contractile deviatoric component) and perpendicular to it when λ < 0 (extensile deviatoric component). We note that the isotropic component of active tension is meaningful here because our active gel is compressible. When order is low (S ≈ 0), active tension is isotropic, Fig. 1b, whereas when order is high, active tension becomes anisotropic, Fig. 1c. We can interpret that active tension along the nematic direction reflects the sliding of antiparallel fibers driven by myosin motors, and whereas active tension perpendicular to it reflects the the out-of-equilibrium binding of bundling proteins or myosins [5, 6571].

Balance of the generalized forces power-conjugate to also includes viscous, elastic-nematic and active contributions, and takes the form

In this expression, ηrot is a nematic viscous coefficient. The second term models alignment induced by the rate of deformation of the flow, e.g. with compression/extension driving alignment perpendicular/parallel to the velocity gradient (Fig. 1e), as experimentally observed in [7]. This term involves the same coefficient β as the last term in Eq. (6) because of Onsager’s reciprocity relations, and the entropy production inequality requires that 2η ηrotβ2 ≥ 0 [64]. The third term is a thermodynamic force driven by F. In agreement with the observations that nematic ordering in actomyosin gels is actively driven, we assume a > 0, and therefore Eqs. (4,8) show that this term tends to restore isotropy. The last term is an active generalized force controlled by activity parameter λ ≥ 0 tending to further align filaments (Fig. 1d) [7]. This term is linear in ρ because in the expansion the constant contribution can be subsumed by the susceptibility parameter a [33]. Thus, the active term acts as a negative density-dependent susceptibility. When c0 = 2aρ0λ < 0, the system can sustain a uniform quiescent state with ρ(x, t) = ρ0, υ(x, t) = 0 and a non-zero nematic order parameter . Even if c0 > 0, and hence the uniform quiescent state is devoid of order, pattern formation can induce density variations such that 2aρλ becomes negative locally and actively favors local nematic order. Physically, the term −ρλq implements the notion that the binding of a bundling protein, which drives active alignment, is more probable when two filaments are in close proximity and nearly aligned, a situation favored by high density and nematic order.

The nonlinearity of the coupled system of partial differential equations given by the balance laws Eqs. (1,7,8), along with the constitutive relations in Eqs. (4,5,6,7), has different sources summarized next. The theory presents nonlinearities intrinsic to transport equations in the advective term of Eq. (1) and in the definition of the Jaumann derivative of the nematic tensor. Furthermore, nonlinearities in q in the constitutive relations results from the standard nematic free-energy adopted here. Our hypothesis that material properties in the gel are proportional to density and our thermodynamically-consistent derivation of the theory [64] result in further nonlinearities involving density in the constitutive relations. Finally, the nonlinearity involving density and the nematic field in the last term of Eq. (8) has been discussed in the previous paragraph.

Our theory has three active parameters, λ, κ and λ, all reflecting the conversion of chemical power into mechanical power in the network. The magnitude and the modes of chemomechanical transduction should depend on the molecular architecture of the network [72, 73], e.g. the stoichiometry of filaments, cross-linkers, and myosins, or the length distribution of filaments. Accordingly, we allow these parameters to vary independently.

By freezing an isotropic state, S = 0, our model reduces to an orientation-independent active gel model, which develops periodic patterns driven by self-reinforcing active flows sustained by diffusion and turnover [29]. In the present model, however, translational, orientational and density dynamics are intimately coupled through the terms involving β, λ, κ and L.

We readily identify the hydrodynamic length , above which friction dominates over viscosity, the Damkolher length above which reactions dominate over diffusion, and the nematic length lq = . Non-dimensional analysis reveals a set of non-dimensional groups that control the system behavior, namely the non-dimensional turnover rate , the Frank constant , the susceptibility parameters and , the drag coefficients and , the nematic activity coefficient , and the active tension parameters and κ (Appendix B). The full list of material parameters for each figure given in Tables I and II and justified in Appendix D.

Results

Onset and nature of pattern formation

To examine the role of nematic order in the emergence of various actin architectures, we performed linear stability analysis of our model particularized to 1D, whose dynamical variables are velocity, density and nematic order, v(x,t), ρ(x, t), and q(x,t), along x, where q > 0 (< 0) corresponds to a nematic orientation n parallel (perpendicular) to the x-axis (Appendix A). We first focused on the case c0 = 2aρ0λ > 0 to examine the loss of stability of a uniform, isotropic, and quiescent steady state (ρ(x, t) = ρ0, v(x,t) = 0, q(x,t) = 0) by increasing the master activity parameter λ and identifying the most unstable modes. This allowed us to determine a threshold activity for pattern formation and the wavelength of the emerging pattern (Appendix C). Since the exact evaluation of such quantities requires solving nonlinear equations, we derived explicit expansions in the limit of small L for the critical contractile activity

where and δ = γDκβ/(2ηc0), and for the the corresponding wavenumber

where .

When κ = 0 or β = 0, and hence δ = 0, we recover the predictions of an active gel model not accounting for network architecture [29], λcrit = λcrit,0 and λcrit = λcrit, 0. The expression for λcrit,0 shows that the instability takes place when activity is large enough to overcome the effect of diffusion and turnover, tending to uniformize density, and that of friction and viscosity, tending to suppress flows. Active tension anisotropy (κ ≠ 0) and flow-induced alignment (β < 0) couple the instability of [29] to nematic order, changing the nature of pattern formation, see Eqs. (9,10) and the expression for δ. This leads to quantitative changes in critical tension and wavenumber, which can be very significant depending on the ratio of hydrodynamic and Damkolher lengths and on the strength and sign of nematic coupling. The nematic corrections increase as c0 → 0, close to the point where the uniform quiescent state develops spontaneous order. We thus studied separately the regime 0 < c0 ≪ 1, finding analogous expansions for the critical tension and wavenumber in terms of δ = γDκβ/(2ηL). Interestingly, Eq. (9) shows that the activity threshold is reduced, and hence pattern formation facilitated, when κ < 0, i.e. when active tension is larger perpendicular to the nematic direction. Besides modifying critical tension and wavenumber, the present model predicts that the dynamical modes with self-reinforcing flows generate patterns where high density co-localizes with high nematic order.

To test the validity of this analysis and further understand the system beyond the onset of pattern formation, we performed fully non-linear finite element simulations in a periodic square 2D domain [64, 74], with a domain size chosen to be 8 pattern where pattern = 2π/vcrit is the pattern lengthscale predicted by linear stability analysis. In these simulations, we increased the activity parameter A beyond the instability starting from a quiescent uniform state. We found that the linear stability analysis very accurately predicts the activity thresholds within two percent across a wide range of parameters. Furthermore, we found that the linear estimate for the wavenumber in Eq. (10) characterizes well the nonlinear patterns, as quantified in Supp. Fig. 1 and illustrated visually by the density patterns in Figs. 2 and 4 showing the entire 8pattern × 8pattern computational domain.. In the nonlinear regime, the exponentially-growing instabilities eventually reach out-of-equilibrium quasi-steady-state patterns involving self-reinforcing flows towards regularly-spaced regions of high density surrounded by a low density matrix, as previously reported for various compressible and density-dependent active gel models [29, 3943]. In the absence of nematic coupling (β = κ = λ = 0), these high density domains are droplets arranged in a regular hexagonal lattice with S = 0 throughout the domain as previously reported [40], Movie 1. In contrast, for a generic parameter set with finite β, κ and λ, high density domains adopt elongated configurations or bands where order is high, surrounded by a low density and low order matrix, Fig. 2b and Movie 2. Thus, our nematic active gel develops out-of-equilibrium localized nematic states starting from an isotropic state through a mechanochemical mechanism, which unlike those in [37, 38], exhibit localization of both density and nematic order, and hence resemble dense nematic structures embedded in a low-density isotropic actin cortex. Furthermore, the shape and internal architecture of dense and nematic phases are qualitatively modified by the nematic coupling.

Our simulations show that self-reinforcing flows develop along the direction of largest active tension, and consequently the pattern architecture depends on the sign of κ, Fig. 2c. For κ < 0, the system self-organizes into high-density and high-order bands, where nematic direction is parallel to their axis, in what we call fibrillar pattern, Fig. 2b(I). Instead, for κ > 0 nematic order is perpendicular to the axis of the bands, in what we call banded pattern, Fig. 2b(III). To systematically study the effect of active tension anisotropy, we varied κ between −0.8 and 0.8 while keeping all other non-dimensional groups fixed and setting λ to be 1.3 times the critical activity. We defined the order parameter ω, Fig. 2a, allowing us to distinguish between banded (ω > 0) and fibrillar (ω < 0) organizations. We found a sharp transition between fibrillar and banded regimes around κ ≈ 0, during which elongated high-density and high-order domains fragment into nematic droplets or tactoids [75, 76], Fig. 2b(II).

Active patterns coupling nematic order and density driven by self-reinforcing flows.

(a) Illustration of the dimensionless parameters characterizing active tension anisotropy (κ) and pattern architecture, quantified by the relative orientation of nematic order and high-density structures . (b) Order parameter of pattern architecture ω as a function of active tension anisotropy κ obtained from nonlinear simulations, showing transition from states with nematic direction parallel to high-density structures (ω < 0, fibrillar patterns) for κ < 0 to states with nematic direction perpendicular to high-density structures (ω > 0, banded patterns with perpendicular nematic organization) for κ > 0. Because |κ| < 1, the active tension is always positive in all directions. (c) Map of density, nematic order S, nematic direction (red segments) and flow field (green arrows) for quasi-steady fibrillar (I) and banded (III) patterns, and for a transition pattern of high density droplets with high nematic order (II) corresponding to nearly isotropic active tension. (d) Illustration of the out-of-equilibrium quasi-steady states, maintained by self-reinforcing flows, diffusion and turnover.

Our results for κ < 0, leading to self-organized dense nematic fibrillar patterns from an isotropic low-density network, are in agreement with evidence suggesting that stress fibers can assemble from the actin cortex without the involvement of stress fiber precursors or actin polymerization at focal adhesions [20]. They also agree with observations showing that actin bundles form a mechanical continuum with the surrounding sparse and isotropic cortex [31]. Their morphology and patterning dynamics is strikingly reminiscent of actin microridges, formed at the apical surfaces of mucosal epithelial cells [77, 78]. As discussed earlier, the condition κ < 0 implies that the deviatoric component of active tension is extensile. In agreement with widely studied incompressible extensile nematic systems, here the material is drawn perpendicular to the nematic direction, but rather than being expelled along the nematic direction as required from incompressibility [36], here it is recycled by disassembly and turnover (Fig. 2c,d left and Eq. (1)). Finally, we note the similarity in terms of density and nematic architecture between our fibrillar patterns and those emerging in other active systems through different mechanisms of self-organization, including polar motile filaments [79, 80] or mean-field models of dry mixtures of microtubules and motors [81].

Requirements for fibrillar and banded patterns

At linear order, our theory shows that the distinctly nematic self-organization requires both flow-induced alignment (β) and active tension anisotropy (κ), whereas no condition is required on nematic activity (λ). We performed further simulations to establish the requirements for fibrillar and banded active patterning in the nonlinear regime. We found that both banded and fibrillar patterns readily form for β = 0 and finite κ, yet a finite value of β enhances fibrillar formation, leading to longer and more stable dense bands, and hinders banded organization, Movie 3. This behavior is expected since the velocity gradients of the self-reinforcing flows tend to align filaments parallel to high-density bands due to βddev in Eq. (8).

Pattern formation for a range of values of anisotropic active parameter κ in the limit λ → 0.

(a) Pattern formation for material parameters used in Fig. 2 except for λ = 0 while leaving c0 unchanged. (b) Here, in addition to λ = 0, we set β2 = 2ηηrot to the largest value allowed by the entropy production inequality. (c) Same parameters as in (b), except for an increase in friction as detailed in Table I.

The active nematic coefficient λ modifies the onset of pattern formation through c0 according to linear stability analysis. Apart from this linear effect, it should also contribute to condensation of nematic order in high-density regions since it appears multiplied by p in Eq. (8). To examine this nonlinear effect, we perform simulations with λ = 0 but keeping c0 fixed. This leads to very different patterns without clear elongated structures and very mild nematic patterning, Fig. 3a. Enhancing nematic patterning by considering the largest thermodynamically allowed value of |β| leads to elongated structures for κ < 0, but rather than high-order co-localizing with high density, the nematic field develops domains with 90° angles between high- and low-density regions, Fig. 3b, as observed in other active nematic systems [37]. This architecture is enhanced by higher friction γ, Fig. 3c. Thus, the architectures found for λ = 0 and κ < 0 are distinct from the fibrillar pattern described previously. Similarly, rather than bands with order perpendicular to their axis, for λ = 0, κ > 0 and high |β| we found patterns of nematic asters (high-density droplets with radial nematic organization around them), Fig. 3b,c.

Together, these results show that active tension anisotropy (κ ≠ 0) and nematic activity (λ = 0) are necessary and sufficient for nematic self-organization into fibrillar or banded patterns, with flow-induced alignment (β < 0) favoring fibrillar organization. Furthermore, they show that while Eqs. (9,10) accurately predict the activity threshold and the wavenumber of the nonlinear patterns, the linear stability analysis does not capture the architecture of the patterns, nor the key role of λ , in the nonlinear regime.

Morphological and dynamical diversity of self-organized fibrillar patterns

Given the morphological and dynamical diversity of nematic bundles in actin gels across cell types, geometric confinement, mechanical environment, or biological and pharmacological treatments [11, 26, 27, 8284], we varied model parameters to examine the architectures predicted by our active gel model, focusing on κ < 0. Significant changes in the effective parameters of our active gel model are reasonable since the active mechanical properties of actomyosin gels strongly depend on micro-architecture both in reconstituted systems and in cells [70, 72].

With our default parameter set, bundle junctions and free ends are unfavorable and reorganize during pattern formation to annihilate as much as possible, Movie 4(I). However, because the initial state of the system is isotropic but fibrillar patterns are not, this coarsening process leads to frustrated labyrinth patterns with domains and defects, which depending on parameters can remain frozen in quasi-steady states as in Fig. 4a(I). We then asked the question of whether an orientational bias, which physically may be caused by cytoskeletal flows, boundaries or directed polymerization [22], could direct the pattern and result in fewer defects. We first slightly modified the system by including a small anisotropic strain-rate, according to which friction is computed relative to an elongating background. This directional bias is sufficient to produce well-oriented defect-free patterns aligned with the direction of elongation, Fig. 4a(II) and Movie 4(II). Alternatively, we considered c0 to be slightly negative, leading to the uniform and nematic steady state ρ(x,t) = ρ0, v(x,t) = 0 and . The linear stability analysis around this state and further nonlinear simulations show that the essential phenomenology of Eqs. (9,10) and Fig. 2 is not altered by the slight pre-existing order (Appendix C). Again, the weak pre-existing order provides sufficient bias to direct pattern orientation, Fig. 4a(III). Thus, our model indicates that an anisotropic bias can guide and anneal nematic fibrillar patterns, in agreement with remarkably oriented patterns of actin bundles in elongated cells [28], as a result of uniaxial cellular stretch [19, 78], or on anisotropically curved surfaces [21, 29].

Control of nematic bundle pattern orientation, connectivity and dynamics.

(a) Effect of orientational bias. (I) A uniform isotropic gel self-organizes into a labyrinth pattern with defects. (II) A small background anisotropic strain-rate efficiently orients nematic bundles. (III) A slight initial network alignment (S0 = 0.05) orients bundles, which later loose stability, bend, and generate/anneal defects. See Movie 4. We recall that the nematic order parameter in the quiescent and uniform equilibrium state is S0 = 0 if c0 ≥ 0 and otherwise, with c0 = 2aρ0λ. (b) Depending on active tension anisotropy, nematic bundles are contractile and straighten (I, κ = −0.2), leading to quasi-steady networks, or extensile and wrinkle (II, κ = −0.8), leading to bundle breaking and recombination, and persistently dynamic networks (III). See Movie 5. The contractility or extensibility of the nematic bundles can be appreciated in the maps of the difference between the stress along the nematic direction (σ||) and the stress perpendicular to it (σ⊥), normalized by the constant σcrit = λcritρ0. (c) Promoting mechanical interaction between bundles. (I) Dynamical pattern obtained by reducing friction, and thereby increasing and . Time sequence in the bottom indicates a typical reconfiguration event during which weak bundles (dashed) become strong ones (solid) and vice-versa. (II) Nearly static pattern obtained increasing , (III) which becomes highly dynamic by further increasing . Time sequence in the bottom indicates illustrate the collapse (black), expansion (purple) and splitting (green) of cells in the network. See Movie 6. We indicate by * dynamical patterns exhibiting spatiotemporal chaos.

At longer times, the nematic bundles in Fig. 4a(III) develop secondary active instabilities leading to coordinated bending, curvature amplification, defect nucleation and annihilation, Movie 4(III), in a behavior reminiscent of active extensile systems [3436]. This possibility is far from obvious because, even though the deviatoric part of active tension is extensile since κ < 0, total active tension is contractile since |κ| < 1, Fig. 1c. To systematically examine this point, we performed simulations at higher active tension anisotropies κ = −0.8, which we compared with our reference κ = −0.2, Fig. 4b. While in our reference system bundles behave like contractile objects tending to straighten, for κ = −0.8 they behave like extensile objects enhancing curvature (see blue insets), which results in continuous defect nucleation as highly bent bundles destabilize and fragment, as well as defect annihilation as pairs of free ends merge to reorganize the network (see second inset in (b,II)), in a behavior akin to active turbulence [85]. See Movie 5 for an illustration and for the slightly extensile case κ = −0.5. To further substantiate this contractile vs extensile behavior of the nematic structures emerging from a contractile active gel with extensile deviatoric behavior, we plotted the difference between the stress along the nematic direction and perpendicular to it, σ||σ. We note that for κ < 0 as required to obtain fibrillar patterns, the active component of this quantity is negative, , Fig. 1c. In the case of bundles with contractile phenomenology, i.e. a tendency to straighten (b,I), we found that σ||σ was positive on the dense bundles and nearly zero in between, indicating that bundles are contractile structures embedded in a largely isotropic matrix. Conversely, in bundles with extensile phenomenology (b,II), we found intricate stress patterns with regions of negative σ||σ. The extensile tension pattern for κ = −0.8 is more clearly appreciated in 1D simulations, Supp. Fig. 2, where the dense nematic structures cannot fragment as a result of the bend-type instability. This figure also shows how the competition between active and viscous stresses dictates the emergent contractile vs extensile behavior of the nematic bundles. In the contractile case (κ = −0.2), even if , the negative viscous stresses due to the self-reinforcing flows are significantly larger in magnitude perpendicular to the nematic direction, resulting in σ||σ. In summary, our theory predicts that a contractile nematic active gel can self-organize into fibrillar patterns with mesoscale bundles that are either contractile or extensile depending on the parameter regime. This complex interplay between the contractility vs extensibility of the base material and that of the mesoscale nematic pattern resonates with recent observations where tissues made of contractile cells behave collectively as extensile nematic systems [45, 86].

Focusing on contractile bundles, we then examined a different parameter regime known to trigger sustained pattern dynamics. For isotropic gels, previous work has shown that reducing friction triggers chaotic dynamics as the distance between high-density regions, ℓpattern = 2π/vcrit, becomes comparable or smaller than the hydrodynamic length scale [29], thus enabling their mechanical interaction. In a model devoid of orientational order, reducing friction is equivalent to increasing . In our model, however, we can either reduce γ, which in non-dimensional terms means increasing and . in concert, or increase while leaving all other non-dimensional parameters fixed. The first of these choices leads to dynamic and hierarchical networks with very dense and highly-aligned bundles, which coexist with perpendicular weak bundles with much lower density enrichment and ordering. These two families of bundles enclose cells of isotropic and sparse gel, Fig. 4c(I). Junctions where two or more dense bundles meet are very unfavorable and short-lived, but junctions of two dense and a weak bundle are much more stable. Because bundles are mechanically coupled, the network actively reorganizes in events where dense bundles become weak bundles and vice-versa, inset and Movie 6(I). We note that here σ|| > σ, Fig. 4c(I) and Supp. Fig. 2d, and hence the persistent dynamics are unrelated to the previously described behavior of extensile bundles.

The second choice to favor mechanical interaction of bundles, increasing only, leads to very different networks with high-density aster-like clusters interconnected by straight actin bundles. Because now is not particularly large, order is low at the core of these clusters, enabling high-valence networks where four bundles often meet at one cluster. For = 10, the network is stable and nearly crystalline, Fig. 4c(II), whereas for = 20, it becomes highly dynamical and pulsatile with frequent collapse of polygonal cells by fusion of neighboring actin clusters and their attached bundles (black polygons) and nucleation of new bundles within large low-density cells (dashed/solid green lines), Fig. 4c(III) and Movie 6(III). This architecture and dynamics resemble those of early C. elegans embryos [87], adherent epithelial cells treated with epidermal growth factor [27] and mouse embryonic stem cells [83]. Recent active gel models accounting for RhoA signaling develop similar pulsatile behaviors in 2D, but do not predict the orientational order of the spatio-temporal patterns of the actomyosin cortex [88].

In summary, our theory maps how effective parameters of the actin gel control the active self-organization of a uniform and isotropic gel into a pattern of high-density nematic bundles embedded in a low-density isotropic matrix, including the activity threshold, the bundle spacing, orientation, connectivity and dynamics.

Assessment of activity parameters κ and λ through discrete network simulations.

(a) Illustration of the computational domain of the discrete network as a uniform representative volume element of the gel. (b) Sketch of model ingredients and setup to compute tension along and perpendicular to the nematic direction. (c) Typical time-signal for parallel and perpendicular tensions following addition of crosslinkers and motors (translucent lines) along with time average (solid lines) for isotropic and anisotropic networks. Tension is normalized by mean tension computed from time-averages and time by actin turnover time. (d) Mean tension as a function of network density for several nematic parameters S0, where both quantities are normalized by their values for the lowest density. With this normalization, Eq. (11) predicts a linear dependence with slope λ = 1 (dashed line). Error bars span two standard deviations. (e) Deviatoric tension as a function of nematic order for different densities. The dashed line is a linear regression to simulation data. (f) Dynamics of nematic order in a periodic network following addition of crosslinkers and motors for three initial values of nematic order. (g) Rate of change of nematic order normalized by turnover rate as a function of initial nematic order at zero and finite temperature.

Microscopic origin of κ < 0 and λ > 0 through discrete network simulations

A central prediction of our model is that the self-organization of nematic bundles, the most prominent emerging organization in actin gels across cell types and length-scales, requires that active tension perpendicular to nematic orientation is larger than along this direction (κ < 0, see Fig. 1c), at least at the onset of pattern formation. This fact may seem counter-intuitive in that dense nematic bundles are associated with large contractile tension along their axis, but as discussed in Fig. 4c(I), even if κ < 0 bundles can be contractile because of the interplay between active and viscous stresses. Another prediction of our model is that the active assembly dense nematic bundles requires active alignment controlled by parameter λ > 0, Fig. 3. Being central to our conclusions, we then sought to examine the plausibility of the conditions κ < 0 and λ > 0 of our phenomenological theory by performing discrete network simulations using the open-source cytoskeleton simulation suite Cytosim [89].

Similarly to reconstituted actomyosin gels, discrete network simulations of the actin cytoskeleton tend to irreversibly collapse into clumps [70], although proper tuning of model parameters can lead to sustained contractile states [90]. However, to our knowledge, these models have not been able to reproduce the heterogeneous nonequilibrium contractile states involving sustained self-reinforcing flows underlying the pattern formation mechanism studied here. For this reason, rather than trying to reproduce the assembly and maintenance of patterns of dense nematic structures, we specifically focused on investigating whether the average behavior of the discrete network is compatible with the conditions κ < 0 and λ > 0 of our continuum theory. For this, we studied the buildup of tension and dynamics of nematic order in uniform representative volume elements of varying density and filament alignment after being brought out of equilibrium, Fig. 5a. Each of our simulated volume elements can represent a uniform system prior to pattern formation, or a material point in the actomyosin gel.

We performed 2D simulations in which semi-flexible filaments interact with crosslinkers and myosin motors, all of which undergo turnover and have a stoichiometry previously used to model the actin cytoskeleton [91], Fig. 5b. See Appendix D for a detailed description of the simulation protocol. Briefly, we modified Cytosim to account for average orientational order in the simulation box, which we evaluated as a sample average of orientations over the ensemble of segments composing the filaments. We further introduced a nematic energy penalty in the network allowing us to restrain average nematic order to a target value S0.

We first prepared a system consisting only of randomly oriented actin fibers, imposed the desired orientational order S0 using the nematic penalty and equilibrated the system. In a first set of simulations, once S0 was reached, we deactivated the nematic penalty and added crosslinkers and myosins, driving the system out-of-equilibrium. The free contraction of the system was prevented by the addition of anchors at the boundary of the representative volume element, which also allowed us to compute anchor forces and hence estimate the effective active tension along the nematic direction σ|| = σijσiσj and perpendicular to it, σ = σijσimj with nimi = 0 and mimi = 1, Fig. 5b.

Addition of crosslinkers and myosins leads to bundling of actin filaments at the microscale, Movie 7, distinct from the mesoscale fibrillar pattern formation emerging from the active gel model. It also leads to out-of-equilibrium tension as measured by the anchors. For isotropic networks (S0 = 0), active tension is isotropic with σ||σ. For anisotropic networks, however, we found that tension becomes anisotropic with σ > σ||, Fig. 5c.

We systematically characterized this behavior varying initial orientational order and network density. According to our active gel model, Eqs. (5,6,7), in the absence of nematic gradients and flow, the tension components σ|| and σ satisfy the following relations in terms of mean and deviatoric tensions

Remarkably, our discrete network simulations closely followed these relations, Fig. 5d,e, which allowed us to estimate κ ≈ −1.6. We robustly found that κ < 0 for perturbations of selected parameters of the discrete network model as long as turnover rates of crosslinkers and myosins were relatively fast.

We then wondered if the discrete network simulations could provide evidence for the orientational activity parameter in our theory, ρλ. For a uniform system with nematic order along a given direction and ignoring the susceptibility parameter b, Eq. (8) becomes

where ηrot is the viscous drag of the filaments in the discrete network simulations, a > 0 the entropic tendency of the model to return to isotropy, ρλ the active forcing of nematic order resulting from crosslinkers and motors, and the last term accounts for the effect of the nematic penalty with coefficient KS. As a first test of this model, we started from an isotropic and periodic network and tracked the athermal dynamics of S under the action of the nematic penalty in the absence of anchors, crosslinkers and myosins. For λ = 0 and a = 0, Eq. (12) predicts an exponential relaxation given by S(t) = S0(1 − e−KSt/ηrot), which very closely matched the simulation data for different values of So and for a single fitting parameter nrot, Supp. Fig.3(III) and Movie 8. We then deactivated the nematic penalty and added crosslinkers and motors, but not anchors, to track unconstrained dynamics of nematic order starting from different values of S0. In agreement with the notion of an active force driving nematic order, we found that S(t) monotonically increased, Fig. 5f. More quantitatively, we tested the short-time prediction of Eq. (12), , by plotting as estimated from our simulations, as a function of S0, Fig. 5g. We found a nearly linear relation with positive slope, hence providing evidence for an active generalized force driving order. In agreement with the theory, in the athermal limit, the tendency to actively increase order is faster as the entropic tendency to disorder is absent (a = 0).

In summary, discrete network cytoskeletal simulations provide a microscopic justification for two key ingredients of our active gel theory, namely that nematic order elicits (1) anisotropic active tensions, which can be larger perpendicular to the nematic direction (κ < 0), and (2) active generalized forces driving further ordering..

Conclusions

We have developed a theory for the active self-organization of initially uniform and isotropic actin gels into various localized dense nematic architectures embedded in an isotropic matrix of low density. This model predicts a variety of emergent patterns involving asters, tactoids, and bands with perpendicular nematic organization. More importantly, it identifies a wide parameter space where the active gel spontaneously develops patterns of dense nematic bundles, the most prominent nematic architecture across scales and cell types. Contrary to previous works on instabilities and pattern formation in active nematics, the mechanism proposed here relies on the advective instability and compressible self-reinforcing flows typical of actomyosin gels. We have characterized how the activity threshold, spacing, geometry, connectivity and dynamics these patterns depends on effective active gel parameters. Because these mesoscale parameters depend on the composition and dynamics of the network at the molecular scale, our results portray actin gels as responsive and reconfigurable active materials with an intrinsic ability to assemble patterns of nematic bundles that cells can finely regulate.

We have further shown that the spontaneous tendency of the gel to assemble bundle patterns can be directed via subtle cues. Thus, a combination of biochemical control of actin dynamics along with geometric, mechanical or biochemical guiding [92, 93] may explain the emergence and context-dependent organization of regular patterns of bundle networks, from sub-cellular to organism scales [12, 14, 19, 21, 22, 2528, 30]. Consistently, perturbations of myosin contractility have been shown to alter, disorder, or even prevent the formation patterns of parallel bundles in C. elegans [19], whereas perturbations of actin polymerization in drosophila embryos impair the robust organization of actin bundle patterns at the cellular and organ scales by disrupting orientation and spacing, but not the intrinsic tendency of the actomyosin cytoskeleton to form patterns of parallel bundles [21]. In adherent cells, the patterns of dense nematic bundles presented here may act as precursors of stress fibers as the actomyosin cytoskeleton interacts with the focal adhesion machinery in a cellular domain with boundary conditions set by the polymerization velocity at the leading edge and the nucleus.

Our theory identifies two key requirements on activity parameters for the self-organization of patterns of nematic bundles, namely active tension anisotropy with larger tension perpendicular to the nematic direction and generalized active forces tending to increase nematic order. Complementarily to our phenomenological theory, we have examined the plausibility of these conditions with discrete network simulations of homogeneous representative volume elements of different density and orientational order, which have confirmed our constitutive assumptions. We expect, however, that in a different parameter regime anisotropic tension may be larger along the nematic direction (κ > 0). For instance, once bundles are dense and maximally aligned, the ability of the active nematic gel to perform active tension perpendicular to the nematic direction may saturate, while myosin motors may contract the gel along the nematic direction more effectively. The regime studied here explains the initial assembly of dense nematic bundles, but not their maturation to become highly contractile or viscoelastic as demonstrated for stress fibers depending on different isoforms of nonmuscle myosin II or on actin regulators such as zyxin [94, 95]. Our work thus suggests further experimental and computational work to establish a comprehensive mapping between molecular and mesoscale properties of the active gel, and how these properties control the emergent network architecture and mechanical properties.

Lengthscale of the pattern.

Illustration of the lengthscale of the pattern close to the onset of the instability (left) and deeper into de nonlinear regime (right). The linear stability estimate of the lengthscale of the pattern is ℓpattern = 2π/vcrit, where Vcrit is given by Eq. (10) in terms of material parameters. The side length of the square periodic computational domain is 8ℓpattern. Close to the onset (left), the pattern along the black arrow, with exactly 8 repeats, shows that the lengthscale in the nonlinear simulations closely follows the theoretical estimate. At activities 30% above the threshold (right), we sampled the typical separation between dense structures (green arrows) to approximate from simulations, finding that (mean and standard deviation). We robustly found across the parameter space that the pattern lengthscale from nonlinear simulations 30% above the threshold was a little smaller than the theoretical estimate ‘pattern, see density maps in Figs. 2 and 4.

Tension distribution along (σ||) and perpendicular to (σ) the dense nematic bundles.

The left column shows the density distribution, the middle column the total tension along (solid) and perpendicular (dashed) to nematic bundles, and the right column the different contributions to the total tension dominated by the active (black) and the viscous (blue) components. In all cases, we consider for convenience fully nonlinear simulations in 1D to easily define the orthogonal directions relative to the self-organized pattern. (a) to (c) show patterns obtained for negative tension anisotropy κ of increasing magnitude and correspond to Fig. 3b in the main text, whereas (d) shows a chaotic pattern resulting from a large hydrodynamic length and corresponds to Fig. 3c(I) in the main text. The model parameters used in the plots are the same as in Figs. 2 and 3 and are described in Tables I and II

Discrete network simulations.

(I) Microstructural modeling approach using cytosim. (a) The nematic ordering of the network, S, measured with respect to a director , is controlled by a system-wide restraining energy. (b) Actin filaments are represented by connected points kept at a fixed distance and have finite bending rigidity κA. Crosslinkers are Hookean springs with stiffness ℓX and resting length KX. Myosin motors are Hookean springs with stiffness M and resting length M; their ends can walk on filaments at a speed υM, which is affected by force application. (II) Simulation protocol to quantify active tension anisotropy. (a) Initial fiber seeding without orientational bias. (b) Athermal imposition of desired order parameter S0 using a restraining potential. (c) Introduction of boundary anchors. (d) Deactivation of orientational restraining potential and addition of crosslinkers and myosin motors, which drive the system out-of-equilibrium. Reactions at anchors allow us to track tensions at the boundary, (e,f). (III) Protocol to quantify the orientational activity parameter ρλ. (a) Initial fiber seeding without orientational bias in a periodic box. (b) Athermal imposition of desired order parameter S0 by restraining potential. (c,d) Activation of temperature, deactivation of orientational restraining potential, and addition of crosslinkers and myosin motors. The dynamics of orientational order of the system driven out-of-equilibrium are then tracked during a time period (c,d). (e) Estimation of the model parameter ηrot, which characterizes the athermal dynamics of network reorientation in the presence of a restraining potential of stiffness KS; note that a single value of ηrot fits the dynamics for all considered S0. (f) The rate of change of nematic order, Ṡ, is estimated from a 10 s time interval.

Model parameters used in figures.

Model parameters used in movies that do not directly reproduce figures.

Global parameters adopted in this study and in previous microstructural models that used cytosim.

Actin filament parameters adopted in this study and in previous microstructural models that used cytosim.

Myosin motor parameters adopted in this study and in previous microstructural models that used cytosim.

Acknowledgements

The authors acknowledge the support of the European Research Council (CoG-681434) and the Spanish Ministry for Science and Innovation (PID2019-110949GB-I00). WM acknowledges the La Caixa Fellowship and the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie action (GA 713637). MP acknowledges the support from the Spanish Ministry of Science and Innovation & NextGenerationEU/PRTR (PCI2021-122049-2B). MA acknowledges the Generalitat de Catalunya (ICREA Academia prize for excellence in research). MDC acknowledges funding from the Spanish Ministry for Science and Innovation through the Juan de la Cierva Incorporation fellowship IJC2018-035270-I. IBEC and CIMNE are recipients of a Severo Ochoa Award of Excellence.

Additional information

Author Contributions

WM, ATS and MA conceived the study and developed the active gel theory. WM, GV and ATS developed and implemented numerical methods. WM performed and analyzed active gel simulations. WM and MP performed and analyzed discrete network simulations. MP developed the code and simulation protocols for discrete network simulations. MDC contributed to linear stability analysis. WM and MA wrote the manuscript. ATS and MA supervised the study.

Additional files

Movie 1: Pattern formation in an active gel model not accounting for nematic order.

Movie 2: Pattern formation for different values of the tension anisotropy coefficient κ, leading to (I) fibrillar patterns for κ ≈ 0, (II) tactoids for κ « 0, and (III) banded patterns for κ > 0. This movie corresponds to Fig. 2b.

Movie 3: Effect of flow-alignment coupling coefficient β on pattern formation for positive and negative tension anisotropy coefficient κ. Flow-alignment favors fibrillar patterns and disfavors banded patterns.

Movie 4: Effect of orientational bias on fibrillar patterns. This movie corresponds to Fig. 3a.

Movie 5: Contractile vs extensile bundles as the magnitude of κ increases, leading to active turbulence of fibrillar patterns. This movie correspond to Fig. 3b.

Movie 6: Changes in network geometry and dynamics following enhanced inter-bundle mechanical communication. This movie correspond to Fig. 3c.

Movie 7: Illustration of simulation protocol to estimate tension anisotropy from discrete network simulations for an isotropic (left) and anisotropic (right) representative volume elements. The phases are (1) alignment of filaments to reach the desired nematic order, (2) addition of boundary anchors (green dots) to prevent network collapse and to measure tension, (3) addition of crosslinkers and motors (blue and red dots) to bring the system out-of-equilibrium while tracking tension along horizontal and vertical directions.

Movie 8: Illustration of simulation protocol to estimate the generalized active tension conjugate to nematic order for an isotropic (left) and anisotropic (right) representative volume elements with periodic boundary conditions. The phases are (1) alignment of filaments to reach desired nematic order, (2) addition of crosslinkers and motors (blue and red dots) to bring the system out-of-equilibrium while tracking the average nematic order of the representative volume element.

Author response video 1

Appendix A: 1D reduced model

To perform the linear stability analysis, we reduce the theoretical model to 1D by considering the flow velocity v(x,t) and the density field as ρ(x, t) along the x-axis. The nematic order tensor is defined as qij = S(ninjδij/2), where n is the local average filament orientation and S is the local degree of alignment. Since it is traceless and symmetric, in general it can be represented by 2 independent degrees of freedom, q11 = −q22 = q1 and q12 = q21 = q2. In the 1D model we assume that n is either along or perpendicular to x, and hence q12 = q21 = 0. We are thus left with one independent degree of freedom q11 = −q22 = q. Positive (negative) q represents alignment along (perpendicular to) the x-axis. The Jaumann derivative of the nematic order parameter reduces in 1D to

Particularizing the general equations given in given in the main text, we present next the 1D governing equations pertinent to the linear stability analysis. Mass conservation reads

balance of linear momentum of the active gel reduces to

where stress along x, σ, is given by

The evolution of nematic order q is governed in the 1D setting by

For c0 = 2aρ0λ⊙ ≥ 0, the uniform steady state of the system is given by ρ(x,t) = ρ0, v(x,t) = 0 and q(x,t) = 0. For c0 < 0, the uniform steady state has spontaneous alignment given by . Finally, we note that, for the dissipation to be positive, the material parameters need to satisfy an entropy production inequality 4ηηrotβ2 ≥ 0 in the reduced one-dimensional model and 2ηηrotβ2 ≥ 0 in the two-dimensional model [64].

Appendix B: Non-dimensionalization

We present next a non-dimensionalization of the governing Eqs. (A2A5). We choose the hydrodynamical lengthscale , above (below) which friction (viscosity) is the dominant dissipative mechanism in the active gel, and the time-scale of diffusion over this length-scale, . Using these characteristic scales, we non-dimensionalize space and time as and . We chose ρ0 as a reference density and thus . Here and elsewhere, overbar denotes non-dimensional quantities.

The dimensionless balance of mass reads

where . The dimensionless balance of linear momentum reads

where and the characteristic stress is given by σ0 = ρ0γD. The non-dimensional constitutive relation for the stress is

where , , and . We note that κ = κaniso/κ and q are already nondimensional. Finally, the dimensionless equation of generalized force conjugate to nematic order is given by

where , , and . Inspection of the Eq. (B4) reveals the characteristic nematic length-scale and the nematic relaxation time-scale .

Appendix C: Linear stability analysis

We perform next the linear stability analysis on the 1D model given by the dimensionless Eqs. (B1B4). For the sake of simplicity of notation, we omit overbars for the dimensionless quantities. We assess the stability of the homogeneous state given by v(x, t) = 0, q(x, t) = q0 and ρ(x, t) = ρ0 by examining the fate of infinitesimal perturbations v (x, t) = δv (x, t), q (x, t) = q0 + δq (x, t) and ρ (x, t) = ρ0 + δρ (x, t), where |δv| ~ |δq| ~ |δρ| ≪ 1. The control parameter in our stability analysis is the activity parameter κ.

Inserting the infinitesimally perturbed homogeneous states into the governing Eqs. (B1B4) and ignoring terms that are quadratic in perturbations, the resulting balance of linear momentum equation is

The linearized form of Eq. (B4) is

If c0 = 2a − ρ0κ > 0, then q0 = 0 and G = c0. If c0 = 2a − ρ0κ < 0, then and . The linearized balance of mass reads

We decompose perturbations into a sum of Fourier terms, i.e., δv = v eαt+ixv, δρ = ρ eαt+ixv and δq = q eαt+ixv, where v′, q′ and p′ are the amplitudes of the perturbations, α is the growth rate and υ is the wavenumber.

Substituting the Fourier representation in Eq. (C1), we obtain

Likewise, Eq. (C3) becomes

Substituting this equation into the Fourier representation of Eq. (C2) we obtain

Finally, combining the three equations above, we obtain the dispersion equation, which, given the material parameters, relates wavenumber υ to growth rate a and takes the form

where

The roots of Eq. (C7) are given by . The system is linearly stable if perturbations decay, i.e. Real (α) < 0, marginally stable if Real (α) = 0, and unstable if Real (α) > 0. We found that for all reasonable material parameters α(υ) is real, and hence the critical condition is given by α = 0, which recalling Eq. (C7) implies that C(υ) = 0. From this condition, we can express the activity parameter as

To find the critical wavenumber vcrit and the associated critical activity parameter κcrit = κ(vcrjt), we minimize this expression with respect to v, or equivalently with respect to x = v2, which requires that

1. Case A: q0 = 0 and L small

When c0 > 0, the spontaneous order in the uniform steady-state vanishes and we have

Given a set of parameters, this equation can be numerically minimized to obtain the critical wavenumber and activity. To obtain workable explicit expressions, we first assume that L = 0, leading to

where

From the condition 0 = ∂κ/∂x we obtain

In the limit of δ → 0 we recover for the largest root the results and of an active gel model without nematic order [29]. Taylor-expanding the equation above in terms of δ, we obtain

and

We note that these expressions are non-dimensional. The dimensional expressions are reported in the main text.

2. Case B: q0 = 0 and c0 small

Using Eq. (C10), we examine the situation in which c0 approaches 0 but is positive so that q0 = 0. Eq. (C10) then becomes

where now

Following the same rationale as before, we find

Expanding up to linear order in δ, we obtain

and

3. Case C: q0 ≠ 0 and L small

Going back to Eq. (C8) and taking L = 0 to obtain workable expressions, we obtain

where

Minimization of κ with respect to x leads to

valid when 1 + τ > 0. Taylor-expanding the equation above in terms of δ we obtain

and

Appendix D: Choice of active gel model parameters in the main text

We consider reference model parameters for the actin cytoskeleton and variations about these values since they can be expected to change significantly between cells and conditions. We take a cortical viscosity of η = 104 Pa s, and hydrodynamic lengths in the order ofls = 0.2icrom [29, 33], although we vary this parameter, as indicated in Tables I and II. Taking kd ~ 0.1s−1 [26, 29, 96] and D ~ 0.04 μm2 s−1 [26, 29], we obtain Damkölher lengths lD in the order of a micron. We consider ηrot/η ~ 1 and β/η ~ −0.2, far from the threshold given by the entropy production inequality. We view κ as the master activity parameter and take values a bit higher than the threshold for pattern formation, which depends on other material parameters as shown in our expressions for κcrit. We vary the nondimensional active tension anisotropy parameter κ. We choose the susceptibility parameters as arot = 5 s−1 and b/ηrot = 20 s−1 so that relaxation of nematic order is significantly faster than the turnover time-scale and choose the Frank constant to be small so that lq is small compared to other lengthscales but not too small so that our computational grid can resolve it. The full list of material parameters for each figure and movie is given in the Tables I and II.

Appendix E: Discrete network simulations

To ascertain the microstructural origin of the anisotropic active tensions underlying the formation of stress fibers, and provide evidence for the orientational activity parameter (ρκ) included in our theory, we establish an agent-based microscopic model of a crosslinked actomyosin network using the open source cytoskeletal simulation suite cytosim [89, 97]. We customized the source code to impose and track nematic order in the system [98]. Below we provide a concise description of the modeling approach and an overview of the simulation parameters.

1. Model description

Actin filaments are represented through a Langevin equation that recreates the Brownian motion of a thermal system with temperature T (thermal energy: kT, where k = 1.38 × 10−5 pN μm K−1 is Boltzmann’s constant) and includes bending elasticity as well as external forces (e.g. those applied by crosslinkers and myosin motors). The filaments are immersed in a medium of viscosity υ and the motion of their points is described by

where F (X , t) contains the forces acting on the points at time t, dB(t) is a stochastic non-differentiable function of time summarizing the random molecular collisions that lead to Brownian motions, and the matrix μ contains the mobility coefficients of the particles that constitute the system.

The simulated system comprises NA actin filaments, which are modeled as inextensible objects with bending rigidity κA and individual filament length lA; both quantities are uniform across the system. Following a previously-proposed approach [99], we represent actin turnover by randomly selecting and completely removing Nt filaments every Ns time steps and replacing them with an equal number of new ones that are placed at randomly-chosen locations and have the same orientations as the removed filaments. The values of Nt and Ns are chosen to ensure that a given fraction of the total filaments in the system be replaced per unit of time, according to a turnover rate rA = Nt/(NsNA).

To control the nematic characteristics of the filament network, we introduce a system-wide restraining energy of stiffness κS that penalizes deviations from a target orientational order parameter S0 , measured with respect to a director, , that is aligned with the Cartesian basis vector , see Supp. Fig. 3I(a),

In this equation, we evaluate the average nematic order of the simulation box as a sample average given as over the ensemble of N segments that constitute the filament network and mI denotes the unit vector corresponding to the I-th segment in the system. The energy defined by Eq. (E2) results in restraining forces that act on the particles comprising the actin filaments, such that a particle of position XP is subjected to a force

In order to write the derivative qs(X)/∂XP explicitly, we denote with mP,− = (XPXP−)/|XPXP−| and mP,+ = (XP+1XP)/|XP+1XP the unit vectors that correspond to the filament segments extending from point XP towards the head (+) and the tail (−) of the filament, see Supp. Fig. 3I(b). In the absence of filament bifurcations, which are not present in our model, mP,+ and mP,− are the only unit vectors in the network that depend on XP . Thus, it is immediate to simplify the derivative on the right hand side of Eq. (E3) and express the restraining force acting on the P-th particle as

We shall also note that Eq. (E4) is simplified when particle P coincides with the head or the tail of a filament, since either mP,+ or mP,− will not exist. For a filament head, the restraining energy introduced in Eq. (E2) will thus result in the force

Similarly, the force acting at a filament tail will be

Finally, assuming that each actin filament is represented by segments of constant length sA, the derivatives of the unit vectors appearing in Eqs. (E4E6) are obtained directly from the definitions of mP,−, and mP,+ as

Crosslinkers and myosin motors are modeled as Hookean springs with finite stiffness KX and KM and resting lengths lX and lM, respectively. For simplicity, the concentration of unattached species is assumed to be uniform across the modeling space and their diffusion is thus not simulated explicitly. The ends of the springs can bind to any discrete location along the filament segments, as long as they fall within finite binding ranges denoted as and . When multiple filament locations fall within the binding range, the springs attach to the closest point on the filament and apply a force that depends linearly on their elongation (uX and uM)

Binding of the species is modeled as a purely stochastic event whose time of occurrence follows an exponential distribution. The expected values for crosslinker and motor binding are and , and being the corresponding binding rates. Unbinding from a filament is not always purely stochastic but follows a rate that can increase when the springs are loaded, according to a relationship known as Bell’s law [100]. Denoting the base unbinding rates for crosslinkers and motors as and , the effective rates under a force of magnitude f are expressed as

where and are constant parameters associated with the bound state. A specific feature of motors is that they can ‘walk’ on actin filaments by displacing their attachment points without necessarily having to unbind; this effect can be modified by load application. Denoting with the speed of the motor in the absence of external loads, the effective speed is given by

where is the stall force, which controls the slope of the speed-force relationship. By convention, we consider that positive (negative) speed values correspond to a motion towards the head (tail) of a filament.

2. Quantification of network-scale active tensions

To determine anisotropic active tensions arising when the network is driven out-of-equilibrium, we examine a square representative volume element of crosslinked actomyosin network at room temperature (kT = 0.0042 pN μm). The modeled systems feature an edge length 5 μm, a filament surface density ρA, a crosslinker density ρX = 15.4 (actin μm)−1, and a myosin motor density ρM = 0.8 (actin μm)−1. The actin density ρA is varied throughout the study and takes the following values: ρA = (78; 156; 234) (actin μm)/μm2, which we denote as ρ, 2ρ, and 3ρ. All simulation parameters adopted in this study are indicated in Tables IIIV, where we also provide an overview of the values used in previously-proposed microstructural models that are comparable to the ones developed here. Throughout the simulations, we adopt a time step of 5 ms and a data acquisition frequency of 40 s−1, i.e. 1 every 5 time steps is written to output for further data processing.

The filaments are initially seeded uniformly throughout the modeling space and without any orientational bias, Supp. Fig. 3II(a), followed by equilibration in the presence of the restraining energy defined by Eq. (E2), which acts with a stiffness κS = 5000 pN/μm to induce a nematic orientation characterized by the ordering parameter S0 , see Supp. Fig. 3II(b). No crosslinkers or myosin motors are present at this stage. After 5 s, Supp. Fig. 3II(c), the filaments are clamped to the system’s boundary using NH anchors. These objects, whose spatial position is fixed, are akin to crosslinkers but have spring stiffness KH = 5000 pN/μm, zero resting length, a binding range = 5 nm, and a binding rate = 100 s−1. Anchor unbinding is completely hindered by setting = 0 s−1 and . The quantity NH is defined by the number of filaments that cross the system’s boundaries. Indicating with , , , and the number of anchoring objects located on the left, right, bottom, and top edge of the system, it follows that .

After having equilibrated the network for an additional 2.5 s, we deactivate the restraining potential by setting κS = 0 and introduce crosslinkers and motors in the system, Supp. Fig. 3II(d), which we let interact with the filaments for 12.5 s in order to drive the system out-of-equilibrium. This results in reaction forces being applied at the anchoring points. At each time step, the total average tensions acting on the system edges that are oriented parallel and perpendicular to the nematic director , σ|| and σ, are measured as

Measuring tensions from, t = 10 s, Supp. Fig. 3II(e), to t = 20 s, Supp. Fig. 3II(f), allows us to quantify tension anisotropy using the deviatoric tension normalized by the mean tension

where < · > denotes the time average over the considered period of interest. Finally, quantifying the above ratio for 12 sets of n =16 model realizations that correspond to values of S0 = (0.0; 0.1; 0.2; 0.3) and ρA = (ρ; 2ρ; 3ρ) allows us to determine the mesoscale activity parameter κ by fitting the ξ versus S0 relation, see Fig. 5e. To this end, we leverage the function ‘stats.linregress’ that is included in the Python package ‘SciPy’ [101].

3. Quantification of the orientational activity parameter

To provide evidence for the orientational activity parameter in our theory (ρλ), we examine the behavior of a periodic square representative volume element of crosslinked actomyosin driven out-of-equilibrium. The absence of anchors in these systems allows us to track the unconstrained nematic order dynamics. We limit our analysis to models with actin density ρA = 2ρ and adopt the same simulation parameters as in Section E 2, unless explicitly stated otherwise.

To highlight the contribution of the restraining potential in Eq. (E2), we initially focus on athermal systems (κT ≈ 0). After seeding NA actin filaments uniformly throughout the modeling space and without any orientational bias, Supp. Fig. 3III(a), we let the restraining potential reorient the system for 2.5 s in the absence of crosslinkers and myosin motors to reach an ordering parameter, S0 Supp. Fig. 3III(b). Following the first 2.5 s of athermal simulation, we increase the temperature to its standard value (kT = 0.0042 pN μm), deactivate the restraining potential by setting κS = 0, and add crosslinkers and myosin motors to the system, which we let interact with the filaments for 30 s in order to drive the system out-of-equilibrium, Supp. Fig. 3III(c,d). We observe that S(t) increases monotonically, allowing us to estimate according to a linear least-squares fit, as implemented in the Python ‘SciPy’ function ‘stats.linregress’ [101]. We notice that in Supp. Fig. 3III(a,b) the time evolution of S is approximately exponential, Supp. Fig. 3III(e), as predicted in the main text. This allows us to leverage our discrete network simulations to estimate the model parameter ηrot ≈ 1583 pN μm s using the nonlinear least-squares approach provided by the function ‘optimize.curvefit’ that is included in the Python package ‘SciPy’ [101]. Note that the value of ηrot is determined by fitting one single evolution curve for S(t)/S0, which is obtained by averaging the system dynamics resulting from the 4 sets of n = 16 model realizations that correspond to values of S0 = (0.1; 0.2; 0.3; 0.4).

Finally, we plot our simulation-based approximation of for several values of S0 , which results in an almost linear relationship with positive slope, implying that ρλ > 2a (see Fig. 5g). Moreover, considering fully athermal simulations (i.e. simulations in which the temperature is not increased after reaching the ordering parameter S0) results in an linear relation between and S0 with slightly larger slope. This finding agrees with the predictions of our theory for a > 0 and demonstrates that temperature, which entropically drives the network towards isotropy, opposes the effect of the orientational activity parameter ρλ.