Active nematics with deformable particles

Abstract

The hydrodynamic theory of active nematics has been often used to describe the spatio-temporal dynamics of cell flows and motile topological defects within soft confluent tissues. Those theories, however, often rely on the assumption that tissues consist of cells with a fixed, anisotropic shape and do not resolve dynamical cell shape changes due to flow gradients. In this paper we extend the continuum theory of active nematics to include cell shape deformability. We find that circular cells in tissues must generate sufficient active stress to overcome an elastic barrier to deforming their shape in order to drive tissue-scale flows. Above this threshold the systems enter a dynamical steady-state with regions of elongated cells and strong flows coexisting with quiescent regions of isotropic cells.

---

1. Introduction

Living systems are inherently active as they use chemical energy from their surroundings to do work. It is becoming increasingly clear that many features seen in cell layers and tissues can be understood using the continuum theory of active nematics, which describes hydrodynamic interactions between active anisotropic particles. Examples include biofilm initiation, topological defects in cell monolayers and epithelial expansion.^1–9

While the theory of active nematics predicts many qualitative features such as short-range orientational order, active turbulence and motile topological defects,^10–14 it is based on liquid crystal hydrodynamics which assumes that the active particles are nematogens with a fixed aspect ratio.¹⁵ While this assumption may be a reasonable description for systems consisting of rod-shaped cells such as fibroblasts or Escherichia coli,¹⁶ many particles, such as MDCK cells or soft colloids, can undergo large shape changes and their aspect ratio can vary significantly due to active forces.^17–19 Simulations have shown that intercellular stresses in monolayers are enhanced by cell deformation, creating a positive feedback loop that affects the collective behaviour of the layer.²⁰

In this paper we extend the continuum equations which describe active nematic fluids by considering not only the magnitude and direction of the nematic order, but also the shape dynamics of the underlying particles. We consider an equation of motion for the aspect ratio of the deformable nematogens which incorporates flow-driven stretching and compression of particles, as well as elastic restoring forces. The aspect ratio of particles in turn affects the orientational dynamics by modifying thermodynamic interactions between particles and the way they align in shear flows.

In Section 2 we present the equations of motion for particle shapes, characterised by the aspect ratio ω, the nematic order tensor Q of the particle orientational distribution, and the associated velocity field u. The following Sections, 3.1 and 3.2, are analytical and numerical investigations of systems consisting of elastic, active particles which are isotropic in the absence of flow. We report a shape instability in which extensile active stress generated by particles drives the formation of regions with highly anisotropic particles with nematic order. Finally, in Section 4 we summarize our results.

2. Equations of motion

We introduce coarse-grained equations of motion to describe the hydrodynamics of elastic, nematic particles. To this end we approximate the shape of particles as ellipsoidal, and define a field ω > 1 to describe their local aspect ratio, Fig. 1(a). The time evolution of the shape distribution follows^21–23The first term in eqn (1) describes advection by a macroscopic flow field u, while the second term accounts for shape changes due to flow gradients. E_‖ = n_iE_ijn_j is the projection of the strain rate tensor E_ij = (∂_iu_j + ∂_ju_i)/2 along the long axis of the particles defined by the unit vector n. Depending on the particle orientation with respect to the extensional flow axis, particles are either stretched or compressed by flows (Fig. 1(b) and (c)). The final term in eqn (1) models particle shape elasticity, where Γ_ω controls the relaxation towards the minimum of a free energy [scr F, script letter F] as defined below.

Fig. 1 Particle deformations driven by (a) free energy minimisation, (b) flow stretching, (c) flow compression, (d) flow alignment.

Nematic ordering of anisotropic particles is described by a tensorial order parameter, Q_ij = 3S/2(n_in_j − δ_ij/3), where S quantifies the magnitude of the nematic alignment, and we restrict the directors to lie within a plane so that n_z = 0 and u_z = 0. The time evolution of Q follows the Beris–Edwards equation,²⁴

(∂_t + u·∇)Q − [scr W, script letter W] = Γ_LCH, (2)

