]contributed equally

]contributed equally

It is widely recognized today that topological defects can play a crucial role in biological contexts, including intracellular cytoskeletal dynamics, tissue growth during embryogenesis, and population-level expansion in bacterial colonies [1, 2, 3, 4, 5, 6, 7]. Two-dimensional (2D) tissues formed by cells reaching confluence are a case of particular interest since they are often precursor of 3D shape. Tissues are often described as active nematics, usually because the elongated shapes of their cells give rise to local nematic order [8, 9]. Such cellular active nematics, like their passive counterparts, can exhibit $\pm\tfrac{1}{2}$ half-integer topological defects, i.e. points in space around which the nematic director turns by $\pm\pi$. Activity endows $+\frac{1}{2}$ defects with an intrinsic velocity, and it is now well-known that cells can be attracted to their core, where they may nucleate additional layers [6, 10] or be extruded [11, 12] or form 3D mounds [13], while they are depleted at the core of $-\frac{1}{2}$ defects [13, 6]

Integer-charge topological defects can also play a role in biological development [14, 15, 16, 17] such as head and foot regeneration in hydra [7] or the growth of plants around their meristems [18]. They do not appear spontaneously in cellular active nematics, essentially for the same ‘energetic’ reasons as in equilibrium liquid crystals [19]. Charge $+1$ defects can, though, be induced either by confinement within small domains  [20, 21, 22], or by gentle large-scale cues imprinted in the underlying substrate [23, 24, 25]. In most cases reported so far, cell accumulation and/or extrusion in 3D has been observed at the core of $+1$ defects.

Models and theories of active nematics have been proposed [2, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35]. Most adopt a ‘top-down’ approach in terms of continuous fields —usually a velocity or polarity field ${\bm{v}}$ and a nematic tensor field ${\bm{Q}}$—, and invoke a contractile/extensile dichotomy. The situation is somewhat confusing: the cells involved are frequently termed contractile. Yet, at the collective level, cellular active nematics are often [13] —but not always [5, 36]— classified extensile from the direction of the intrinsic motion of $+\tfrac{1}{2}$ defects and the flow field around them. These comet-like structures can indeed be found to move ‘toward their head’ like in standard theories with an (extensile) active stress of negative coefficient $-\zeta{\bm{Q}}$ ($\zeta>0$) [37, 2]. Progress toward understanding this conundrum was recently made by invoking the effect of ‘polar fluctuations’ [38, 36, 39, 40].

Our understanding of cell accumulation around the cores of $+1$ defects is also unsatisfactory. In recent experiments on neural progenitor cell (NPC) layers grown on patterned substrates with $+1$ topological defects, large-scale flows toward defect cores leading to cell accumulation were observed irrespective of the specific type of defect considered, be they asters, spirals, or targets (see [25], Fig. 1(d,e), and Fig. 1(a-c) for schematic $+1$ defect types). Current theory, however, predicts cell depletion for targets and strongly chiral spirals.

Here, we adopt a bottom-up approach to build a theory motivated by these experimental observations on NPC active nematics. Following Patelli et al. [41], we first study a particle-level model incorporating the basic ingredients at play, including an external field standing for the guiding patterns, and derive from it a continuous theory. We show that both model and theory account qualitatively for the main experimental result, i.e. accumulation of cells at the core of any type of $+1$ defect. We then argue that cell accumulation is essentially due to two generally ignored anti-symmetric active forces that overcome the standard linear force $-\zeta\nabla\cdot{\bm{Q}}$. We finally discuss the relevance and consequences of our findings in the context of other cellular active nematics experiments and previously proposed theories.

The simplest theory used to describe cellular active nematics [13, 6] usually includes a ‘force balance’ equation

with the coefficient $\gamma$ sometimes replaced by an anisotropic friction $\gamma\to\gamma_{0}(I-\varepsilon\bm{Q})$, where the coefficient $0\leq\varepsilon\leq 1$ quantifies larger friction perpendicular to cell alignment [6].

