Chiral stresses in nematic cell monolayers

Abstract

Recent experiments on monolayers of spindle-like cells plated on adhesive stripe-shaped domains have provided a convincing demonstration that certain types of collective phenomena in epithelia are well described by active nematic hydrodynamics. While recovering some of the hallmark predictions of this framework, however, these experiments have also revealed a number of unexpected features that could be ascribed to the existence of chirality over length scales larger than the typical size of a cell. In this article we elaborate on the microscopic origin of chiral stresses in nematic cell monolayers and investigate how chirality affects the motion of topological defects, as well as the collective motion in stripe-shaped domains. We find that chirality introduces a characteristic asymmetry in the collective cellular flow, from which the ratio between chiral and non-chiral active stresses can be inferred by particle-image-velocimetry measurements. Furthermore, we find that chirality changes the nature of the spontaneous flow transition under confinement and that, for specific anchoring conditions, the latter has the structure of an imperfect pitchfork bifurcation.

---

I. Introduction

Multicellular systems of prokaryotes, such as suspensions of planktonic bacteria, have historically played a pivotal role in the development of the hydrodynamics of active fluids, since the early work of Batchelor on the stress distribution in suspensions of microswimmers.¹ By contrast, multicellular systems of eukaryotes have entered only recently into the realm of active hydrodynamics, following a number of inspiring experimental works on epithelial and mesenchymal cell layers and tissues (see e.g.ref. 2–8). Among these, monolayers of clotured spindle-like cells, such as NIH 3T3 mouse embryo fibroblasts,⁴ murine neural progenitor cells (NPCs),⁶ human bronchial epithelial cells (HBEC),⁷ Retinal Pigment Epithelial (RPE1) cells⁸ and C2C12 mouse myoblasts⁸ represent an especially promising class of model systems, because of their connection with active nematic liquid crystals (see e.g.ref. 9).

First identified as a broken symmetry in certain types of cell cultures,¹⁰ and later exploited to decipher their static¹¹ and dynamical properties,^2,4–8 nematic order has surged as one of the central themes in collective cell dynamics. In layers of spindle-like cells, where the local orientation can be unambiguously identified, nematic order is marked by the presence of ±1/2 disclinations,^4,6i.e. point-like singularities about which the average cellular orientation rotates by ±π. Consistently with the predictions of active nematic hydrodynamics,^12,13 these defects self-propel and pairwise annihilate until cell crowding freezes the system into a jammed configuration. Before dynamical arrest, the collective motion of the cells gives rise to a decaying turbulent flow at low Reynolds number, whose statistics, spatial organization and spectral structure are in exceptional agreement with the hydrodynamic picture¹⁴ (but see ref. 15 for an alternative theoretical picture based on glassy dynamics).

Another remarkable demonstration of active hydrodynamic behavior in eukaryotic cell layers has been recently reported by Duclos et al., upon confining spindle-like RPE1 and C2C12 cells within adhesive stripe-shaped domains.⁸ Depending on the width of the stripe the system was found either in a stationary state, with the cells parallel to the longitudinal direction of the confining region, or in a collectively flowing state characterized by a spontaneous tilt of the cells toward the center of the stripe. The latter picture, often referred to as spontaneous flow transition, had been anticipated for over a decade by Voituriez et al.¹⁶ and represents one of the hallmarks of active liquid crystals. While confirming this seminal prediction, however, Duclos et al. have also highlighted a number of unexpected features that could be ascribed to the existence of chirality over length scales larger than the typical size of a cell.

The notion of chirality is not new in active matter and has been theoretically explored well before the interest around collective cell dynamics in eukaryotes had started to blossom. Fürthauer et al., for instance, demonstrated that microscopic torque dipoles, such as those arising from rotating molecular motors or flagella, give rise to antisymmetric stresses and angular momentum fluxes, which, in turn, drive rotating flows and other chiral patterns on the large scale.^17,18 More recently, Banerjee et al. showed that rotational motion at the microscopic scale further enriches the spectrum of hydrodynamical behaviors or chiral active fluids by giving rise to non-dissipative “odd” viscosity,¹⁹ analogous to that found in quantum Hall fluids.²⁰ Whereas undoubtedly interesting and relevant for a broad class of biological and synthetic systems, these mechanisms appear however unsuited to account for the chirality observed in the experiments by Duclos et al., because of the manifest lack of rotational motion at the scale of individual cells.

In this article, we show that macroscopic chirality can arise in nematic cell layers as a consequence of a misalignment between the cell's local orientation and active forces, even in the absence of microscopic rotational motion (Section II). Collectively, this gives rise to a chiral and yet symmetric stress tensor that, in two dimensions, complies with the symmetries of the nematic phase. Next, we explore the effect of such a chiral stress on the active flow generated by ±1/2 disclinations and identify a characteristic signature of chirality from which the ratio between chiral and non-chiral active stresses can be experimentally estimated (Section III). Finally, following Duclos et al.,⁸ we investigate the hydrodynamic stability of a chiral nematic cell monolayer confined on adhesive stripes and subject to various boundary conditions and classify all possible scenarios arising from the interplay between the geometry of the confining region, the extensile/contractile stresses and chirality (Section IV).

II. Chiral stresses

Let us consider a two-dimensional volume element and let n be the nematic director representing the average direction of the enclosed cells (Fig. 1a). The most generic form of the stress tensor associated with such a volume element can be expressed in the basis of the nematic director n^⊥ and its normal vector n^⊥, namely:

σ^a = σ_‖nn + σ_⊥n^⊥n^⊥ + τ(nn^⊥ + n^⊥n). (1)