where Γ_LC controls the relaxation towards an equilibrium quantified by the molecular field H = −δ[scr F, script letter F]/δQ + (I/3)Tr(δ[scr F, script letter F]/δQ). Nematogens are not only advected by the fluid, but also rotated by gradients in the flow field, which gives rise to the the co-rotational term²⁴where Ω_ij = (∂_ju_i − ∂_iu_j)/2 is the vorticity tensor, the antisymmetric part of the velocity gradient tensor W_ij = ∂_iu_j. The flow-alignment parameter ξ quantifies the rotation of anisotropic particles in shear flows (Fig. 1(d)). Its value depends on the shape of the particles and it is zero for isotropic particles.²⁵ Therefore we assume ξ = ξ₀(ω − 1) where ξ₀ sets the flow-alignment scale.

The total free energy of the system, , consists of two contributions. The first contribution is an elastic energy associated with particle shape deformations,

where ω₀ is the particle aspect ratio at equilibrium and A_ω, are elastic deformation parameters. The second represents the liquid crystal free energy which emerges from particle–particle interactions,where K_LC penalizes distortions in the director field in the one-elastic-constant approximation and A_LC sets the bulk-energy scale. Highly anisotropic particles tend to align nematically due to excluded volume effects whereas there is no alignment interaction between circular particles. This is reflected in the free energy which has a minimum at a reference value, S = S₀, for ω → ∞ while favouring S = 0 for ω = 1.

Finally, the velocity field u of the fluid obeys the incompressible Navier–Stokes equations,

∇·u = 0, (6)

ρ(∂_t + u·∇)u = ∇·Π, (7)

where ρ is the fluid density and the stress tensor Π = Π^passive + Π^active includes passive and active contributions. The passive stress tensor consists of viscous dissipation and elastic stress arising from liquid crystal hydrodynamics,²⁶

Π^viscous = 2ηE, (8)

where η is the viscosity, p is the bulk pressure and . Active dipolar forces produced by individual particles give rise to an active stress,²⁷

Π^active = −ζQ, (10)

where ζ quantifies the magnitude of active forces. ζ > 0, corresponds to extensile activity, where fluid is pushed outwards from the particles along their direction of elongation, and pulled inwards along the perpendicular axis. ζ < 0 corresponds to contractile activity, where the flow direction is reversed. We consider extensile systems in what follows.

We use a hybrid lattice Boltzmann–finite difference method to solve the equations of motion (1), (2), (6), (7).²⁸ The lattice size is 200 × 200 with a lattice spacing Δx = 1, LB timestep Δt = 1 and periodic boundary conditions. Unless stated otherwise, the simulation parameters are Γ_LC = 0.1, A_LC = 0.1, S₀ = 1, K_LC = 0.015, ξ₀ = 0, Γ_ω = 0.01, A_ω = 0.04, = 0.003, ω₀ = 1, ζ = 0.001, ρ = 1, η = 2/3 The relevant inverse time-scales are as follows: Γ_ωA_ω, the rate of relaxation of the particles to their preferred shape, Γ_LCA_LC, the rate of relaxation of the order parameter towards the free energy minimum, and ζ/η, the rate of active stresses injected into the system. In the absence of flow gradients, the regime where Γ_ωA_ω ≫ Γ_LCA_LC has nematic alignment responding slowly to changes in shape, creating a time delay associated with how particles will align after their shape changes. In the inverse case, nematic alignment responds quickly to shape changes, which is the regime we are considering. The relevant length-scales are the active length-scale and a passive length-scale , which sets the scale of shape gradients. Simulations are initialised with random noise in ω around a finite value of ω = 1.3 and S = (ω − 1)/ω. We use a runtime of 2 × 10⁵ LB timesteps. Initial conditions become irrelevant after a runtime of ≈10⁵ timesteps. If ω falls below unity, we rotate the director field to point towards the extensile flow axis. This is physically motivated by considering a particle that is instantaneously circular; any subsequent elongation will occur along the local extensile axis of the active flow field.

3. Results