In the NPC experiments of interest here (see, e.g., Fig. 1(d,e) and [25]), even though shallow, the ridges on the substrate largely impose that the orientation of the nematic field ${\bm{Q}}$ is simply given by that of the $+1$ imprinted defect pattern. There exists a continuous family of $+1$ topological defects parameterized by the tilt angle $\theta_{0}=\theta-\phi$ that their nematic orientation $\theta$ makes with the reference angle $\phi$ of the polar coordinate frame $(r,\phi)$ centered at the defect core (Fig. 1b). In nematic systems, $\theta_{0}\in[-\tfrac{\pi}{2},\tfrac{\pi}{2}]$, asters corresponds to $\theta=\phi$ hence $\theta_{0}=0$, while targets are defined by $\theta_{0}=\pm\tfrac{\pi}{2}$. All other $\theta_{0}$ values correspond to spiral patterns (Fig. 1(a-c)). Writing axisymmetric ${\bm{Q}}$ fields as ${\bm{Q}}=S(r)(\begin{smallmatrix}\cos(2\theta)&\sin(2\theta)\\ \sin(2\theta)&-\cos(2\theta)\end{smallmatrix})$, the ‘active force’ ${\bm{f}}^{a}$ in (1) reads [42]

where $S^{\prime}(r)=\frac{dS(r)}{dr}$, and $\hat{\bm{e}}_{r}$ ($\hat{\bm{e}}_{\phi}$) is the unit vector in the radial (azimuthal) direction. If $\zeta>0$ (“extensile” case), the radial component is inward for $0\leq|\theta_{0}|<\frac{\pi}{4}$ and outward for $\frac{\pi}{4}<|\theta_{0}|\leq\frac{\pi}{2}$ since $S^{\prime}(r)\geq 0$ holds around the core. Arrows in Fig. 1(a-c) represent the direction of force in typical cases. Accumulation will thus occur at the cores of $+1$ defects with $|\theta_{0}|\leq\frac{\pi}{4}$, but the cores of spirals with $|\theta_{0}|\geq\frac{\pi}{4}$ and target cores will be depleted, in contradiction with the experiments, which found cell accumulation in all cases. We further note that anisotropic friction is insufficient to make the target attractive [42].

We now turn to describing the collective dynamics reported in layers of NPCs by a simple particle-level model. When confluent, these cells take elongated shapes, giving rise to local nematic order, and move along their body axis, stochastically reversing their velocity at a relatively short period of 1-2h (compared to experiments that last several days) [13]. These are the ingredients of the Vicsek-style model for dense active nematics introduced and studied in [41], that we complement here by some external field inducing $+1$ topological defects. The positions $\bm{r}_{i}$ and polarity angle $\theta_{i}$ of point particles obey the overdamped discrete-time dynamics

where ${\bm{e}}_{i}$ is the unit vector with orientation $\theta_{i}$, $v_{0}$ is the self-propulsion force, $\epsilon_{i}^{t}$ is a sign reversing with probability $k$ at each timestep, and $\eta$ is the strength of the white uniform angular noise $\chi_{i}^{t}\in[-\frac{\pi}{2},\frac{\pi}{2}]$. Averages $\langle\cdot\rangle_{j}$ in interaction terms are taken over neighboring particles $j$ within unit distance. The first stands for nematic alignment of polarities, and includes the central particle $i$. The second interaction term, of strength $C_{r}$, codes for soft repulsive torques ($\hat{\bm{r}}_{ji}^{t}$ is the unit vector along $\bm{r}_{i}^{t}-\bm{r}_{j}^{t}$). The external field term, of strength $C_{p}$, represents nematic alignment with the fixed defect pattern: the unit vector $\bm{e}_{i}^{p}$ is along the local polar angle $\theta_{i}^{p}=\phi_{i}^{t}+\theta_{0}$ where $\phi_{i}^{t}$ is the azimuthal angle of the $+1$ defect pattern at location ${\bm{r}}_{i}^{t}$.