Here σ_‖ and σ_⊥ represent the stresses experienced by the volume element in the direction parallel and perpendicular to n, whereas τ is the shear stress. By construction, σ^a is invariant under the transformation n → −n, thus is consistent with the symmetry of the nematic phase, but is evidently chiral for any non-zero τ value, as the stress distribution depicted in Fig. 1a does not coincide with its mirror image. Furthermore, taking n^⊥n^⊥ = I − nn, with I the identity tensor, allows one to cast the active stress as the sum of an isotropic contribution, equivalent to an active pressure, a deviatoric term, common to all active nematic liquid crystals regardless of their chirality,^21,22 and a chiral term, namely:where:Microscopically, the chiral stress τ might originate from the fact that the force exerted by an individual cell is tilted with respect to the cell orientation. To illustrate this concept let us consider an individual cell whose major and minor axes are parallel to the unit vectors ν_c and ν^⊥_c (Fig. 1b). Following Lau and Lubensky,²³ one can express the stress tensor as σ^a = Σ_cd_cδ(r − r_c), where the index c runs over all the cells in the system, r_c is the position of the c-th cell andis a force–dipole tensor. In eqn (4), R_c is the distance between the cell's surface and center of mass, −dAf_c the force exerted by the cell's area element dA on the surrounding medium and the integral is calculated over the surface Σ_c of the c-th cell. A complete derivation can be found in ref. 23. For an effectively two-dimensional system, such as the one considered here, dA = d[small script l]h, with h the thickness of the cell along the z-direction, here assumed to be uniform throughout the sample, and [small script l] the one-dimensional arc-length.

Fig. 1 Schematic representation of stresses and forces in cellular nematic monolayers. (a) A volume element whose faces are conventionally oriented in the direction of the nematic director n and its normal vector n^⊥. In the most general setting, the volume element is subject to three independent stresses: the two normal stresses σ_‖ and σ^⊥ and the shear stress τ. (b) The chiral stress τ arises at the microscopic level when the force exerted by an individual cell is neither parallel nor perpendicular to the cell axis.

Now, the magnitude of the active stresses P^a, α, and τ depend exclusively upon the distribution of the forces exerted by the cells along their contour. The simplest approximation of the force density field f_c consists then of a dipole of the form:

f_c = F_cδ(R_c − aν_c) − F_cδ(R_c + aν_c), (5)

where a is the cell's major semi-axis (Fig. 1b and ref. 22). To make progress, we express the cellular forces in the {ν_c,ν^⊥_c} basis, namely F_c = F_‖ν_c + F_⊥ν^⊥_c, with F_‖ and F_⊥ the longitudinal and transverse components of the force exerted by the c-th cell, here assumed for simplicity uniform throughout the system. Replacing this in eqn (5) and coarse-graining σ^a over the length scale of a volume element Ω(r) centered at r (Fig. 1a) yields:

σ^a = aρ[2F_‖〈ν_cν_c〉 + F_⊥〈ν_cν^⊥_c + ν^⊥_cν_c〉], (6)

where 〈⋯〉 is the spatial average within Ω(r) and ρ is the local cell number density. Finally, calculating the averages and comparing with eqn (2) readily gives:

P^a = −aρF_‖, α = 2aρSF_‖, τ = aρSF_⊥, (7)

where S = 2〈|ν_c·n|²〉 − 1 is the local nematic order parameter in two dimensions. As expected, in the absence of nematic order (i.e. S = 0), the active forces exerted by the cells result exclusively in an effective pressure, while both the deviatoric and chiral stress vanishes identically. When nematic order is not uniform throughout the system, hence S varies in space, the active stress tensor of eqn (2) is more conveniently expressed in terms of tensor order parameter Q_ij = S(n_in_j − δ_ij/2), namely:

σ^a_ij = −P^aδ_ij + α₀Q_ij − 2τ₀ε_ikQ_kj, (8)

where α₀ = 2aρF_‖ and τ₀ = aρF_⊥ are constants independent on the nematic order parameter S and ε_ij is the two-dimensional Levi-Civita tensor (i.e. ε_xx = ε_yy = 0 and ε_xy = −ε_yx = 1).

Some comments are in order. Although eqn (5) is only a rudimental approximation of the force field generated by an irregularly-shaped cell, considering a more involved force distribution does not change the qualitative picture with respect to the emergence of chiral stresses, as long as this is asymmetric with respect to the cell's longitudinal direction. To illustrate this concept, we discuss in Appendix A the case of a quadrupolar force distribution. Whereas the exact origin of this asymmetry is beyond the scope of the present article, the biophysical literature is not scarce of examples where chirality can be detected at the single-cell level. For instance, various mammalian cells, when plated on micropatterns, can break the left-right symmetry by suitably positioning their internal organelles with respect to the cell body.²⁴ Analogously, chirality can emerge at the scale of the entire cell from the self-organization of the actin cytoskeleton.^25,26 The broken symmetry can furthermore propagate over the mesoscopic scale and bias the cell's collective migratory motion.²⁷

III. Defect motion

The motion of ±1/2 disclinations has become a hallmark of active nematic liquid crystals. As it was theoretically predicted^12,13 and experimentally verified in both microtubules-kinesin^28–31 and actin-myosin suspensions,³² ±1/2 disclinations in active nematics drive a flow, whose structure and speed is strictly related to the geometry of the defect. For +1/2 disclinations, in particular, this flow has a Stokeslet structure that results in a propulsion of the defect in the direction of its symmetry axes p or −p depending on whether the system is contractile or extensile (Fig. 2 and ref. 33). Exceptionally, the same mechanism has been experimentally identified in various types of cell cultures, including spindle-shaped NIH 3T3 mouse embryo fibroblasts,⁴ murine neural progenitor cells (NPCs),⁶ Madin Darby canine kidney cells (MDCKs)⁵ and human bronchial epithelial cell (HBEC).⁷ Unlike in suspensions of rod-like cytoskeletal filaments, however, the direction of motion of +1/2 disclinations is not necessarily parallel to ±p and also the reconstructed flow around the defects is less symmetric than that found in cytoskeletal suspensions. Whereas it is not unlikely for such a feature to originate from statistical errors, here we demonstrate that chiral active stresses affect the dynamics of +1/2 defects precisely by rotating their direction of motion with respect to p. In fact, the angle between p and the direction of motion of the defects could be used to indirectly measure the relative magnitude of the chiral stress τ compared to the deviatoric stress α. This phenomenon shares some similarities with recent results by Maitra and Lenz about the dynamics of +1/2 disclinations in rotating active nematics.³⁴