To illustrate the coupled time evolution of particle orientation n and aspect ratio ω we consider a passive system in an external extensional flow and solve eqn (2) for the director field in the Ericksen–Leslie limit (constant S)²⁴ for a finite flow-aligning parameter which scales linearly with the aspect ratio, coupled to eqn (1) with Γ_ω = 0. Depending on the position in phase space, there are two distinct re-orientation mechanisms: If the angle θ between the compressive flow axis and particle orientation is sufficiently large, particles undergo rotation where their shape stays anisotropic and their orientation rotates until it aligns with the extensional flow axis (red trajectory in Fig. 2). If θ is small, however, elongated particles cannot rotate fast enough and get squeezed by the flow until they become isotropic, ω = 1. Subsequently their orientation aligns along the extensional flow axis, Q_ij ∼ E_ij, and their aspect ratio increases again (blue trajectory in Fig. 2).

Fig. 2 Left: phase space trajectories of the evolution of ω subject to an extensional flow at an angle θ to the director. Right: rotation and squeezing of a particle in an extensional flow.

3.1 Linear stability analysis

To investigate the behaviour of an active system consisting of nearly circular particles, we perform a stability analysis around a quiescent isotropic system by adding small perturbations to the aspect ratio of particles, ω = 1 + δω. The shape perturbation δω ∼ exp(iq·r + λt) is applied along a wavevector q = q(cos θ, sin θ, 0) and we assume particles deform along their long axis. The flow fields which result from activity give rise to the dispersion relation,

For a solution with λ > 0 to exist for some value of q, we require the coefficient of the quadratic term to be negative. Thus the critical value of the activity below which an isotropic system will be stable is given by:

ζ_c is clearly minimised when θ = π/4, (Fig. 3). At this angle, the parallel and perpendicular components of the flow work in tandem to create the most efficient elongation of the particles (ESI†). Previously, activity was found to create nematic order due to flow aligning effects as long as ξζ > 0.²⁹ Here we provide a mechanism where activity may induce nematic order without the need for flow alignment. Starting from an isotropic system of circular particles, activity causes elongation and then thermodynamic forces described by the free energy may allow for nematic order to be established.

Fig. 3 Dispersion relation of λ(q) from eqn (11) plotted for θ = π/4 for different values of {ζ, K_LC }. Red/Blue curves correspond to activity values smaller/greater than the threshold activity, ζ_c respectively.

The competition becomes apparent if we explicitly write down the wavevector below which the instability grows:

This depends on a balance between the active length-scale, , and the passive length-scale, .

3.2 Active turbulence

For values of the extensile activity above which flows are maintained the system enters an active turbulent regime. Fig. 4 shows the dynamical steady state corresponding to fully developed active turbulence in the system of deformable particles (Movie, ESI†). At any given point in time the system is characterised by distinct regions of elongated particles separated by areas where there is little or no particle extension. There is strong nematic ordering in the regions where the particles are extended, and it is possible to identify the motile topological defects characteristic of active turbulence. In active turbulence with particles of fixed length, nematic domains are broken up by the active instability and by the passage of topological defects. Here there is an additional mechanism which results in regions of circular particles. Within a nematic domain the gradients in Q are small, and thus active forces are diminished. Consequently, the magnitude of E_‖ drops, and it is no longer able to overcome the free energy cost of deformation, eqn (4). The particles within these domains will therefore begin to contract. Fig. 5 shows an idealised sketch of this process.

Fig. 4 Fully developed active turbulence with deformable particles. (a) Aspect ratio, ω of particles. (b) Velocity field, u. Colour scale indicates the magnitude of the velocity. (c) Nematic director field n and scalar order parameter S. Colour scale indicates the magnitude of S. Red/blue symbols correspond to the core of ±1/2 defects respectively. +1/2 defects are oriented such that the tail is along the line shown.

Fig. 5 Creation and contraction of nematic domains. Nematic elasticity smooths out gradients in the director field, resulting in a domain of aligned directors. The active force vanishes and particles contract due to the elastic free energy cost.

