Theory of active self-organization of dense nematic structures in the actin cytoskeleton

The actomyosin cytoskeleton can be understood as a compressible, active, and viscous fluid gel with orientational order undergoing turnover (Salbreux et al., 2012; Balasubramaniam et al., 2022). 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 (Hannezo et al., 2015; Bois et al., 2011; Kumar et al., 2014; Callan-Jones and Voituriez, 2013). This kind of advective instability has been invoked to explain cell polarization and amoeboid motility (Ruprecht et al., 2015; Callan-Jones and Voituriez, 2013; Bergert et al., 2015), the formation of the cytokinetic ring (Mietke et al., 2019), the formation of periodic dense actin structures during tracheal morphogenesis in Drosophila (Hannezo et al., 2015), or the emergence of self-sustained dynamical states in actomyosin gels extracted from cells (Krishna et al., 2024; Malik-Garbi et al., 2019), and has been recently reproduced to some degree in confined reconstituted gels from purified proteins (Sciortino et al., 2025). In all of these examples, the actomyosin gel develops sustained compressible flows converging toward regions of high density. Furthermore, the compressive strain rate induced by these active flows has been shown to drive nematic order (Reymann et al., 2016; Salbreux et al., 2009). Finally, observations on adherent cells show that active contractility is required for actin bundle formation (Wirshing and Cram, 2017; Lehtimäki et al., 2021; Hotulainen and Lappalainen, 2006). 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 (Lehtimäki et al., 2021; Vignaud et al., 2021), rather than thermodynamically enforcing high nematic order everywhere except at topological defects (Doostmohammadi et al., 2018; Gennes and Prost, 1993; Beris and Edwards, 1994; Soares e Silva et al., 2011).

Previous models for dry and dilute aligning active matter develop density patterns (Zumdieck et al., 2005; Chaté, 2020; Putzig et al., 2016) but fail to capture the hydrodynamic interactions of actomyosin gels, whereas models for active nematic fluids either ignore density (Santhosh et al., 2020; Marenduzzo et al., 2007; Giomi, 2015; Jülicher et al., 2018; Pearce, 2020; Metselaar et al., 2019; Srivastava et al., 2016; Pokawanvit et al., 2022) or account for the density of active particles suspended in an incompressible flow (Aditi Simha and Ramaswamy, 2002; Ramaswamy and Rao, 2007; Hatwalne et al., 2004; Giomi et al., 2014), 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 ignore either hydrodynamics (Zumdieck et al., 2005) or the density and flow compressibility characteristic of actomyosin gels (Santhosh et al., 2020), 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 Mirza et al., 2025, 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 Jülicher et al., 2018, or as a nematic generalization of the isotropic theory used in Hannezo et al., 2015. As elaborated in Mirza et al., 2025, it is possible to develop an alternative compressible and density-dependent active nematic theory based on the Beris-Edwards formalism (Santhosh et al., 2020; Beris and Edwards, 1994; Marenduzzo et al., 2007; Giomi, 2015; Pearce, 2020; Metselaar et al., 2019; Giomi et al., 2014).