Fig. 2 Pressure (a and b) and velocity (c and d) fields in proximity of ±1/2 defects obtained from the analytical solutions of eqn (9) for an extensile chiral active nematic with α = 4τ < 0. The configuration of the director is indicated by white lines for the case of +1/2 (a) and −1/2 (b) defects. For +1/2 defects, the velocity field at the core, thus the direction of motion of the defect, is tilted by an angle θ_tilt = arctan(1/2) ≈ 27° with respect to the defect polarity direction p = x̂. For both chiral and achiral +1/2 defects, the pressure is anisotropic and larger toward the direction of motion of the defects.

An analytical approximation of the flow driven by the active stresses in the surrounding of a disclination can be obtained by solving the incompressible Stokes equation with a body force resulting from the active stress associated with an isolated ±1/2 defect. Namely:

η∇²v + f^a_± = ∇P, ∇·v = 0, (9)

where η is the shear viscosity of the tissue, here assumed isotropic for simplicity, and f^a_± = ∇·σ^a is the body force arising from spatial variations of the active stress of eqn (2) in the presence of a ±1/2 defect. This approach, introduced in ref. 13, does not account for feedback of the flow on the orientation of the director and hence the structure of the defect. Nevertheless it provides a simple and faithful approximation of the defective flow as well as an estimate of the defect velocity. A generic solution of eqn (9) can be expressed in the form v = v_± + v₀, where v₀ is a solution of the homogeneous Stokes equation and is dictated by the boundary conditions and v_± is given by:where G is the two-dimensional Oseen tensor given by:with a constant set by the boundary conditions. Analogously, taking the divergence of eqn (9) allows to express the pressure field as the solution of the following Poisson equation:

∇²P_± − ∇·f^a_± = 0. (12)

Now, in the presence of a disclination of strength s = ±1/2 located at the origin of the (x,y)-plane and oriented in the direction p = (cos ψ,sin ψ), the nematic director n = (cos θ,sin θ) has local orientation θ = sϕ + (1 − s)ψ, with ϕ = arctan(y/x) the polar angle.³³ The body force f^a_± is then readily found in the form:where p^⊥ = (−p_y,p_x) and r is the distance form the defect core. The effect of chirality is most dramatic for +1/2 disclinations. In achiral active nematics, τ = 0 and the active force f^a₊ is purely longitudinal. For non-zero τ, f^a₊ acquires a transverse component resulting in a tilt in the direction of motion of the defect by an anglewith respect to the orientation p (Fig. 2). As anticipated, eqn (14) can in principle be used in combination with experimental reconstruction of defect trajectories in order to estimate the relative magnitude of the chiral and deviatoric stresses in nematic cell monolayers.

To calculate the flow velocity in proximity of a ±1/2 defect, we set, without loss of generality, ψ = 0, thus p = x̂ and p^⊥ = ŷ and we assume the defect at the center of a circular domain of radius R. The exact velocity of the flow at the boundary of such a domain, hence the homogeneous solution v₀, is not relevant for the purpose of this discussion. In practice, this will be determined by the chemistry of the substrate and the possible presence of other topological defects in the same region.⁴ Under these assumptions, and carrying out algebraic manipulations as those in ref. 13, the flow velocity caused by a ±1/2 defect can be found from eqn (10), (11) and (13) in the form:

These flows are illustrated in Fig. 2c and d for a specific choice of the angle θ_tilt. The corresponding pressure field is readily found from eqn (12) with ∇·f^a₊ = −(α cos ϕ + 2τ sin ϕ)/(2r²) and ∇·f^a₋ = 3(α cos 3ϕ + 2τ sin 3ϕ)/(2r²). This yields:where P⁽⁰⁾_± are harmonic functions depending on the boundary conditions. Interestingly, the pressure field given by eqn (16) is independent on the distance from the defect core, but varies with the angle ϕ and, depending on the sign of the stresses α and τ, is maximal or minimal at specific directions relatively to the polarity vector p. For instance, for extensile systems (i.e. α < 0) with negligible chirality (i.e. τ ≈ 0), eqn (16a) predicts a pressure maximum in the −p direction, thus toward the “head” of the comet-like structure characteristic of +1/2 defects. Despite our analysis revolving around incompressible systems, we can expect this trend to persist in the presence of small spatial variations of the density field ρ. In this case, and under the assumption that ρ ∼ P, we expect a higher density of extensile (contractile) cells in the front (back) of +1/2 defects, consistent with the experimental observation of Kawaguchi et al.⁶ Finally, the lack of spatial dependence in eqn (16) originates for the specific parametrization of the orientation field θ in proximity of the defects, i.e. θ(ϕ) = ±ϕ/2. In a more realistic setting, θ(r,ϕ) = ±ϕ/2 + θ_far(r,ϕ), where θ_far is a non-singular function that determines the far-field configuration of the nematic director and vanishes at the defect core. Accounting for the far-field configuration of the director results into a dependence of the pressure on the distance from the defect core.

In summary, the presence of a symmetric chiral active stress, such as that embodied by the parameter τ, affects the flow generated by ±1/2 disclinations by stretching and rotating the velocity field in the surrounding of the defects (Fig. 2c and d). Most prominently, this results in a tilt in the direction of motion of +1/2 defects: i.e.v_self = v₊(r = 0) = R/(4η)(αp + 2τp^⊥). Thus +1/2 defects self-propel at an angle θ_tilt with respect to their orientation p [see eqn (14)]. Such an angle could in principle be measured in experiments on two-dimensional cell cultures, thus providing a direct measurement of the relative magnitude of the chiral stress. The same behavior has been reported in the case of actively rotating +1/2 defects, in the limit of vanishing angular velocity.³⁴