To investigate how physical parameters affect the dynamical state of the system of deformable particles, we measure the time and space averaged aspect ratio, 〈ω〉 as the activity parameter ζ, the elastic deformation parameter A_ω, the elastic constant, K_LC and the flow alignment parameter, ξ₀ are varied. In all cases, activity was high enough to keep the system in active turbulence. Fig. 6(a) shows that 〈ω〉 increases linearly with the activity for extensile systems. In active turbulence where the evolution of the flow field is dominated by driving from the active stress, the magnitude of the velocity, and consequently the strain rate, increase with the activity parameter, resulting in a higher degree of elongation. Conversely, Fig. 6(b) shows that 〈ω〉 decreases linearly with A_ω. As the elastic energy cost of deforming particles increases, the average aspect ratio will decrease. Next, as K_LC is increased, gradients in the director field are smoothed out, reducing the overall active stress and consequently reducing elongation. This effect is relatively small compared to the other parameters, as Fig. 6(c) shows. Finally, 〈ω〉 also increases linearly with the value of the flow aligning parameter, ξ₀, shown in Fig. 6(d). As ξ₀ increases, particles will align more efficiently with the extensional flow axis, thus increasing their elongation.

Fig. 6 Average aspect ratio as a function of (a) the activity, ζ, (b) elastic deformation parameter, A_ω, (c) the elastic constant, K_LC and (d) the flow alignment parameter ξ₀. A linear relationship using a least-squares fit has been plotted in each case. Error bars are given by the standard deviation of the spatial average of the aspect ratio over 10⁵ timesteps.

Next, we look at any emergent length-scales in the system. As we have already shown, in a system of fully developed active turbulence, distinct domains of elongated, anisotropic particles coexist with domains of isotropic particles. To quantify the size of these domains, we calculate the normalised correlation function defined as C_ω(r) = 〈(ω(0,t) − [small omega, Greek, macron](t))(ω(r,t) − [small omega, Greek, macron](t))〉/〈(ω(0,t) − [small omega, Greek, macron](t))²〉 where [small omega, Greek, macron] is the spatial average of the aspect ratio at each point in time and 〈·〉 denotes a space and time average over all points separated by a distance r. Fig. 7 shows how C_ω varies with distance, scaled by the active length-scale, for different choices of the activity and elasticity parameters ζ and K_LC. The collapse of the data to a single curve shows that the size of the domains is governed by L_act, as in active turbulence with fixed length nematogens, reflecting that flow gradients, which occur on this length scale are responsible for the correlations in the aspect ratio of particles.

Fig. 7 Aspect ratio correlation, C_ω as a function of distance, scaled by the active length-scale, . Different coloured curves correspond to differing values of L_act, generated by varying ζ and K_LC independently in systems exhibiting fully developed active turbulence.

We also measure the density of topological defects throughout the system as the activity ζ, the elastic constant K_LC and the cell elasticity A_ω are varied. To obtain this data we track defects^30,31 and also set a cut-off value ω = 1.1 below which topological defects cannot be observed.

As Fig. 8 shows, the defect density scales as where the modified active length-scale is defined as , where ζ* is the threshold activity required to keep the system in fully developed active turbulence – this is taken to be the activity value below which circular cells dominate the system and consequently, topological defects cannot exist. An estimate for this threshold activity is given below.

Fig. 8 Average number of topological defects as a function of the inverse modified active length-scale, multiplied by L_nem. The plotted line shows a least-squares fit through all the data points. Data points were collected by varying the activity parameter, ζ for different realisations of the parameters {K_LC, A_ω}. Error bars are given by the standard deviation of the spatial average of the defect number over 10⁵ timesteps.

This is an empirical result and is a fundamentally different scaling than fixed length active nematics where the defect density scales with the inverse active length-scale squared, L_act⁻² = ζ/K_LC.³² Indeed it seems very reasonable that other length-scales should be involved. One clear difference from the usual active nematic models is that the active stresses depend on the local magnitude of the order parameter and therefore vary across the system, particularly at interfaces, and this will influence the defect count. A plausible explanation for this scaling is as follows: plays a similar role to L_act in traditional active nematics, with a constant offset in the activity due to a minimum threshold required to achieve active turbulence. Secondly, the deformation energy A_ω will tend to reduce the magnitude of the aspect ratio, ω, and therefore the order parameter, S. In turn, this will reduce the effectiveness of the nematic elasticity which depends on S. There is then a reduced, effective elastic constant ∼K_LC/A_ω so the defect density scales with .