This model (with $C_{p}=0$) was studied at fixed, moderate, repulsion strength $C_{r}$ in [41]. In order to single out the effects of the imposed $+1$ defect field, we ran simulations in the region where the homogeneous nematic (HN) ordered state is stable [41]. Applying various $+1$ defect patterns, varying the field strength, we found that, typically, values as low as $C_{p}=0.1$ insure that nematic orientation follows the imposed pattern. Figure 2(a-c) shows steady-state radial number density profiles $\rho(r)$ obtained for $C_{p}=0.1$ For all types of $+1$ defects, particle accumulation is observed in the central region, provided $C_{r}$ is neither too small nor too large. We scanned various parameter planes and always found large regions where accumulation occurs. Results for the $(v_{0},C_{r})$ parameter plane are shown in Fig. 2(d-f). (Complementary results can be found in Appendix A, Fig. 4).

Our particle-level model is thus able to account, at least qualitatively, for the systematic accumulation of cells observed in the NPC experiments. We now follow [41], and derive a continuous theory from it. (Details about this derivation can be found in Appendix B and [42].) The starting point is a Boltzmann equation governing the one-body distribution function $f(\bm{r},\theta,t)$ that incorporates a field term inducing nematic alignment with the local orientation $\theta_{p}$ of the imposed pattern:

where $\bm{e}(\theta)$ is a unit vector along $\theta$, $I_{col}$ is the collision integral coding for alignment and repulsion, $a$ is the velocity reversal rate, $\langle\cdot\rangle_{\psi}$ is the average over a noise distribution $P(\psi)$ (with rms $\eta$). This equation is then truncated and closed using a scaling ansatz. The result is a set of coupled equations governing the first three angular Fourier modes $f_{k}(\bm{r},t)=\int_{-\pi}^{\pi}f(\bm{r},\theta,t)e^{ik\theta}d\theta$ of $f$:

where $\triangledown\equiv\partial_{x}+i\partial_{y}$ and $\triangledown^{*}\equiv\partial_{x}-i\partial_{y}$. While $f_{0}$ is real and identical to the density field $\rho$, $f_{1}$ and $f_{2}$ are complex-valued, representing density-weighted polarity/velocity and nematic order fields

All coefficients in Eqs. (6-8) depend on the microscopic-level parameters (see [42] for explicit expressions). These equations are identical to those derived in [41] except for the last term of (7) and (8) which involves the external patterned field.

Equations (6-8) can be simulated as is, and their solutions have been shown to be qualitatively faithful to the particle-level model (in the absence of external field) [41]. We have run simulations with imposed $+1$ defect patterns, again using parameters mostly in the stable HN ordered state region, and found that accumulation at defect cores occurs even for small $C_{p}$ values. Yet, with their many terms, Eqs. (6-8) do not easily allow a deeper understanding of the features at the origin of accumulation at $+1$ defect cores. They can be simplified rather straightforwardly. First of all, as seen in simulations, $|f_{1}|$ is always much smaller than $|f_{2}|$ in the regimes under consideration here where nematic order is well developed almost everywhere. This suggests to neglect all terms containing $f_{1}$ twice. We also found, in simulations with an imposed field pattern, that canceling the $f_{2}\triangledown^{*}\rho$ term in (7) and the $\kappa_{1}$ and $\kappa_{2}$ terms in (8) does not impact much the behavior of the system. Our simplified Eqs. (6-8) thus read, using tenso-vectorial notations:

where $\tilde{\bm{Q}}^{p}$ is the nematic field of the imposed pattern, $\tilde{S}=\rho S$, $\bm{E}$ is the symmetric strain rate-like tensor defined by $E_{ij}=\frac{1}{2}(\partial_{i}w_{j}+\partial_{j}w_{i}-\delta_{ij}\nabla\cdot\bm{w})$, $\bm{\Omega}$ is the antisymmetric rate of rotation tensor $\Omega_{ij}=\frac{1}{2}(\partial_{i}w_{j}-\partial_{j}w_{i})$, and $\mho_{ij}=w_{j}\partial_{i}-w_{i}\partial_{j}$. Note that Eq. (12) is similar to the Beris–Edwards equation [43, 44, 27, 35]