IV. Spontaneous flow on adhesive stripes

In this section we revisit a classic problem of the hydrodynamic stability and spontaneous flows of active nematics in a quasi-one-dimensional channel (Fig. 3). First discussed in a seminal paper by Voituriez et al.,¹⁶ later elaborated by many others^35–37 and recently observed in experiments on cell monolayers⁸ and suspensions of microtubules and kinesin,³⁸ this phenomenon consists of a continuous transition between a stationary and uniformly oriented configuration to a state characterized by a spontaneous distortion of the nematic director coupled to an internally driven shear flow. The transition, in many aspects similar to the Fréedericksz transition in passive liquid crystals,³⁹ results as a consequence of two different mechanisms. First, a distortion of the nematic director drives a shear flow as illustrated by eqn (9) for the time-independent case; second, the nematic director rotates in a shear flow. As a consequence, when the hydrodynamic torque driven by the active stresses outweighs the elastic restoring torque, the uniformly oriented configuration becomes unstable to splay or bending deformations, depending on the sign of the active stress α and other material parameters. Roughly speaking, this occurs when the active length scale , defined by the ratio between the passive and active torques, with K the Frank elastic constant, becomes comparable to the width L of the channel.¹⁴

Fig. 3 Schematic representation of the channel of infinite length in x-direction and width L in y-direction.

In the following, we extend and generalize the theoretical analysis by Duclos et al.⁸ by considering various experimentally relevant scenario in terms of boundary anchoring and flow. The hydrodynamic equations governing the dynamics of an incompressible (i.e. ∇·v = 0) active nematic liquid crystal are given by:^16,36,37

where D/Dt = ∂_t + v·∇ is the material derivative, u_ij = (∂_iv_j + ∂_jv_i)/2 and ω_ij = (∂_iv_j − ∂_jv_i)/2 are respectively the strain-rate and vorticity tensor and h = −δF/δn is the so-called molecular field, governing the relaxational dynamics of the nematic director, with F the Frank free energy. In one-elastic-constant approximation, the latter is simply given by and h = K∇²n. The material parameters λ and γ are respectively the flow-alignment parameter, controlling the tendency of the nematic director to rotate in a shear flow, and the rotational viscosity of the nematic fluid. In eqn (17b), the pressure P includes the active pressure P^a given in eqn (2) and (7), but, due to incompressibility, is not an independent field and σ^e is the elastic stress resulting from the departure of the director configuration from the ground state of the Frank free energy. Although it does affect the onset of the spontaneous flow instability, the latter is often unimportant for the phenomenology of active nematics and will be neglected here for simplicity.

We solve eqn (17) in a rectangular channel which is infinitely long in x-direction and has width L in y-direction (Fig. 3). Because the channel is invariant for translations along the longitudinal direction, we assume for simplicity n and v to be independent of x, both while stationary and spontaneously flowing. Furthermore, incompressibility and mass conservation demand the y-component of the velocity to vanish identically, i.e. v_y = 0. Thus taking n = (cos θ,sin θ) and v = (v_x,0), the eqn (17) reduce to:

whereas from the y-component of eqn (17b) we obtain the following expression for the pressure field:with P₀ a constant. Next, we look for stationary solutions of eqn (18) subject to different boundary conditions in terms of orientation of the nematic director at the boundary and whether or not the cells are allowed to slide along the channel walls, thus: θ = θ(y), v_x = v_x(y). In order for the fluid to be stationary, the shear stress must be uniform across the channel. Thus, from eqn (18b):Solving eqn (20) with respect to ∂_yv_x and substituting this into eqn (18b) yields a single homogeneous equation for θ, namely:Before analyzing specific cases in detail, we consider the generic situation in which the nematic director is anchored to the channel walls at an arbitrary angle θ₀ ≠ arccos(1/λ)/2, thus:

θ(0) = θ(L) = θ₀. (22)

In the following, we will separately analyze the scenarios in which the cells are stationary at the boundary of the channel (Section IV.A) and when, on the other hand, they are able to slide while keeping their orientation fixed (Section IV.B).

Our analysis is complemented by numerical solutions of eqn (18) with various boundary conditions. For this purpose we rescale time by the viscous time scale τ_ν = ρL²/η, length by the channel width L and stress by the viscous stress scale σ_ν = ρL²/τ_ν², i.e. t → t/τ_ν, y → y/L, σ → σ/σ_ν. All the other quantities in Fig. 4 and 5 are rescaled accordingly.

Fig. 4 Bifurcation diagram of the spontaneous flow transition obtained from numerical (dots) and analytical (lines) solutions of eqn (18) for K/γ = L = 1, λ = −0.5 and various τ values. (a) No-slip boundary conditions and parallel anchoring (see Section IV.A.1). The chiral stress τ does not influence the critical α value, but weakly affect the post-transitional configuration of the nematic director. (b) Stress-free boundary conditions and special boundary anchoring θ(0) = θ(L) = −θ_tilt/2 (see Section IV.B.3). The chiral active stresses, embodied by the parameter τ, explicitly break the clock-counterclockwise symmetry of the lowest free-energy configuration rendering the pitchfork bifurcation “imperfect”. In this case, only one of the two branches of the bifurcation diagram is connected to the trivial solution, which may then be the only one observed experimentally.

Fig. 5 Numerical solutions of eqn (18) in the dimensionless quantities defined above and for K/γ = 1, λ = −0.5 as well as different values of τ and α for three sets of boundary conditions. (a) and (b) are numerical solutions for the eqn (18) with parallel anchoring as well as no-slip boundary conditions [panel (a) and analyzed analytically in Section IV.A.1] and stress-free boundary conditions [panel (b) and Section IV.B.1], respectively. (c) Numerical solutions for stress-free walls and the stationary solution θ₀ = −θ_tilt/2 = −arctan(2τ/α)/2 imposed at the boundary, see Section IV.B.3.