In our model, the local state of the system is described by the areal density of cytoskeletal material $\rho (t,\mathit{x})$ and by the network architecture given by the symmetric and traceless nematic tensor $\mathit{q}(t,\mathit{x})$ (see Figure 1a), both of which depend on time t and position $\mathit{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 $q_{ij}=S\left(n_i{n}_j-{\delta }_{ij}/2\right)$, where $\mathit{n}$ is the average molecular alignment, $S=\sqrt{2q_{ij}{q}_{ij}}$ the degree of local alignment about $\mathit{n}$, and ${\delta }_{ij}$ is the identity. We denote by $\mathit{v}(t,\mathit{x})$ the velocity field of the gel. The rate-of-deformation tensor $\mathit{d}=\frac{1}{2}(\mathrm{\nabla }\mathit{v}+\mathrm{\nabla }{\mathit{v}}^T)$ measures the local rate of distortion of the fluid, whereas $\mathit{w}=\frac{1}{2}(\mathrm{\nabla }\mathit{v}-\mathrm{\nabla }{\mathit{v}}^T)$ measures its local spin, where ∇ is the nabla differential operator. The deviatoric part of the rate-of-deformation tensor is defined by $d_{ij}^{\mathrm{d}\mathrm{e}\mathrm{v}}=d_{ij}-(d_{kk}/2){\delta }_{ij}$. The rate of change of q relative to a frame that translates and locally rotates with the flow generated by $\mathit{v}$ is given by the Jaumann derivative $\hat{\mathit{q}}=\mathrm{\partial }\mathit{q}/\mathrm{\partial }t+\mathit{v}\cdot \mathrm{\nabla }\mathit{q}-\mathit{w}\cdot \mathit{q}+\mathit{q}\cdot \mathit{w}$ (Gennes and Prost, 1993).

Figure 1

Download asset Open asset

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 (Hannezo et al., 2015)

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 ${\rho }_0$ is the steady-state areal density and $k_d$ is the depolymerization rate. Note that for a uniform, steady-state, and quiescent gel, the first three terms vanish and $\rho (t,\mathit{x})={\rho }_0$.

Force balance in the gel can be expressed as

where $\gamma >0$ models a viscous drag with the surroundings (e.g. the plasma membrane) and $\mathit{\sigma }={\mathit{\sigma }}^{\mathrm{n}\mathrm{e}\mathrm{m}}+{\mathit{\sigma }}^{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{s}}+{\mathit{\sigma }}^{\mathrm{a}\mathrm{c}\mathrm{t}}$ 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 ${\mathit{\sigma }}^{\mathrm{n}\mathrm{e}\mathrm{m}}$ follows from a standard derivation adapted here to a density-dependent material. It derives from the free energy $\mathcal{F}={\int }_{\mathrm{\Omega }}\rho f(\mathit{q},\mathrm{\nabla }\mathit{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 (Mirza et al., 2025)

and to the explicit form of the nematic nonsymmetric stress tensor

The dissipative part of the stress

includes a viscous stress controlled by the gel viscosity parameter η (Salbreux et al., 2009) and a stress resulting from changes in the nematic field controlled by the coupling parameter β < 0 (Reymann et al., 2016). 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

The activity parameter λ controls the isotropic tension and is contractile for λ > 0, as assumed here. The additional activity parameter ${\lambda }_{\mathrm{a}\mathrm{n}\mathrm{i}\mathrm{s}\mathrm{o}}$, or equivalently $\kappa ={\lambda }_{\mathrm{a}\mathrm{n}\mathrm{i}\mathrm{s}\mathrm{o}}/\lambda$, 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 $|\kappa |<1$, then the total active tension remains positive in all directions, with a larger magnitude parallel to the nematic direction when $\kappa >0$ (contractile deviatoric component) and perpendicular to it when $\kappa <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\approx 0$), active tension is isotropic (Figure 1b), whereas when order is high, active tension becomes anisotropic (Figure 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 out-of-equilibrium binding of bundling proteins or myosins (Blanchoin et al., 2014; Harris et al., 2006; Courson and Rock, 2010; Schuppler et al., 2016; Li et al., 2017; Nandi, 2018; Ennomani et al., 2016; Chen et al., 2020).

Balance of the generalized forces power conjugate to $\hat{\mathit{q}}$ also includes viscous, elastic-nematic, and active contributions and takes the form

In this expression, ${\eta }_{\text{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 (Figure 1e), as experimentally observed in Reymann et al., 2016. This term involves the same coefficient β as the last term in Equation 6 because of Onsager’s reciprocity relations, and the entropy production inequality requires that $2\eta {\eta }_{\text{rot}}-{\beta }^2\ge 0$ (Mirza et al., 2025). The third term is a thermodynamic force driven by $\mathcal{F}$. In agreement with the observations that nematic ordering in actomyosin gels is actively driven, we assume a > 0, and therefore Equation 4, Equation 8 show that this term tends to restore isotropy. The last term is an active generalized force controlled by the activity parameter ${\lambda }_{\odot }\ge 0$ tending to further align filaments (Figure 1d; Reymann et al., 2016). This term is linear in $\rho$ because in the expansion ${\overline{\lambda }}_{\odot }+\rho {\lambda }_{\odot }$, the constant contribution ${\overline{\lambda }}_{\odot }$ can be subsumed by the susceptibility parameter a (Salbreux et al., 2009). Thus, the active term acts as a negative density-dependent susceptibility. When $c_0=2a-{\rho }_0{\lambda }_{\odot }<0$, the system can sustain a uniform quiescent state with $\rho (\mathit{x},t)={\rho }_0$, $\mathit{v}(\mathit{x},t)=0$ and a nonzero nematic order parameter $S_0^2=-c_0/b$. Even if $c_0>0$, and hence the uniform quiescent state is devoid of order, pattern formation can induce density variations such that $2a-\rho {\lambda }_{\odot }$ becomes negative locally and actively favors local nematic order. Physically, the term $-\rho {\lambda }_{⨀}\mathit{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 in Equation 1, Equation 7, Equation 8, along with the constitutive relations in Equation 4, Equation 5, Equation 6, Equation 7, has different sources summarized below. The theory presents nonlinearities intrinsic to transport equations in the advective term of Equation 1 and in the definition of the Jaumann derivative of the nematic tensor. Furthermore, nonlinearities in q in the constitutive relations result 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 (Mirza et al., 2025) result in further nonlinearities involving density in the constitutive relations. Finally, the nonlinearity involving density and the nematic field in the last term of Equation 8 has been discussed in the previous paragraph.

Our theory has three active parameters, λ, $\kappa$ and ${\lambda }_{⨀}$, 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 (Chugh, 2017; Koenderink and Paluch, 2018), e.g., the stoichiometry of filaments, crosslinkers, 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 (Hannezo et al., 2015). In the present model, however, translational, orientational, and density dynamics are intimately coupled through the terms involving β, ${\lambda }_{\odot }$, $\kappa$ and L.

We readily identify the hydrodynamic length ${\ell }_s=\sqrt{\eta /\gamma }$, above which friction dominates over viscosity, the Damköhler length ${\ell }_D=\sqrt{D/k_d}$ above which reactions dominate over diffusion, and the nematic length ${\ell }_q=\sqrt{L/\left|2a-{\lambda }_{\odot }{\rho }_0\right|}$. Nondimensional analysis reveals a set of nondimensional groups that control the system behavior, namely the nondimensional turnover rate ${\overline{k}}_d={\ell }_s^2/{\ell }_D^2$, the Frank constant $\overline{L}=L/(\eta D)$, the susceptibility parameters $\overline{a}=a/(\gamma D)$ and $\overline{b}=b/(\gamma D)$, the drag coefficients ${\overline{\eta }}_{\mathrm{r}\mathrm{o}\mathrm{t}}={\eta }_{\mathrm{r}\mathrm{o}\mathrm{t}}/\eta$ and $\overline{\beta }=\beta /\eta$, the nematic activity coefficient ${\overline{\lambda }}_{\odot }={\rho }_0{\lambda }_{\odot }/(\gamma D)$, and the active tension parameters $\overline{\lambda }=\lambda /(\gamma D)$ and κ (Appendix B). The full list of material parameters for each figure is given in Appendix 4—table 1 and Appendix 4—table 2 and justified in Appendix D .

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)$, $\rho (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 $c_0=2a-{\rho }_0{\lambda }_{\odot }>0$ to examine the loss of stability of a uniform, isotropic, and quiescent steady state ($\rho (x,t)={\rho }_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 ${\lambda }_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t},0}={\left(\sqrt{\gamma D}+2\sqrt{k_d\eta }\right)}^2=\gamma D(1+2\sqrt{{\ell }_s/{\ell }_D}{)}^2$ and $\delta =\gamma D\kappa \beta /(2\eta c_0)$, and for the corresponding wavenumber

where ${\nu }_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t},0}^2={\left[k_d\gamma /(4\eta D)\right]}^{1/2}=1/(2{\ell }_s{\ell }_D)$.

When $\kappa =0$ or β = 0, and hence δ = 0, we recover the predictions of an active gel model not accounting for network architecture (Hannezo et al., 2015), ${\lambda }_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}={\lambda }_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t},0}$ and ${\nu }_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}={\nu }_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t},0}$. The expression for ${\lambda }_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t},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 ($\kappa \ne 0$) and flow-induced alignment ($\beta <0$) couple the instability of Hannezo et al., 2015, to nematic order, changing the nature of pattern formation (see Equation 9, Equation 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 Damköhler lengths and on the strength and sign of nematic coupling. The nematic corrections increase as $c_0$, close to the point where the uniform quiescent state develops spontaneous order. We thus studied separately the regime $0<c_0\ll 1$, finding analogous expansions for the critical tension and wavenumber in terms of $\delta =\gamma D\kappa \beta /(2\eta L)$. Interestingly, Equation 9 shows that the activity threshold is reduced, and hence pattern formation is facilitated, when $\kappa <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 nonlinear finite element simulations in a periodic square 2D domain (Mirza et al., 2025; Mirza, 2025), with a domain size chosen to be $8{\ell }_{\mathrm{p}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}}$, where ${\ell }_{\mathrm{p}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}}=2\pi /{\nu }_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}$ is the pattern lengthscale predicted by linear stability analysis. In these simulations, we increased the activity parameter λ beyond the instability starting from a quiescent uniform state. We found that the linear stability analysis very accurately predicts the activity thresholds within 2% across a wide range of parameters. Furthermore, we found that the linear estimate for the wavenumber in Equation 10 characterizes well the nonlinear patterns, as quantified in Figure 2—figure supplement 1 and illustrated visually by the density patterns in Figure 2 showing the entire $8{\ell }_{\mathrm{p}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}}\times 8{\ell }_{\mathrm{p}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}}$ computational domain. In the nonlinear regime, the exponentially growing instabilities eventually reach out-of-equilibrium quasi-steady-state patterns involving self-reinforcing flows toward regularly spaced regions of high density surrounded by a low-density matrix, as previously reported for various compressible and density-dependent active gel models (Hannezo et al., 2015; Bois et al., 2011; Kumar et al., 2014; Callan-Jones and Voituriez, 2013; Mietke et al., 2019; Barberi and Kruse, 2023). In the absence of nematic coupling ($\beta =\kappa ={\lambda }_{\odot }=0$), these high-density domains are droplets arranged in a regular hexagonal lattice with S = 0 throughout the domain as previously reported (Kumar et al., 2014, Figure 2—video 1). In contrast, for a generic parameter set with finite β, κ, and ${\lambda }_{\odot }$, high-density domains adopt elongated configurations or bands where order is high, surrounded by a low-density and low-order matrix (Figure 2b, Figure 2—video 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 Zumdieck et al., 2005; Santhosh et al., 2020, 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.

Figure 2 with 3 supplements see all

Download asset Open asset

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 κ (Figure 2c). For $\kappa <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 (Figure 2bI). Instead, for κ > 0, nematic order is perpendicular to the axis of the bands, in what we call banded pattern (Figure 2bIII). To systematically study the effect of active tension anisotropy, we varied κ between –0.8 and 0.8 while keeping all other nondimensional groups fixed and setting λ to be 1.3 times the critical activity. We defined the order parameter ω (Figure 2a), allowing us to distinguish between banded ($\omega >0$) and fibrillar ($\omega <0$) organizations. We found a sharp transition between fibrillar and banded regimes around $\kappa \approx 0$, during which elongated high-density and high-order domains fragment into nematic droplets or tactoids (Weirich et al., 2017; Weirich et al., 2019; Figure 2bII).

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 (Lehtimäki et al., 2021). They also agree with observations showing that actin bundles form a mechanical continuum with the surrounding sparse and isotropic cortex (Vignaud et al., 2021). Their morphology and patterning dynamics is strikingly reminiscent of actin microridges, formed at the apical surfaces of mucosal epithelial cells (Depasquale, 2018; van Loon et al., 2020). 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 (Doostmohammadi et al., 2018), here, it is recycled by disassembly and turnover (Figure 2c, d, left, and Equation 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 (Huber et al., 2018; Denk and Frey, 2020) or mean-field models of dry mixtures of microtubules and motors (Maryshev et al., 2019).

At linear order, our theory shows that the distinctly nematic self-organization requires both flow-induced alignment (β) and active tension anisotropy ($\kappa$), whereas no condition is required on nematic activity (${\lambda }_{\odot }$). 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 $\kappa$, yet a finite value of β enhances fibrillar formation, leading to longer and more stable dense bands and hindering banded organization (Figure 3—video 1). This behavior is expected since the velocity gradients of the self-reinforcing flows tend to align filaments parallel to high-density bands due to $\beta {\mathit{d}}^{\mathrm{d}\mathrm{e}\mathrm{v}}$ in Equation 8.

The active nematic coefficient ${\lambda }_{\odot }$ modifies the onset of pattern formation through $c_0$ according to linear stability analysis. Apart from this linear effect, it should also contribute to the condensation of nematic order in high-density regions since it appears multiplied by $\rho$ in Equation 8. To examine this nonlinear effect, we perform simulations with ${\lambda }_{\odot }=0$ but keeping $c_0$ fixed. This leads to very different patterns without clear elongated structures and very mild nematic patterning (Figure 3a). Enhancing nematic patterning by considering the largest thermodynamically allowed value of $|\beta |$ 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 (Figure 3b), as observed in other active nematic systems (Zumdieck et al., 2005). This architecture is enhanced by higher friction γ (Figure 3c). Thus, the architectures found for ${\lambda }_{\odot }=0$ and κ < 0 are distinct from the fibrillar pattern described previously. Similarly, rather than bands with order perpendicular to their axis, for ${\lambda }_{\odot }=0$, $\kappa >0$, and high $|\beta |$, we found patterns of nematic asters (high-density droplets with radial nematic organization around them) (Figure 3b, c).

Figure 3 with 1 supplement see all

Download asset Open asset

Together, these results show that active tension anisotropy ($\kappa \ne 0$) and nematic activity (${\lambda }_{\odot }\ne 0$) are necessary and sufficient for nematic self-organization into fibrillar or banded patterns, with flow-induced alignment ($\beta <0$) favoring fibrillar organization. Furthermore, they show that while Equations 9 and 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 ${\lambda }_{\odot }$, in the nonlinear regime.

Given the morphological and dynamical diversity of nematic bundles in actin gels across cell types, geometric confinement, mechanical environment, or biological and pharmacological treatments (Dudin et al., 2019; Yolland et al., 2019; Jalal et al., 2019; Alert et al., 2022; Gupta et al., 2015; Xia et al., 2019), we varied model parameters to examine the architectures predicted by our active gel model, focusing on $\kappa <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 (Ennomani et al., 2016; Chugh, 2017).

With our default parameter set, bundle junctions and free ends are unfavorable and reorganize during pattern formation to annihilate as much as possible (Figure 4—video 1I). 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 Figure 4aI. We then asked the question of whether an orientational bias, which physically may be caused by cytoskeletal flows, boundaries, or directed polymerization (Hotulainen and Lappalainen, 2006), 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 (Figure 4aII, Figure 4—video 1II). Alternatively, we considered $c_0$ to be slightly negative, leading to the uniform and nematic steady state $\rho (x,t)={\rho }_0$, $v(x,t)=0$, and $q(x,t)=q_0=\pm \frac{1}{2}\sqrt{-c_0/b}$. The linear stability analysis around this state and further nonlinear simulations show that the essential phenomenology of Equation 9, Equation 10, and Figure 2 is not altered by the slight preexisting order (Appendix C). Again, the weak preexisting order provides sufficient bias to direct pattern orientation (Figure 4aIII). 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 (Dinwiddie et al., 2014), as a result of uniaxial cellular stretch (Wirshing and Cram, 2017; van Loon et al., 2020), or on anisotropically curved surfaces (Öztürk-Çolak et al., 2016; Hannezo et al., 2015).

Figure 4 with 4 supplements see all

Download asset Open asset

At longer times, the nematic bundles in Figure 4aIII develop secondary active instabilities leading to coordinated bending, curvature amplification, defect nucleation, and annihilation (Figure 4—video 1III) in a behavior reminiscent of active extensile systems (Aditi Simha and Ramaswamy, 2002; Ramaswamy and Rao, 2007; Doostmohammadi et al., 2018). This possibility is far from obvious because, even though the deviatoric part of active tension is extensile since $\kappa <0$, total active tension is contractile since $|\kappa |<1$ (Figure 1c). To systematically examine this point, we performed simulations at higher active tension anisotropies $\kappa =-0.8$, which we compared with our reference $\kappa =-0.2$ (Figure 4b). While in our reference system, bundles behave like contractile objects tending to straighten, for $\kappa =-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 Figure 4bII, in a behavior akin to active turbulence Verkhovsky et al., 1997). See Figure 4—video 2 for an illustration and for the slightly extensile case $\kappa =-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, ${\sigma }_{\parallel }-{\sigma }_{\perp }$. We note that for $\kappa <0$, as required to obtain fibrillar patterns, the active component of this quantity is negative, ${\sigma }_{\parallel }^{\mathrm{a}\mathrm{c}\mathrm{t}}-{\sigma }_{\perp }^{\mathrm{a}\mathrm{c}\mathrm{t}}<0$ (Figure 1c). In the case of bundles with contractile phenomenology, i.e., a tendency to straighten Figure 4bI, we found that ${\sigma }_{\parallel }-{\sigma }_{\perp }$ 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 Figure 4bII, we found intricate stress patterns with regions of negative ${\sigma }_{\parallel }-{\sigma }_{\perp }$. The extensile tension pattern for $\kappa =-0.8$ is more clearly appreciated in 1D simulations (Figure 4—figure supplement 1), 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 ($\kappa =-0.2$), even if ${\sigma }_{\parallel }^{\mathrm{a}\mathrm{c}\mathrm{t}}<{\sigma }_{\perp }^{\mathrm{a}\mathrm{c}\mathrm{t}}$, the negative viscous stresses due to the self-reinforcing flows are significantly larger in magnitude perpendicular to the nematic direction, resulting in ${\sigma }_{\parallel }>{\sigma }_{⟂}$. 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 (Balasubramaniam et al., 2022; Saw et al., 2017).

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, ${\ell }_{\mathrm{p}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}}=2\pi /{\nu }_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}$, becomes comparable to or smaller than the hydrodynamic lengthscale (Hannezo et al., 2015), thus enabling their mechanical interaction. In a model devoid of orientational order, reducing friction is equivalent to increasing ${\overline{k}}_d$. In our model, however, we can either reduce γ, which in nondimensional terms means increasing ${\bar{k}}_d$, a, b, and ${\overline{\lambda }}_{\odot }$ in concert, or increase ${\bar{k}}_d$ while leaving all other nondimensional 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 (Figure 4cI). 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 Figure 4—video 3I). We note that here ${\sigma }_{\parallel }>{\sigma }_{⟂}$ (Figure 4cI, Figure 4—figure supplement 1d), 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 ${\bar{k}}_d$ only, leads to very different networks with high-density aster-like clusters interconnected by straight actin bundles. Because now ${\bar{\lambda }}_{\odot }$ 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 ${\overline{k}}_d=10$, the network is stable and nearly crystalline (Figure 4cI), whereas for ${\overline{k}}_d=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) (Figure 4cIII, Figure 4—video 3III). This architecture and dynamics resemble those of early C. elegans embryos (Munro et al., 2004), adherent epithelial cells treated with epidermal growth factor (Jalal et al., 2019), and mouse embryonic stem cells (Gupta et al., 2015). Recent active gel models accounting for RhoA signaling develop similar pulsatile behaviors in 2D but do not predict the orientational order of the spatiotemporal patterns of the actomyosin cortex (Staddon et al., 2022).

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.

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 lengthscales, requires that active tension perpendicular to nematic orientation is larger than along this direction ($\kappa <0$, see Figure 1c), at least at the onset of pattern formation. This fact may seem counterintuitive in that dense nematic bundles are associated with large contractile tension along their axis, but as discussed in Figure 4cI, even if $\kappa <0$, bundles can be contractile because of the interplay between active and viscous stresses. Another prediction of our model is that the active assembly of dense nematic bundles requires active alignment controlled by parameter ${\lambda }_{\odot }>0$ (Figure 3). Being central to our conclusions, we then sought to examine the plausibility of the conditions $\kappa <0$ and ${\lambda }_{\odot }>0$ of our phenomenological theory by performing discrete network simulations using the open-source cytoskeleton simulation suite cytosim (Nedelec and Foethke, 2007).