Numerical simulations of Eqs. (10-12) and of Eqs. (6-8), performed at similar parameters, yield very similar results ¹¹1One quantitative difference is that a smaller $C_{p}$ value is sufficient, for the simplified equations, to fully impose the external pattern and in particular to avoid that the $+1$ defect splits into two $+\tfrac{1}{2}$ ones.. We now show features of the typical axisymmetric steady-state solutions of Eqs. (10-12) when imposing weak patterns ($C_{p}=0.1$). In Fig. 3(a), where a target pattern is induced, radial density and velocity profiles reveal that inward flows and accumulation in the core region are only replaced by outward flow and depletion for large enough nominal speed $v_{0}$. Figure 3(b,c) display phase diagrams in a parameter plane equivalent to that of Fig. 1(d-f) for both target and aster patterns. Showing large regions of accumulation near defect cores, they are very similar to those obtained at particle level.

To obtain some analytical insights and identify the roles of active forces around $+1$ defect, we make further simplifications. First of all, we neglect the $\rho$ dependence of coefficients ($\alpha,\mu,\pi_{0},\pi_{1},$ and $\zeta$ are impacted). In the steady state, Eq. (Integer Topological Defects Reveal Anti-Symmetric Forces in Active Nematics) yields

When $\tilde{\bm{Q}}$ is well aligned with the external field in the steady state and $\bm{w}$ is negligible compared with $\tilde{\bm{Q}}$ in Eq. (12), $[\mu-\xi\tilde{S}^{2}]\bm{\tilde{Q}}+C_{p}\bm{Q}^{p}\rho=0$ holds. Writing $\tilde{\bm{Q}}=\rho S{\bm{Q}}^{p}$ yields $[\mu-\xi\rho^{2}S^{2}]S+C_{p}=0$. Solving for $\rho S$ assuming small deviations from $\sqrt{\mu/\xi}$ ($\mu,\xi>0$), we obtain $\rho S\sim\sqrt{\mu/\xi}+\frac{\rho C_{p}}{2\mu}+\mathcal{O}(C_{p}^{2})$. Injecting this expression in $\Gamma$ yields $\Gamma(\tilde{\bm{Q}},{\bm{Q}}^{p})=\Gamma(\tilde{\bm{Q}})=\gamma_{0}(\bm{I}-\varepsilon\tilde{\bm{Q}}).$ where $\gamma_{0}=\alpha+\beta\tilde{S}^{2}>0$ and $\varepsilon=[\sigma+C_{p}/(2\rho S)]/\gamma_{0}$. ($\alpha>0,\beta>0$) In [42], we argue that, for axisymmetric solutions, the anisotropy coefficient $\varepsilon$ (in $\Gamma(\tilde{\bm{Q}})$) cannot reverse the sign of the radial component of ${\bm{w}}$. Thus, this radial flow vanishes (steady-state condition) when $f_{r}=\pi_{0}d\rho/dr$ where $f_{r}$ is the radial component of ${\bm{f}}^{a}$. Consequently, accumulation or depletion, governed by the sign of $d\rho/dr$, is decided by the sign of $f_{r}$.

We also show in [42] that $f_{r}$, when a target or aster pattern is imposed, can be approximated by:

where the $+$ sign is for a target, and the $-$ sign is for an aster. We thus conclude that under our approximations, the accumulation-depletion transition is given by $f_{r}=0$. On the transition line, we can safely assume that $\rho=\rho_{0}$. For weak fields ($C_{p}\ll 1$) and away from the defect core, this gives $\pm\zeta-(\gamma_{1}-\gamma_{2})\sqrt{\mu(\rho_{0})/\xi}=0$. The corresponding lines, indicated by the solid lines in the phase diagrams of Fig. 3(b,c), match very well the numerical results for Eqs.(9-11) presented there. On the other hand, the transition lines $\zeta=0$ (the dashed lines) given by considering only the linear force ($\gamma_{1}=\gamma_{2}=0$), are unable to account for the accumulation-depletion transition.