A. No-slip boundary conditions

In this section, we consider the case in which the cells are unable to slide along the boundary of the channel, i.e. v_x(0) = v_x(L) = 0, which, in turn, experience a non-vanishing stress σ_xy ≠ 0 resulting from the cellular forces. In this case, θ(y) = θ₀ is always a trivial solution of eqn (21) with the boundary condition given by eqn (22), and the monolayer admits a stationary and uniformly aligned configuration. Because of the internal stresses, however, such a uniform state can become unstable with respect to splay or bending deformations for sufficiently large active stresses or channel widths. To illustrate this point, we take θ(y) = θ₀ + δθ(y), with δθ(0) = δθ(L) = 0, and linearize eqn (21) about δθ = 0. This yields:

∂_y²δθ + q²δθ = 0, (23)

where we have introduced the constant:The solution of eqn (23) is readily found to be:

δθ = C sin(qy), qL = nπ, (25)

with a constant C and n ∈ [Doublestruck Z] an arbitrary integer. By virtue of eqn (20), the corresponding velocity field is given by:

η∂_yv_x = δθ(2τ sin 2θ₀ − α cos 2θ₀). (26)

whose solution with no-slip boundary conditions is given by:with m ∈ [Doublestruck Z] another arbitrary integer. Upon comparing eqn (25) and (27) we find that the first mode to be excited is (n,m) = (2,1), thus the trivial solution θ(y) = θ₀ becomes unstable when q = q_c = 2π/L or, equivalently, when L = L_c = 2π/q. To be more specific, we consider, in the following, two practically relevant cases, where θ₀ = 0 (parallel anchoring) and θ₀ = π/2 (homeotropic anchoring).

1. Parallel anchoring. If the nematic director is parallel to the channel walls θ₀ = 0, thus:As in non-chiral active nematics, the instability is triggered by splay deformations (i.e. transverse to the nematic director) and uniquely depends on the non-chiral active stress α. Furthermore, as in non-chiral active nematics,¹³ such a splay instability affects flow-aligning systems (i.e. λ > 1) in the presence of extensile active stresses (i.e. α < 0), and flow-tumbling systems (i.e. λ < 1) in the presence of contractile active stresses (i.e. α > 0). A critical α value is readily found in the form: α_c = 8π²ηK/[γL²(1 − λ)].

At the onset of the transition, the constant C can be calculated upon expanding eqn (21) up to the third order in δθ, Then, using the solution of the linearized equation yields a cubic equation in C. Solving the latter gives (see Appendix B):

The spontaneous flow instability consists, therefore, of a standard pitchfork bifurcation whose relevant fields, θ and v_x, scale like (L − L_c)^1/2 at criticality, see Fig. 4a. Despite the fact that the chiral stress τ does not affect the instability of the stationary state, it leaves a clear signature on the post-transitional behavior of the flowing monolayers. This can be seen in Fig. 5a, showing numerical solutions of eqn (18) in the flowing state for various α and τ values. The most prominent effect of chirality, in this case, is evidently to render both the distortion of the nematic director and the associated flow asymmetric with respect to the channel centerline.

2. Homeotropic anchoring. The case in which the nematic director is perpendicularly anchored to the channel walls, θ₀ = π/2, yields:In this case the instability is triggered by bending deformations (i.e. parallel to the nematic director). In contrast with the scenario of Section IV.A.1, flow-aligning systems are unstable in the presence of extensile active stresses, whereas strongly flow-tumbling systems (i.e. λ < −1), are unstable in the presence of contractile active stresses. The critical α value is readily found in the form: α_c = −8π²ηK/[γL²(1 + λ)].

To conclude this section, we observe that for both parallel and homeotropic anchoring, the spontaneous flow instability crucially relies on the flow-alignment behavior of the system, governed by the phenomenological parameter λ. As for molecular liquid crystals, where λ depends upon the molecules shape and interactions, we expect the flow-alignment parameter to be affected by the cellular shape, which in turn is not fixed, and by the passive and active processes underlying the cell–cell and cell–substrate interactions.

B. Stress-free boundary conditions

In this subsection we consider the scenario in which the cells are allowed to slide along the boundary, while keeping a fixed orientation θ₀ with respect to the channel walls. As in the adhesive stripes used in ref. 8, the channel walls do not comprise a real physical barrier, but represent instead the interface between two regions of the substrate with different coating. As the walls now do not exert any force on the cells σ_xy(0) = σ_xy(L) = 0. Mechanical equilibrium [i.e.eqn (20)] thus implies σ_xy = 0 everywhere.

The most striking difference with respect to the case discussed in Section IV.A, as well as the most prominent consequence of the chiral stress τ, is that the stationary and uniformly aligned configuration (i.e. θ = θ₀ and v_x = 0) is not a trivial solution of eqn (21) for non-vanishing α and τ values, unless the chiral and non-chiral active stresses cancel each other identically. Before considering this latter case (see Section IV.B.3), we find approximated expressions for the local orientation θ and the velocity v_x in the limits in which the active stresses are either very small or very large.

For very large active stresses, the active terms in eqn (21) overweight the elastic term on the left hand side. As a consequence, the equilibrium configuration of the nematic monolayer consists of a region in the bulk of the channel where the cells have uniform orientation θ = −θ_tilt/2 = −arctan(2τ/α)/2 and two boundary layers, whose size is roughly , with σ₀ = (α/2)sin 2θ₀ + τ cos 2θ₀, where the director interpolates between the bulk and boundary orientation (see Fig. 5b). This phenomenon closely resembles flow-alignment in nematics (see e.g.ref. 40) with −θ_tilt/2 playing the role of the Leslie angle θ_L = arccos(1/λ)/2. Whereas passive flow-alignment, however, requires λ > 1 (e.g. flow-aligning nematics), such an active flow-alignment occurs at any finite value of α and τ, provided the elastic boundary layer is sufficiently small to have a clear distinction between bulk and boundary alignment.