Similarly to reconstituted actomyosin gels, discrete network simulations of the actin cytoskeleton tend to irreversibly collapse into clumps (Ennomani et al., 2016), although proper tuning of model parameters can lead to sustained contractile states (Mak et al., 2016). 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 $\kappa <0$ and ${\lambda }_{\odot }>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 (Figure 5a). Each of our simulated volume elements can represent a uniform system prior to pattern formation, or a material point in the actomyosin gel.

Figure 5 with 3 supplements see all

Download asset Open asset

We performed 2D simulations in which semiflexible filaments interact with crosslinkers and myosin motors, all of which undergo turnover and have a stoichiometry previously used to model the actin cytoskeleton (Cortes et al., 2020, Figure 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 the average nematic order to a target value S₀.

We first prepared a system consisting only of randomly oriented actin fibers, imposed the desired orientational order S₀ using the nematic penalty, and equilibrated the system. In a first set of simulations, once S₀ 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 boundary forces and hence estimate the effective active tension along the nematic direction ${\sigma }_{||}={\sigma }_{ij}{n}_i{n}_j$ and perpendicular to it, ${\sigma }_{\mathrm{\perp }}={\sigma }_{ij}{m}_i{m}_j$ with $n_i{m}_i=0$ and $m_i{m}_i=1$ (Figure 5b).

Addition of crosslinkers and myosins leads to bundling of actin filaments at the microscale (Figure 5—video 1) 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 ($S_0=0$), active tension is isotropic with ${\sigma }_{||}\approx {\sigma }_{\mathrm{\perp }}$. For anisotropic networks, however, we found that tension becomes anisotropic with ${\sigma }_{\mathrm{\perp }}>{\sigma }_{||}$ (Figure 5c).

We systematically characterized this behavior by varying initial orientational order and network density. According to our active gel model, Equation 5, Equation 6, Equation 7, in the absence of nematic gradients and flow, the tension components ${\sigma }_{||}$ and ${\sigma }_{\perp }$ satisfy the following relations in terms of mean and deviatoric tensions

Remarkably, our discrete network simulations closely followed these relations (Figure 5d, e), which allowed us to estimate $\kappa \approx -1.6$. We robustly found that $\kappa <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, $\rho {\lambda }_{\odot }$. For a uniform system with nematic order along a given direction and ignoring the susceptibility parameter b, Equation 8 becomes

where ${\eta }_{\mathrm{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, $\rho {\lambda }_{\odot }$ 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 ${\mathcal{K}}_S$. 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 ${\lambda }_{\odot }=0$ and a = 0, Equation 12 predicts an exponential relaxation given by $S(t)=S_0(1-e^{-{\mathcal{K}}_{S}t/{\eta }_{\mathrm{r}\mathrm{o}\mathrm{t}}})$, which very closely matched the simulation data for different values of S₀ and for a single fitting parameter ${\eta }_{\mathrm{rot}}$ (Figure 5—figure supplement 1III, Figure 5—video 2). 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 S₀. In agreement with the notion of an active force driving nematic order, we found that $S(t)$ monotonically increased (Figure 5f). More quantitatively, we tested the short-time prediction of Equation 12, ${\eta }_{\mathrm{r}\mathrm{o}\mathrm{t}}\dot{S}=(\rho {\lambda }_{\odot }-2a)S_0$, by plotting $\dot{S}$ as estimated from our simulations, as a function of S₀ (Figure 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 ($\kappa <0$), and (2) active generalized forces driving further ordering.

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 $\rho (x,t)$ along the x-axis. The nematic order tensor is defined as $q_{ij}=S(n_i{n}_j-{\delta }_{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 two independent degrees of freedom, $q_{11}=-q_{22}=q_1$ and $q_{12}=q_{21}=q_2$. In the 1D model, we assume that n is either along or perpendicular to x, and hence $q_{12}=q_{21}=0$. We are thus left with one independent degree of freedom $q_{11}=-q_{22}=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 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 $c_0=2a-{\rho }_0{\lambda }_{\odot }\ge 0$, the uniform steady state of the system is given by $\rho (x,t)={\rho }_0$, $v(x,t)=0$, and $q(x,t)=0$. For $c_0<0$, the uniform steady state has spontaneous alignment given by $q(x,t)=q_0=\pm \frac{1}{2}\sqrt{-c_0/b}$. Finally, we note that, for the dissipation to be positive, the material parameters need to satisfy an entropy production inequality $4\eta {\eta }_{\mathrm{r}\mathrm{o}\mathrm{t}}-{\beta }^2\ge 0$ in the reduced 1D model and $2\eta {\eta }_{\mathrm{r}\mathrm{o}\mathrm{t}}-{\beta }^2\ge 0$ in the 2D model (Mirza et al., 2025).

We present next a nondimensionalization of the governing Equations A2–A5 . We choose the hydrodynamical lengthscale ${\ell }_s=\sqrt{\eta /\gamma }$, above (below) which friction (viscosity) is the dominant dissipative mechanism in the active gel, and the timescale of diffusion over this lengthscale, $t_0={\ell }_s^2/D=\eta /(\gamma D)$. Using these characteristic scales, we nondimensionalize space and time as $\overline{x}=x/{\ell }_s$ and $\overline{t}=t/t_0$. We chose ${\rho }_0$ as a reference density and thus $\overline{\rho }=\rho /{\rho }_0$. Here and elsewhere, overbar denotes nondimensional quantities.

The dimensionless balance of mass reads

where ${\bar{k}}_d=k_d{t}_0$. The dimensionless balance of linear momentum reads

where $\overline{\sigma }=\sigma /{\sigma }_0$ and the characteristic stress is given by ${\sigma }_0={\rho }_0\gamma D$. The nondimensional constitutive relation for the stress is

where $\overline{L}=L/(\gamma D{\ell }_s^2)$, $\overline{\beta }=\beta /\eta$, and $\overline{\lambda }=\lambda /(\gamma D)$. We note that $\kappa ={\lambda }_{\mathrm{a}\mathrm{n}\mathrm{i}\mathrm{s}\mathrm{o}}/\lambda$ and q are already nondimensional. Finally, the dimensionless equation of generalized force conjugate to nematic order is given by

where ${\overline{\eta }}_{\text{rot}}={\eta }_{\text{rot}}/\eta$, $\overline{a}=a/\left(\gamma D\right)$, $\overline{b}=b/\left(\gamma D\right)$, and ${\overline{\lambda }}_{\odot }={\lambda }_{\odot }{\rho }_0/(\gamma D)$. Inspection of Equation B4 reveals the characteristic nematic lengthscale ${\ell }_q^2=L/\left|2a-{\lambda }_{\odot }{\rho }_0\right|$ and the nematic relaxation timescale $t_q=({\ell }_q^2{\eta }_{\mathrm{r}\mathrm{o}\mathrm{t}})/L={\eta }_{\mathrm{r}\mathrm{o}\mathrm{t}}/\left|2a-{\lambda }_{\odot }\rho \right|$.

We perform next the linear stability analysis on the 1D model given by the dimensionless Equations B1–B4. 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)=q_0$, and $\rho (x,t)={\rho }_0$ by examining the fate of infinitesimal perturbations $v\left(x,t\right)=\delta v\left(x,t\right)$, $q\left(x,t\right)=q_0+\delta q\left(x,t\right)$, and $\rho \left(x,t\right)={\rho }_0+\delta \rho \left(x,t\right)$, where $\left|\delta v\right|\sim \left|\delta q\right|\sim \left|\delta \rho \right|\ll 1$. The control parameter in our stability analysis is the activity parameter λ.

Inserting the infinitesimally perturbed homogeneous states into the governing Equations B1–B4 and ignoring terms that are quadratic in perturbations, the resulting balance of linear momentum equation is

The linearized form of Equation B4 is

If $c_0=2a-{\rho }_0{\lambda }_{\odot }>0$, then $q_0=0$ and $G=c_0$. If $c_0=2a-{\rho }_0{\lambda }_{\odot }<0$, then $q_0^2=-c_0/(4b)$ and $G=8b{q}_0^2$. The linearized balance of mass reads

We decompose perturbations into a sum of Fourier terms, i.e., $\delta v=v^{\prime }{e}^{\alpha t+\mathrm{i}x\nu }$, $\delta \rho ={\rho }^{\prime }{e}^{\alpha t+\mathrm{i}x\nu }$, and $\delta q=q^{\prime }{e}^{\alpha t+\mathrm{i}x\nu }$, where $v^{\prime }$, $q^{\prime }$, and ${\rho }^{\prime }$ are the amplitudes of the perturbations, α is the growth rate, and $\nu$ is the wavenumber.

Substituting the Fourier representation in Equation C1, we obtain

Likewise, Equation C3 becomes

Substituting this equation into the Fourier representation of Equation C2, we obtain

Finally, combining the three equations above, we obtain the dispersion equation, which, given the material parameters, relates wavenumber $\nu$ to growth rate α and takes the form

where

The roots of Equation C7 are given by $\alpha (\nu )=\left[-B(\nu )\pm \sqrt{B(\nu {)}^2-4A(\nu )C(\nu )}\right]/(2A(\nu ))$. The system is linearly stable if perturbations decay, i.e., $\text{Real}(\alpha )<0$, marginally stable if $\text{Real}(\alpha )=0$, and unstable if $\text{Real}(\alpha )>0$. We found that for all reasonable material parameters $\alpha (\nu )$ is real, and hence the critical condition is given by $\alpha =0$, which recalling Equation C7 implies that $C(v)=0$. From this condition, we can express the activity parameter as

To find the critical wavenumber ${\nu }_{\mathrm{crit}}$ and the associated critical activity parameter ${\lambda }_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}=\lambda ({\nu }_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}})$, we minimize this expression with respect to $\nu$, or equivalently with respect to $x={\nu }^2$, which requires that

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 $\eta ={10}^4$ Pa s, and hydrodynamic lengths in the order of ${\ell }_s=0.2\mathrm{i}\mathrm{c}\mathrm{r}\mathrm{o}\mathrm{m}$ (Hannezo et al., 2015; Salbreux et al., 2009), although we vary this parameter, as indicated in Appendix 4—Tables 1 and 2. Taking $k_d\sim 0.1{\mathrm{s}}^{-1}$ (Yolland et al., 2019; Hannezo et al., 2015; Callan-Jones et al., 2008) and $D\sim 0.04$ μm² s⁻¹ (Yolland et al., 2019; Hannezo et al., 2015), we obtain Damkölher lengths ${\ell }_D$ in the order of a micron. We consider ${\eta }_{\mathrm{r}\mathrm{o}\mathrm{t}}/\eta \sim 1$ and $\beta /\eta \sim -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 ${\lambda }_{\mathrm{crit}}$. We vary the nondimensional active tension anisotropy parameter $\kappa$. We choose the susceptibility parameters as $a/{\eta }_{\mathrm{r}\mathrm{o}\mathrm{t}}=5$ s⁻¹ and $b/{\eta }_{\mathrm{r}\mathrm{o}\mathrm{t}}=20$ s⁻¹ so that relaxation of nematic order is significantly faster than the turnover timescale and choose the Frank constant to be small so that ${\ell }_q$ 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 Appendix 4—table 1 and Appendix 4—table 2.

To ascertain the microstructural origin of the anisotropic active tensions underlying the formation of stress fibers and provide evidence for the orientational activity parameter ($\rho {\lambda }_{\odot }$) included in our theory, we establish an agent-based microscopic model of a crosslinked actomyosin network using the open-source cytoskeletal simulation suite cytosim (Nedelec and Foethke, 2007; Nedelec, 2022). We customized the source code to impose and track nematic order in the system (Pensalfini, 2025). Below, we provide a concise description of the modeling approach and an overview of the simulation parameters.