Our results suggest a novel classification of cellular active nematics. From Eq. (16), we see that $\gamma_{1}-\gamma_{2}$ is a crucial quantity to decide whether $+1$ defects will lead to accumulation or depletion. This comes in addition to the sign and value of $\zeta$, the coefficient of the standard linear force term. The plane $(\zeta,\gamma_{1}-\gamma_{2})$ can be divided into 4 quadrants classifying various responses to the presence of $+1$ defects. Sketched in Fig. 3(d), it shows how the traditional extensile-contractile dichotomy, governed by the sign of $\zeta$, is unable to account alone for the phenomena at stake. Our bottom-up approach yields expressions for these coefficients, and shows that $\gamma_{2}<0$ holds for the systems with repulsive interactions, while $\gamma_{1}>0$ in most cases (see Appendix B), so that $\gamma_{1}-\gamma_{2}>0$ for the NPC layers of interest here, which appear to be located in zone II. We confirmed further the validity of our conclusion by varying $\zeta$ and $\gamma_{1}-\gamma_{2}$ independently of the coefficients derived from the particle model, considering also the case of $-1$ defects (see Appendix C).

Finally, to grasp the intuitive meaning of the $\gamma_{1}$ and $\gamma_{2}$ terms, we rewrote them assuming $S=1$ and uniform density. They can then be recombined, together with the usual $\zeta$ active stress, into classic bend and splay terms:

where $A=-2\zeta$ and $B=2(\gamma_{1}-\gamma_{2})$, and $\bm{n}=(n_{x},n_{y})$ is the director field (see [42]). The classification diagram (Fig. 3(d)) can then be re-expressed as in Fig. 3(e). In this diagram, the horizontal axis is pure splay, and the vertical axis is pure bend. As long as the anti-symmetric term exists, $\gamma_{1}-\gamma_{2}\neq 0$, the one-constant approximation is not valid. Thus, there are always situations where the behavior near integer defects cannot be classified extensile (I) or contractile (III).

To summarize, we have demonstrated the importance of the two generally ignored anti-symmetric forces and shown that the conventional active stress term alone cannot account for the accumulation towards the core of $+1$ topological defects in cellular active nematics such as the NPC layers of interest here. All in all, we believe our results point to the importance of nonlinear terms often absent from theories of 2D active nematics systems.

Appendix A: Complementary results at particle level.— Phase diagrams similar those presented in Fig. 2(d-f), but in various other parameter planes are presented in Fig. 4. For each parameter plane shown, the top panel is for the target pattern, with the bottom panel for the aster pattern.

Appendix B: Complementary information about continuous theory.— Here we give some important points in the derivation of the continuous theory (6-8) that is detailed in [42] and follows [41]. In the Boltzmann equation (5), the collision integral reads:

where $\Psi$ is the $\pi$-periodic nematic alignment function ($\Psi(\theta_{1},\theta_{2})=\frac{1}{2}(\theta_{2}-\theta_{1})$ for $-\frac{\pi}{2}<\theta_{2}-\theta_{1}<\frac{\pi}{2}$), and $K$ is the kernel function

which codes for the repulsive torques between particles, depends on the relative position $\mathbf{r^{\prime}}\!-\!\mathbf{r}\equiv s\mathbf{e}(\phi)$ of the two colliding particles, and incorporates the distance-dependent weighting function $g(s)$.

Upon expanding $f(\mathbf{r}\!+\!s\mathbf{e}(\phi),\theta_{2})$ at first order in $s$, and Fourier transforming (5) into a hierarchy for the fields $f_{k}$, $g(s)$ appears only via the first two of its moments $b_{n}=\int_{0}^{s}s^{n+1}g(s)ds$: $b_{0}$ that is set to unity in all generality, and $b_{1}$ that is thus the only microscopic-level coefficient controlling the strength of repulsive torques, akin to $C_{r}$ in the particle-level model.