For small α and τ values, we can postulate that the nematic director will depart only slightly from its orientation at the boundary. Thus, taking again θ(y) = θ₀ + δθ(y) and linearizing eqn (21) around δθ = 0, we obtain:

∂_y²δθ + q²(δθ + δθ₀) = 0, (30)

with q² given, as before, by eqn (24) and:For a general anchoring angle θ₀, a solution of eqn (30) with boundary conditions δθ(0) = δθ(L) = 0 is given by:The associated velocity field can then be found from a direct integration of the linearized equation:The lack of a boundary condition for eqn (33) can be compensated with a global constraint on the total momentum, namely .

In the following, we provide explicit approximated expression for the velocity field in the special cases where θ₀ = 0 (parallel anchoring) and θ₀ = π/2. Furthermore, we will investigate the stability of the trivial solution of eqn (21) obtained when θ₀ is such that the chiral and non-chiral active stresses cancel each other identically.

1. Parallel anchoring. For θ₀ = 0 eqn (28) and (31) yield:The corresponding velocity field is then readily obtained by integrating eqn (33). This gives:where the wave number q_‖ is that given in eqn (28). Numerical solutions for this case are displayed in Fig. 5b for various α and τ values.

We stress that, whereas the flowing configurations resulting from the instability of the stationary state are left-right and clock-counterclockwise symmetric (i.e. the cells are equally likely to flow toward the negative or positive x-direction and, correspondingly, to tilt clock- or counterclockwise, see Section IV.A.1), in this case the direction of the tilt as well as that of the flowing monolayer is set by the signs of the constants α and τ.

2. Homeotropic anchoring. For θ₀ = π/2 from eqn (28) and (31) we find that the amplitude δθ₀ is given, once again, by eqn (34). Thus, the expressions for δθ and v_x are formally identical to those given in eqn (32) and (35), but with wave number q_⊥ as given in eqn (29).

3. Stationary solution. To conclude this subsection, we consider a special situation where the orientation of the cells at the boundary is fixed, as before, but such that the chiral and non-chiral stresses cancel each other identically. Thus: θ₀ = −θ_tilt/2 = −arctan(2τ/α)/2. In this case the orientation of the nematic director in the bulk of the channel, determined by the balance between the chiral and non-chiral active stress, is equal to that at the boundary. As a consequence, the boundary layer described in Section IV.B disappears and the system can achieve a stationary and uniformly aligned configuration. As those described in Section IV.A, however, the latter is unstable for sufficiently large active stresses or channel width.

Using the same algebraic manipulations adopted in Section IV.A, one can show that the perturbation δθ is again of the form given in eqn (25) with:

Analogously, the velocity is given by eqn (27), but with no constraint on the phase qL, because of the stress-free boundary conditions. As a consequence, the first mode to be excited is n = 1, thus the stationary state becomes unstable when q = q_c = π/L or, equivalently, when L = L_c = π/q. Some numerical solution of eqn (18), in this regime, are shown in Fig. 5c.

Notably, the transition from stationary to flowing is, in this case, no longer left-right and clock-counterclockwise symmetric, as in the examples discussed in Section IV.A, for any τ ≠ 0. This is well illustrated by the bifurcation diagram of Fig. 4b, showing the departure in the director orientation from the boundary value at the center of the channel [i.e. θ(1/2) − θ₀, with L = 1]. The dots have been obtained from a numerical integration of eqn (18), whereas the solid lines correspond to analytical solutions obtained by solving a third order equation for the constant C in eqn (25), as in Section IV.A.1. We find:

where and . For α > α_c, the solution consists of two branches, of which only one is connected with the stationary solution θ(y) = θ₀. Furthermore, the gap between the two branches increases monotonically with τ. If the instability is triggered upon applying a small random perturbation to the stationary state, this will always select the closest branch, thus the one connected to the trivial solution. As a consequence, a chiral cellular monolayer driven out of the stationary state by a small perturbation, will systematically tilt and flow in the same direction, which is in turn determined by the sign of the chiral active stress τ. The transition described above is known in bifurcation theory as a perturbed or imperfect pitchfork bifurcation and occurs when a standard pitchfork bifurcation, whose normal form is θ³ − μθ = 0, is biased by a small symmetry-breaking perturbation: i.e. θ³ − μθ + P_L + P_Sθ² = 0, where μ, P_L, and P_S are constant parameters. If P_L = P_S = 0, the equation is invariant under θ → −θ. Thus, for μ > 0, the trivial solution is unstable and the transition is supercritical, while for μ < 0, only the trivial solution is stable, and the bifurcation is subcritical. By contrast, for non-vanishing P_L and P_S, the equation is no longer invariant under θ → −θ and the bifurcation is no longer symmetric (see ref. 41 and 42 for an overview).

In our case, the role of the symmetry-breaking perturbation is played by the chiral stress τ. Thus, in the unperturbed scenario, τ = 0 and the stationary solution is θ = 0, with the critical α value being α_c = 2π²ηK/[γL²(1 − λ)]. When τ ≠ 0, on the other hand, expanding eqn (21) around θ = 0 and using ∂_y²θ = −q²θ one finds:

Evidently, this coincides with the normal form of a perturbed bifurcation for any finite τ value. For τ = 0, on the other hand, one recovers the normal form of the symmetric pitchfork bifurcation. Upon increasing τ, the bifurcation is shifted toward smaller α values, until, for τ = π²/ηK(γL²), α_c = 0 and the system is never stationary.

V. Discussion and conclusions

