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 ρ(t,x) and by the network architecture given by the symmetric and traceless nematic tensor q(t,x) (see Figure 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, S=2qijqij 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 d=12(∇v+∇vT) measures the local rate of distortion of the fluid, whereas w=12(∇v−∇vT) measures its local spin, where ∇ is the nabla differential operator. The deviatoric part of the rate-of-deformation tensor is defined by dijdev=dij−(dkk/2)δij. The rate of change of q relative to a frame that translates and locally rotates with the flow generated by v is given by the Jaumann derivative q^=∂q/∂t+v⋅∇q−w⋅q+q⋅w (Gennes and Prost, 1993).

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 β.

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)

(1)

∂ρ∂t+∇⋅(ρv)−DΔρ+kd(ρ−ρ0)=0,

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.

Force balance in the gel can be expressed as

(2)

ργv=∇⋅σ,

where γ>0 models a viscous drag with the surroundings (e.g. the plasma membrane) and σ=σnem+σ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=∫Ωρf(q,∇q)dS, where

(3)

f(q,∇q)=12aS2+18bS4+12L|∇q|2

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)

(4)

hij=−δFδqij=−ρ(2a+bS2)qij+L∇k(ρ∇kqij),

and to the explicit form of the nematic nonsymmetric stress tensor

(5)

σijnem=−ρ∂f∂∇jqlk∇iqlk+qikhjk−qjkhik=L[−ρ∇iqkl∇jqkl+qik∇l(ρ∇lqjk)−qjk∇l(ρ∇lqik)].

The dissipative part of the stress

(6)

σijdiss=ρ[2η(dij+dkkδij)+βq^ij],

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

(7)

σijact=ρ(λδij+λanisoqij)=ρλ(δij+κqij).

The 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 (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 q^ also includes viscous, elastic-nematic, and active contributions and takes the form

(8)

ηrotq^+βddev−1ρh−ρλ⨀q=0.

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 (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ηηrot−β2≥0 (Mirza et al., 2025). 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 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 λ⊙≥0 tending to further align filaments (Figure 1d; Reymann et al., 2016). This term is linear in ρ because in the expansion λ¯⊙+ρλ⊙, the constant contribution λ¯⊙ can be subsumed by the susceptibility parameter a (Salbreux et al., 2009). 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, v(x,t)=0 and a nonzero nematic order parameter S02=−c0/b. 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 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, λ, κ 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 (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 β, λ⊙, κ and L.

We readily identify the hydrodynamic length ℓs=η/γ, above which friction dominates over viscosity, the Damköhler length ℓD=D/kd above which reactions dominate over diffusion, and the nematic length ℓq=L/|2a−λ⊙ρ0|. Nondimensional analysis reveals a set of nondimensional groups that control the system behavior, namely the nondimensional turnover rate k¯d=ℓs2/ℓD2, the Frank constant L¯=L/(ηD), the susceptibility parameters a¯=a/(γD) and b¯=b/(γD), the drag coefficients η¯rot=ηrot/η and β¯=β/η, the nematic activity coefficient λ¯⊙=ρ0λ⊙/(γD), and the active tension parameters λ¯=λ/(γ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 .