An estimate for the threshold activity ζ* may be obtained by using the velocity field around +1/2 defects³³ to calculate the strain rate, E⁺_‖. For a director given by n = (cos θ/2, sin θ/2,0) and assuming S = (ω − 1)/ω outside the core, the strain rate is given in polar coordinates by

By substituting (14) into eqn (1), neglecting advection and quartic contributions in the free energy, a simple calculation yields a lower bound to the activity required for particles to elongate along the axis of symmetry of +1/2 defects:

ζ* = 2Γ_ωA_ωη. (15)

This estimate ignores elasticity: it is lower than the numerical values by between 10% for low values of the elastic constant and 50% for high values of the elastic constant.

4. Discussion

We have extended the continuum theory of active nematics to allow for deformable particles. This was achieved by introducing an equation of motion (1) for the aspect ratio ω of particles, which couples to the equations for the nematic order parameter (2) and velocity fields (7) through the flow alignment parameter ξ and the nematic bulk free energy f_LC.

By analysing the linear stability of isotropic systems with initially circular nematogens, we showed that they must overcome an activity threshold in order to extend sufficiently to drive flows. Above this threshold the system enters a dynamical steady state characterised by coexisting regions of elongated particles, which tend to align nematically and are associated with topological defects and spatiotemporal chaotic flows, and quiescent regions consisting of isotropic particles. While isotropic particles get elongated by tissue-scale flow gradients from neighbouring nematic regions, there is little extensional flow deep in the nematic regions, where the aspect ratio of particles subsequently decreases due to elastic shape energy. In steady-state, particles continuously get stretched by flow gradients and subsequently contract in the absence of flow gradients, leading to large variations of cell shape. Similar cell shape variations have been observed in simulations of the vertex model for unjammed tissues³⁴ and interpreted in terms of a continuum model.³⁵ Mean field theories of cell shape variations based on the vertex cellular Potts models have also been presented.^36,37

We have focused on extensile materials. Contractile systems which are initially circular will not show any dynamical behaviour because active forces just reinforce the thermodynamic forces restoring the particles to circular. We note, however, that our model is restricted to nematic driving: polar forces which can lead to coherent motion of circular cells are also observed in cell monolayers, and the resultant force on a cell does not necessarily act along its nematic axis.^38,39

Our model was motivated by the observation of active turbulence in confluent cell layers where cells can be deformed. In the light of our results it would be interesting to further investigate the spatial and temporal variation of the aspect ratio of cells within confluent layers that are on average circular, such as the MDCK cell line.^8,20

Author contributions

All authors worked on defining the project, analysing data and writing the paper. I. H. ran the simulations.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This project was funded by the Gould & Watson Scholarship and the European Commission's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 812780.

Notes and references