In this article we have investigated how a chiral and yet symmetric stress tensor might arise microscopically in nematic cell monolayers and how such a chiral stress influences some of the hallmark phenomena of active nematics. In Section II we proposed a microscopic model where a misalignment of the active force dipole and the cell's orientation is demonstrated to lead to a macroscopic chiral active stress tensor of the form eqn (2). In Section III we showed how the presence of chiral active stresses tilts the flow around ±1/2 disclinations, thereby leading to a misalignment between the defect polarity and the direction of motion, by an angle θ_tilt = arctan(2τ/α), with τ and α the chiral and non-chiral active stress respectively. In Section IV we investigated the spontaneous flow transition in a quasi-one-dimensional channel for both no-slip and stress-free boundary conditions as well as for various types of anchoring. For no-slip boundaries (Section IV.A), we recovered the classic pitchfork bifurcation first discussed by Voituriez et al.¹⁶ In this case the chirality does not affect the transition itself, but does leave a signature on the post-transitional configurations of the nematic director and velocity field, in the form of asymmetry with respect to the channel centerline. In case of stress-free boundaries (Section IV.B), we found that chirality renders the stationary and uniformly aligned configuration incompatible with most of the anchoring conditions. As a consequence, the cellular monolayer is always in motion, for any non-vanishing chiral and non-chiral active stress. For very large active stresses, in particular, we found an active analog of flow-alignment in nematics, with the bulk orientation (analogous to the Leslie angle⁴⁰) set by the ratio between chiral and non-chiral active stresses, i.e. −θ_tilt/2. Finally, in the special case in which the nematic director is anchored at an angle −θ_tilt/2 at the channel walls, we found that the spontaneous flow transition becomes asymmetric, i.e. only one of the two branches of the pitchfork bifurcation is connected to the trivial solution, which may then be the only one observed experimentally. This latter result could potentially explain the experimental observations by Duclos et al.,⁸ who found that NIH 3T3 cells are more likely to tilt clockwise then counterclockwise once the spontaneous flow transition sets up. In ref. 43 and 44 it was found that at the interface between an active nematic phase, with negligible distortions of the director field, and an isotropic phase there is an active anchoring angle that is set by the activity. This discussion can easily be repeated for the chiral active stress tensor in eqn (2) with the result being that the angle is changed due to the presence of chirality. Furthermore, considering a confined nematic phase without isotropic phase but with distorted director field, the case considered in Section IV being one example, it is readily found that the active anchoring angle is equal to the angle θ_tilt introduced above. Thus, the active anchoring angle considered in ref. 43 and 44 (no distortions of director field but region with active nematic and isotropic phase) and the anchoring angle θ_tilt derived above (distorted director field but only a nematic phase) can be seen as being the two special cases of the more general case with a nematic–isotropic phase interface and non-negligible distortions of the director field.

Further experimental investigations into the influence of chirality would be interesting. In particular, the tilt of the flow around ±1/2 disclinations has, according to our knowledge, not yet been observed. Thus, measurements of the tilt angle and experimental investigations of the flow field are needed to compare the theory with real-life cell monolayers. Additionally, as mentioned, the tilt angle opens a possibility to determine the relative magnitude of the chiral stress directly by particle-image-velocimetry measurements. Furthermore, since the cells used in ref. 8 were only weakly chiral the effects of chirality were not as pronounced. Performing similar experiments with cells with stronger chirality and for different boundary conditions would enable further tests of the presented theory.

Conflicts of interest

There are no conflicts to declare.

Appendix A: stress owing to a force quadrupole

The derivation of the active stresses given in Section II can be straightforwardly generalized to account for a more complex force distribution. For illustrative purposes we consider here the case of a quadrupole consisting of two force dipoles applied at the ends of the cell in longitudinal and transverse directions. In this case the force density field f_c is given by:

f_c = F^(a)_cδ(R_c − aν_c) − F^(a)_cδ(R_c + aν_c) + F^(b)_cδ(R_c − bν^⊥_c) − F^(b)_cδ(R_c + bν^⊥_c), (A1)

where a and b are, respectively, the major and minor semiaxis. Taking F⁽ⁱ⁾_c = F⁽ⁱ⁾_‖n + F⁽ⁱ⁾_⊥n^⊥, with i = a, b, and coarse-graining over the scale of a volume element, we find: P^a = −(aF^(a)_‖ + bF^(b)_⊥), α = 2(aF^(a)_‖ − bF^(b)_⊥) and τ = aF^(a)_⊥ + bF^(b)_‖.

Appendix B: nonlinear expansion

Expanding eqn (21) about θ₀ = 0 up to third order in δθ yields
which can be written as

The stress σ_xy is determined by the no-slip boundary condition and force balance, i.e.,

An expansion up to third order in δθ = C sin(2πy/L) around θ = 0 yields the conditionand thus σ_xy = τ − 2τC² sin²(2πy/L) = τ(1 − 2δθ²). This determines σ_xy and can be used in eqn (B2) to findComparing with the normal form given in Section IV.B.3 we find that, for L > L_c, μ = 3(λ − 1)/(2(4λ − 1)) is positive for all λ < 1/4 and λ > 1 and vanishes for λ = 1. With δθ = C sin 2πy/L solving for C at y = L/2 and expanding about L ∼ L_c to leading order yields the trivial solution C = 0 andwhich has a singularity at λ = 1/4. As remarked on in the general discussion in Section IV.B.3 the system displays a supercritical bifurcation for μ > 0 but a subcritical bifurcation for μ < 0. Thus the character of the bifurcation is determined by the sign of μ and in the present case this sign changes at the singularity λ = 1/4 and at λ = 1, where the system transitions between flow-tumbling (λ < 1) and flow-aligning (λ > 1).

Acknowledgements

We are grateful to Mattia Serra for illuminating discussions about Section IV. This work was supported by funds from the Netherlands Organisation for Scientific Research (NWO/OCW), as part of the Frontiers of Nanoscience program (L. G.), the Vidi (L. G. and L. A. H.), and Vici (R. M. H. M.) schemes and the Leiden/Huygens fellowship (K. S.).

References

1. G. K. Batchelor, J. Fluid Mech., 1970, 41, 545.
2. G. Duclos, S. Garcia, H. G. Yevick and P. Silberzan, Soft Matter, 2014, 10, 2346.
3. S. Garcia, E. Hannezo, J. Elgeti, J. F. Joanny, P. Silberzan and N. S. Gov, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 15314.
4. G. Duclos, C. Erlenkämper, J. F. Joanny and P. Silberzan, Nat. Phys., 2016, 13, 58.
5. 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.
6. K. Kawaguchi, R. Kageyama and M. Sano, Nature, 2017, 545, 327.
7. C. Blanch-Mercader, V. Yashunsky, S. Garcia, G. Duclos, L. Giomi and P. Silberzan, Phys. Rev. Lett., 2018, 120, 208101.
8. G. Duclos, C. Blanch-Mercader, V. Yashunsky, G. Salbreux, J.-F. Joanny, J. Prost and P. Silberzan, Nat. Phys., 2018, 14, 7.
9. A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans and F. Sagués, Nat. Commun., 2018, 9, 3246.
10. R. Kemkemer, D. Kling, D. Kaufmann and H. Gruler, Eur. Phys. J. E: Soft Matter Biol. Phys., 2000, 1, 215.
11. R. Kemkemer, V. Teichgräber, S. Schrank-Kaufmann, D. Kaufmann and H. Gruler, Eur. Phys. J. E: Soft Matter Biol. Phys., 2000, 3, 101.
12. L. Giomi, M. J. Bowick, X. Ma and M. C. Marchetti, Phys. Rev. Lett., 2013, 110, 228101.
13. L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek and M. C. Marchetti, Philos. Trans. R. Soc., A, 2014, 372, 20130365.
14. L. Giomi, Phys. Rev. X, 2015, 5, 031003.
15. S. Henkes, K. Kostanjevec, J. M. Collinson, R. Sknepnek and E. Bertin, 2019, arXiv:1901.04763.
16. R. Voituriez, J.-F. Joanny and J. Prost, Europhys. Lett., 2005, 70, 404.
17. S. Fürthauer, M. Strempel, S. W. Grill and F. Jülicher, Eur. Phys. J. E: Soft Matter Biol. Phys., 2012, 35, 89.
18. S. Fürthauer, M. Strempel, S. W. Grill and F. Jülicher, Phys. Rev. Lett., 2012, 110, 048103.
19. D. Banerjee, A. Souslov, A. G. Abanov and V. Vitelli, Nat. Commun., 2017, 8, 1573.
20. J. E. Avron, R. Seiler and P. G. Zograf, Phys. Rev. Lett., 1995, 75, 697.
21. T. J. Pedley and J. O. Kessler, Annu. Rev. Fluid Mech., 1992, 24, 313.
22. R. A. Simha and S. Ramaswamy, Phys. Rev. Lett., 2002, 89, 058101.
23. A. W. C. Lau and T. C. Lubensky, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 011917.
24. L. Q. Wan, K. Ronaldson, M. Park, G. Taylor, Y. Zhang, J. M. Gimble and G. Vunjak-Novakovic, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 12295.
25. Y. H. Tee, V. Thiagarajan, R. F. Hariadi, K. L. Anderson, C. Page, N. Volkmann, D. Hanein, S. Sivaramakrishnan, M. M. Kozlov and A. D. Bershadsky, Nat. Cell Biol., 2015, 17, 445.
26. K. Schakenraad, J. Ernst, W. Pomp, E. H. J. Danen, R. M. H. Merks, T. Schmidt and L. Giomi, 2019, arXiv:1905.09805.
27. K. E. Worley, D. Shieh and L. Q. Wan, Integr. Biol., 2015, 7, 580.
28. T. Sanchez, D. N. Chen, S. J. DeCamp, M. Heymann and Z. Dogic, Nature, 2012, 491, 431.
29. F. C. Keber, E. Loiseau, T. Sanchez, S. J. DeCamp, L. Giomi, M. J. Bowick, M. C. Marchetti, Z. Dogic and A. R. Bausch, Science, 2014, 345, 1135.
30. P. Guillamat, J. Ignés-Mullol and F. Sagués, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 5498.
31. P. Guillamat, J. Ignés-Mullol and F. Sagués, Nat. Commun., 2017, 8, 564.
32. R. Zhang, N. Kumar, J. L. Ross, M. L. Gardel and J. J. de Pablo, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, E124.
33. A. J. Vromans and L. Giomi, Soft Matter, 2016, 12, 6490.
34. A. Maitra and M. Lenz, Nat. Commun., 2019, 10, 920.
35. D. Marenduzzo, E. Orlandini, M. E. Cates and J. M. Yeomans, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 031921.
36. L. Giomi, M. C. Marchetti and T. B. Liverpool, Phys. Rev. Lett., 2008, 101, 198101.
37. S. A. Edwards and J. M. Yeomans, Europhys. Lett., 2009, 85, 18008.
38. J. Hardoüin, R. Hughes, A. Doostmohammadi, J. Laurent, T. Lopez-Leon, J. M. Yeomans, J. Ignés-Mullol and F. Sagués, 2019, arXiv:1903.01787.
39. V. Fréedericksz and A. Repiewa, Z. Phys., 1927, 42, 532.
40. P. G. de Gennes and J. Prost, The Physics of Liquid Crystals, Oxford University Press, Oxford, 1993.
41. M. Golubitsky and D. G. Schaeffer, Singularities and Groups in Bifurcation Theory, Springer-Verlag, New York, 1985, vol. I.
42. P. M. Kitanov, W. F. Langford and A. R. Willms, Dynam. Cont. Dis. Ser. A, 2013, 20, 197.
43. M. L. Blow, S. P. Thampi and J. M. Yeomans, Phys. Rev. Lett., 2014, 113, 248303.
44. A. Doostmohammadi, S. P. Thampi and J. M. Yeomans, Phys. Rev. Lett., 2016, 117, 048102.

---