1. C. Blanch-Mercader and J. Casademunt, Soft Matter, 2017, 13, 6913–6928.
2. C. Blanch-Mercader, V. Yashunsky, S. Garcia, G. Duclos, L. Giomi and P. Silberzan, Phys. Rev. Lett., 2018, 120, 208101.
3. T. B. Saw, W. Xi, B. Ladoux and C. T. Lim, Adv. Matter., 2018, 30, 1802579.
4. Y. I. Yaman, E. Demir, R. Vetter and A. Kocabas, Nat. Commun., 2019, 10, 1–9.
5. P. Pearce, B. Song, D. J. Skinner, R. Mok, R. Hartmann, P. K. Singh, H. Jeckel, J. S. Oishi, K. Drescher and J. Dunkel, Phys. Rev. Lett., 2019, 123, 258101.
6. G. Duclos, C. Erlenkämper, J.-F. Joanny and P. Silberzan, Nat. Phys., 2017, 13, 58–62.
7. D. Pearce, J. Nambisan, P. Ellis, A. Fernandez-Nieves and L. Giomi, Phys. Rev. Lett., 2021, 127, 197801.
8. T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans and B. Ladoux, Nature, 2017, 544, 212–216.
9. K. Copenhagen, R. Alert, N. S. Wingreen and J. W. Shaevitz, Nat. Phys., 2021, 17, 211–215.
10. A. Ardaševa and A. Doostmohammadi, Nat. Rev. Phys., 2022, 4, 354–356.
11. R. Mueller, J. M. Yeomans and A. Doostmohammadi, Phys. Rev. Lett., 2019, 122, 048004.
12. S. P. Thampi and J. M. Yeomans, Eur. Phys. J. Special Topics, 2016, 225, 651–662.
13. R. Alert, J. Casademunt and J. Joanny, Annu. Rev. Condens. Matter Phys., 2022, 13, 143–170.
14. A. Doostmohammadi, S. P. Thampi and J. M. Yeomans, Phys. Rev. Lett., 2016, 117, 048102.
15. P.-G. De Gennes and J. Prost, The Physics of Liquid Crystals, Oxford University Press, 1993.
16. D. Dell’Arciprete, M. Blow, A. Brown, F. Farrell, J. S. Lintuvuori, A. McVey, D. Marenduzzo and W. C. Poon, Nat. Commun., 2018, 9, 1–9.
17. D. Vlassopoulos and M. Cloitre, Curr. Opin. Colloid Interface Sci., 2014, 19, 561–574.
18. T. Lecuit and P.-F. Lenne, Nat. Rev. Molecular Cell Bio., 2007, 8, 633–644.
19. C.-P. Heisenberg and Y. Bellache, Cell, 2013, 153, 948–962.
20. L. Balasubramaniam, A. Doostmohammadi, T. B. Saw, G. H. N. S. Narayana, R. Mueller, T. Dang, M. Thomas, S. Gupta, S. Sonam and A. S. Yap, et al. , Nat. Mater., 2021, 20, 1156–1166.
21. B. A. Bilby and M. Kolbuszewski, Proc. Roy. Soc. London. A, 1977, 355, 335–353.
22. T. Gao, H. H. Hu and P. P. Castañeda, J. Fluid Mech., 2011, 687, 209–237.
23. T. Gao, H. H. Hu and P. P. Castañeda, J. Fluid Mech., 2013, 715, 573–596.
24. A. N. Beris and B. J. Edwards, Thermodynamics of flowing systems with internal microstructure, Oxford University Press, New York; Oxford, 1994.
25. G. B. Jeffery, Proc. Roy. Soc. London. A, 1922, 102, 161–179.
26. D. Marenduzzo, E. Orlandini and J. M. Yeomans, Phys. Rev. Lett., 2007, 98, 118102.
27. R. A. Simha and S. Ramaswamy, Phys. Rev. Lett., 2002, 89, 058101.
28. D. Marenduzzo, E. Orlandini, M. Cates and J. M. Yeomans, Phys. Rev. E., 2007, 76, 031921.
29. S. Santhosh, M. R. Nejad, A. Doostmohammadi, J. M. Yeomans and S. P. Thampi, J. Stat. Phys., 2020, 180, 699–709.
30. M. L. Blow, S. P. Thampi and J. M. Yeomans, Phys. Rev. Lett., 2014, 113, 248303.
31. A. J. Vromans and L. Giomi, Soft Matter, 2016, 12, 6490–6495.
32. L. Giomi, Phys. Rev. X, 2015, 5, 031003.
33. L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek and M. Cristina Marchetti, Philos. Trans. Roy. Soc., A, 2014, 372, 20130365.
34. D. Bi, X. Yang, M. C. Marchetti and M. L. Manning, Phys. Rev. X, 2016, 6, 021011.
35. A. Hernandez and M. C. Marchetti, Phys. Rev. E, 2021, 103, 032612.
36. S. Sadhukhan and S. Nandi, eLife, 2022, 11, e76406.
37. M. Czajkowski, D. Bi, M. L. Manning and M. C. Marchetti, Soft Matter, 2018, 14, 5628–5642.
38. B. Ladoux and R.-M. Mège, Nat. Rev. Mol. Cell Bio., 2017, 18, 743–757.
39. J. E. Ron, J. d’Alesandro, V. Cellerin, R. Voituriez, B. Ladoux and N. S. Gov, bioRxiv, 2023, 07.02, 547405.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3sm00627a

---