Stress-shape misalignment in confluent cell layers

This study investigates the relationship between cell shape and cell-generated stresses in confluent cell layers. Using simultaneous measurements of cell shape orientation and cell-generated contractile forces in MDCK and LP-9 colonies, we report the emergence of correlated, dynamic domains in which misalignment between the directors defined by cell shape and by contractile forces reaches up to 90$^o$, effectively creating extensile domains in a monolayer of contractile cells. To understand this misalignment, we develop a continuum model that decouples the orientation of cell-generated active forces from the orientation of the cell shapes. This challenges the prevailing understanding that cells throughout a tissue create either contractile or extensile forces, and the validity of the usual active nematic models of cell motility where active forces are strictly slaved to cell shape orientation.


Main
Cells are the fundamental building blocks of life, and their ability to collectively generate active forces plays a crucial role in physiological processes from morphogenesis [1], tissue growth and repair [2] to apoptosis [3], tumour development [4] and metastasis [5].Confluent cell monolayers plated on adhesive substrates are widely used as model systems in investigations aimed at understanding collective cell motility [6].There is growing evidence that the dynamics of such confluent cell layers can often be well described by the theories of active nematics [7][8][9][10][11].
Describing cells as active emphasises that they continuously take energy from their surroundings and use it to initiate life processes [12,13].Nematic particles are elongated in shape, and nematic ordering occurs when their long axes tend to align parallel (Fig. 1a).Such ordering has frequently been observed in long, thin cells such as fibroblasts [7,14] or LP-9 [15] but is more surprising in cell types that are on average isotropic, such as the Madin-Darby Canine Kidney (MDCK) cell line.Here extensions in cell shape, driven by active forces, are locally correlated to give nematic order [16][17][18].
The combination of activity and nematic ordering leads to striking collective behaviours which are mirrored between active nematic models and cell monolayers.These include active turbulence, characterised by cell velocities that are chaotic with regions of high vorticity [19], spontaneous directed flow in confinement [20], and the identification of motile topological defects [17], long-lived cell configurations at which domains of different cell orientations meet and the nematic order is frustrated.Fig. 1b,c show the cell orientations around the +1/2 and -1/2 defects which predominate in 2D cell layers.
The direction of the long axis of a cell is an obvious way to define a preferred nematic shape axis.The local axis of principal stress gives an alternative way to choose the local axis of any nematic ordering.In almost all work to date it has been assumed that the two definitions are equivalent, meaning that the stress and shape axes are tightly coupled [16,[21][22][23].Under this assumption, differences between the axes of stress and shape would be modest and attributable to biological noise.Here we challenge this assumption and show that there are dynamic, correlated regions in cell layers, where the stress and shape axes are systematically misaligned.Introducing the possibility of such misalignment in a continuum model of active nematics reproduces the temporal and spatial correlations and misalignment angle observed in our experiments and emphasises the key role of active flows in driving the misalignment.
Our experiments were performed on confluent MDCK layers of diameter 1 mm plated on polyacrylamide substrates with a Young's modulus of 6 kPa (Fig. 1d).We define the cell shape orientation in the tissue by assigning each cell a director (headless vector) n which lies along the long axis of the anisotropic cell shape as shown in Fig. 1e [24].Monolayer Stress Microscopy was used to measure the stress tensor σ, from which we computed the orientation of the first principal stress, which defines the local orientation m along which contractile forces are generated (Fig. 1f).By interpolating between individual cells we obtain continuous director fields n and m which respectively describe cell shape orientation and the principal axis of contractile stress throughout the tissue (see Methods).
We analysed the local cell shape and stress orientations over the course of 12 hours in 11 MDCK tissue samples, in time lapse experiments taking data every 15 minutes.We define θ as the misalignment angle between the local cell shape orientation n and the principal axis of contractile stress m in the tissue (Fig. 2a).The distribution of θ is shown in Fig. 2b.While most cells create contractile stresses along their cell shape axis (θ ≈ 0), there is a large number of cells in which the axis of contractile stress is significantly misaligned with respect to shape orientation.If the misalignment angle reaches θ ≈ 90 • , cells create contractile stresses perpendicular to the cell shape orientation, thereby pulling inward not along their long shape axis but rather along their short shape axis.In the following we will refer to cells with large misalignment (θ > 45 • ) as extensile and cells with small misalignment (θ < 45 • ) as contractile following the usual terminology in the mathematics and active matter literature, and we distinguish these ranges of θ as blue and red in Fig. 2b.
We now investigate the spatial and temporal correlation of the shape and the stress orientations in the MDCK monolayers.Fig. 2c (left) shows a tissue snapshot where the cell shape orientation field n is shown as black lines and the color map again indicates whether θ is greater (extensile, blue) or less (contractile, red) than 45 • .Fig. 2c (right) shows the time evolution of a region of the tissue with snapshots taken at 15-minute intervals.Similar data is presented dynamically in Movie 1.It is evident that the misalignment angle forms evolving spatio-temporal patterns where extensile cells form small, dynamic clusters in a mostly contractile background.The extensile clusters grow, shrink and coalesce over time.The time-averaged area fraction of extensile cells is 27 ± 4% (Fig. 2d).
To further quantify the spatial patterns we calculated the spatial and time correlation functions of the cell shape orientation, the cell stress orientation, and the mismatch angle θ.These are defined as where ψ x represents the shape director angle, stress angle, or the mismatch angle θ, and ⟨. ..⟩ r0,t0 denotes an average over space (a circle of diameter 312µm in the centre of the island to avoid edge effects) and time.The spatial correlations, shown in Fig. 3a, indicate a length-scale ∼ 50µm for the cell stress orientation, and a longer length scale ∼ 100µm for the cell shape orientation.From the time correlation functions, in Fig. 3b, we identify a time-scale for the decay of the extensile patches ∼ 300 minutes.
Many cells contain bundles of actomyosin, termed stress fibres, that tend to form along the long axes of cells and are the primary source of contractile stresses [25].We hypothesise that the regions of large misalignment angle are due to active flows disturbing the natural alignment of the stress fibres with the long axis of the cell due to different responses of the shape axis n and the stress axes m to flows.This leads to the formation of extensile regions which have a large mismatch between the shape and stress axes.A cell's stress axis and shape axis then gradually relax towards each other.
Continuum tissue models, based on the equation of motion of active nematics, have been very successful in explaining cell motility on a coarse-grained level [7][8][9][10]17].However, the assumption has always been that the principal axis of contractile stress m and the shape axis n are indistinguishable.Therefore, to investigate the consequences of our hypothesis, we extend the continuum modelling to decouple m and n.
We describe the shape of the cells and the stress using nematic order parameters [24].The nematic order parameters encode the magnitude of nematic order in cell shape S n or in the stress S m , and the director field associated with cell shape, n or the stress m.We assume that the flows are created by contractile active forces that act along the direction of stress fibres m.The active flows advect the cells, and the vorticity of the flows rotates both the shape and the stress director fields.The experimental spatial correlation functions show that the nematic order of the shape director n has a longer length scale than that of the stress director m which we model by choosing different elastic constants in the free energy.This also means that the shape and the stress directors respond in different ways to the vortical flows leading to misalignment between m and n.We include a term in the free energy which acts to slowly relax the misalignment.In the continuum simulations, we have the freedom to choose a length and time scale, and we do this by matching the length and time scales of the correlation functions (Eq. 1) between simulation and experiment, as shown in Figs.3a and 3b.See Methods for further details of the modelling.
We compare the simulation results to the experiments in Fig. 2. In agreement with the experiments, spatially correlated domains of extensile cells in a contractile background emerge in the simulations (Fig. 2f and Movie 2).We also obtain a quantitative match to the probability distribution for the misalignment angle θ if we use a realignment time scale of 25 minutes (Fig. 2b).
Active turbulence is also characterised by topological defects and we noticed that topological defects in the cell orientation field tend to sit at the interfaces between extensile and contractile domains.This is illustrated in Figs.2c (experiment) and 2f (simulations) where +1/2 defects are indicated by a yellow arrow and −1/2 defects by a green trefoil.To quantify the results, we consider the points that are at a distance smaller than r < 5.2µm (equal to the grid size after interpolation) from the interface as the interface domain while the rest of the tissue is defined as the bulk domain.A spatial and temporal average shows that interface domains make up about 18 ± 2% of the tissue area fraction in the experiments, and about 15 ± 3% of the area fraction in simulations.However, about 86 ± 2% (86 ± 5%) of defects lie within interface domains in experiments (simulations).These results are shown in Fig. 2e.
Our interpretation leads to the prediction that in the extensile regions the stress fibres will be in the process of re-forming to realign with the new direction of cell elongation, and are therefore less efficient in producing contraction.In Fig. 3c we plot the magnitude of the contractile stress as a function of θ showing a clear decrease.Moreover, to check our interpretation visually we fixed cells for fluorescent imaging of actin fibres in a selection of samples.We found that it was rarely possible to visually ascertain an unambiguous direction of the stress fibres in cells with large θ whereas stress fibres were much clearer in cases when they were aligned along the long axis of the cell (see Supplementary Material).
We next, as a comparison, performed similar experiments on the human mesothelial cell line LP-9 (Fig. 4a).These cells, which have an elongated morphology with a high aspect ratio, showed a behaviour which contrasted with the MDCK islands.A very small number of topological defects were present at the beginning of the experiments.These persisted and remained in approximately the same position throughout the experiments (40 hours) and no new defects were created, indicating that the cell layer was behaving primarily as a passive nematic, with any active flows not sufficiently strong to create defects or substantially change the cell orientation.
There was, however, still a population of extensile cells with θ > 45 • , but it was far smaller than in the MDCK monolayers.The extensile cells formed small (3.4% of the area of the tissue) clusters (Fig. 4b) adjacent to the defects in the cell shape director field.We found no evidence that the extensile regions disappear within the time scale of the experiment.
As shape director n does not change we modelled the LP-9 cells by fixing a defect in cell shape n at the center of a circular cellular island and allowed the stress field m to relax.An extensile region was indeed formed next to the defect in cell shape n, as shown in Fig. 4c.This dynamical steady state is a result of the balance between the elastic energy, which favours nematic alignment of the stress directors m, and the term which encourages m to align with n.The exact position and size of the extensile region relative to the defect varies, depending on the initial condition for m, indicating the existence of metastable solutions.
We have shown that active flows and the different elastic properties of cell shapes and active stress filaments lead to long-lived regions of large shape-stress misalignment within cell monolayers.Experimental results on MDCK and LP-9 cells are interpreted in terms of an active nematic model which corrects the usual assumption that the shape and stress axes are identical, but rather allows them to respond differently to flows causing misalignment between the shape and stress axes.In MDCK monolayers the shape and stress slowly realign, thus introducing a time-scale that must be taken into account when describing the monolayer dynamics.Our results question the common practice of using the direction of defect velocities to identify cell monolayers as contractile or extensile.
Time lapse imaging and analysis: Polyacrylamide gels embedded with fluorescent particles (580/605, diameter 0.5 µm, carboxylate modified; Life Technologies) were fabricated with Young's moduli of 6 kPa and thickness of 150 µm using methods described previously [26,27].Polydimethysiloxane (PDMS, Sylgard 184) was cured in 200 µm thick sheets.The sheets were cut into 20×10 mm squares and then 1 mm holes were cut using a 1 mm biopsy punch.The PDMS masks were adhered to the gels using previous methods [26], and the 1 mm circular hole was coated with 0.01 mg/ml type I rat tail collagen I (BD Biosciences) with the covalent crosslinker sulfo-SANPAH (Pierce Biotechnology).MDCK and LP-9 cells were seeded onto 1 mm islands 24 hr before imaging.The cells and particles were imaged using an Eclipse Ti-E microscope (Nikon, Melville, NY) with a 10× numerical aperture 0.5 objective (Nikon) and an Orca Flash 4.0 digital camera (Hamamatsu, Bridgewater, NJ) running on Elements Ar software (Nikon).All imaging was done at 37 • C and 5% CO 2 .The cells were imaged every 15 or 20 min for 24 hr.After imaging, the cells were removed by incubating in 0.05% trypsin for 1 hr, and images of the fluorescent particles in the substrate were captured for a traction-free reference state for traction force microscopy [28].To this end, Fast Iterative Digital Image Correlation was used [29] followed by Fourier transform traction microscopy [30] accounting for finite substrate thickness [31,32].The displacements were computed using 32×32 pixel subsets (21×21 µm 2 ) with a spacing of 8 pixels (5.2 µm).Stresses within the monolayer were computed with monolayer stress microscopy [33,34].Cell orientations were determined using the ImageJ plugin OrientationJ.
Imaging and analysis for fixed cells: Imaging actin stress fibers required fixing the cell monolayers, which had to be done after collecting data for traction and stress measurements [36,37].To this end, reference images of the √ N , where L is system size and N is the number of cells in the snapshot.c) Linear interpolation is used to find the director field on a lattice with a smaller mesh.This step is required to find the position and orientation of defects correctly.d) Using the director on the new lattice, we then use the defect-identifying algorithm introduced in Ref. [35] to find the direction and orientation of ±1/2 defects.
fluorescent particles were collected before seeding the cells.Polyacrylamide gels were made as described above, and reference images of the fluorescent particles in the reference (undeformed) state were collected.Microscopy was the same as described above with the exception of using a 20× numerical aperture 0.5 objective (Nikon).MDCK cells were seeded and allowed to come to confluence overnight.The cell culture medium was changed 1 hr prior to imaging, and cells were imaged every 10 min for 1 hr.Cell monolayers were fixed using chilled 4% paraformaldehyde solution for 20 min.Cells were then stained for actin using ActinRed 555 ReadyProbe Reagent (Invitrogen, catalog number R37112) according to manufacturer instructions and images.To analyse the orientation of the fixed cells, the cells were manually segmented based on images of the cell cortex.The orientation of actin stress fibers was identified manually for cells having visually clear stress fibers.Stresses were computed using monolayer stress microscopy, as described above, but for these experiments, the stresses at the boundaries were unknown, meaning that the recovered stresses represented not the full stress tensor but rather deviations from the average of the stress tensor at the boundaries.Our prior experiments with this cell type show that the average stress tensor is nearly isotropic (indicating small shear stresses) [27,38].If the stress state at the boundary is isotropic, then there is no error in computing the orientation of principal stresses and the traceless stress tensor, both of which are reported in the main text for experiments using fixed imaging.
Construction of the director fields for cell orientation and contractile stress: Fig. 5 summarise the process of defining smooth director fields and identifying the position and orientation of defects.First, cell positions and orientations are mapped onto a square lattice with a lattice constant dl = L/ √ N , where L and N are the system length and number of the cells in the snapshot, respectively.To find the orientation of the cells on a lattice, each cell is mapped to the closest lattice site.The result is shown in Fig. 5(b).In the original lattice built by the cells, the lattice spacing is large and that makes it impossible to find the position and orientation of the defects accurately.As a result, we need to find the cell orientation on a lattice with a smaller mesh size.We construct this lattice by linear interpolation between lattice points (Fig. 5c).Using the interpolated director field, we can then use a defect finding algorithm to find the position and orientation of ±1/2 defects (Fig. 5d).
The stress matrix measured in experiments has three independent components σ xx , σ xy and σ yy .It has a non-zero trace, and we first make it traceless by adding a constant c to the diagonal elements, so that σ yy + σ xx + 2c = 0. To find the orientation of the contractile stress, we find the two mutually perpendicular axes which are parallel to the orientations of positive and negative principal stress.These axes can be found by a rotation of the stress matrix through an angle θ p such that the shear stress σ xy becomes zero.We note that there are two directions over which the shear stress becomes zero, θ p and θ p + π/2.We define the orientation of contractile cell-generated stress to be along the orientation of the positive stress (pointing outwards).

The continuum model for active cell monolayers:
We describe the motion of cells within the monolayer by a continuous velocity field u.The local orientation of cell shapes is described by a tensor field Q n = S n (nn − I/2) and the orientation of the stress fibres by a second tensor field Q m = S m (mm−I/2).These nematic order parameters encode the magnitude of nematic order in the cell shape, S n , or the magnitude of the contractile stress, S m , and the director field associated with the cell shape, n, or the stress, m [24].This description differs from previous active nematic continuum theories in that it allows for a finite misalignment angle θ = cos −1 (n • m) between the elongation axis of cell shape and the axis along which contractile forces are generated by stress fibres (see Fig. 2a).
Following empirical arguments, we assume that in equilibrium the cell monolayer is governed by the following free energy density: In the absence of activity, the first and the second terms lead to a phase with nematic order in cell shape and cell stress.The third and the fourth terms penalise gradients in the shape and stress orientations, respectively.These are motivated by the observations of nematic ordering in both n and m in the experiments.The final term tends to align shape and stress orientations.The shape orientation, n, and the stress orientation, m, change in response to active flows.Since they have different elastic constants, they respond differently to the active flows which can lead to a mismatch between their directions.The dynamics of the nematic tensors is governed by [39] where γ is the rotational diffusivity, is the fluid vorticity, ∆ m/n is a uniform noise in the rotation, the molecular field H x = −( δf δQ x − I 2 Tr δf δQ x ) shows the relaxation of the orientational order to the minimum of the free energy, and we have set the flow tumbling parameter equal to zero.
We assume that the flows observed in confluent cell monolayers can be well approximated by where ρ is mass density and the stress tensor Π = Π passive + Π active includes a passive and an active contribution, where the passive, viscous terms Π passive are well known from liquid crystal hydrodynamics [39,40].Flows in confluent cell layers are predominantly driven by active dipolar forces created on the single-cell level by stress fibers which convert chemical energy into mechanical work.This gives an active term in the stress tensor, where we choose ζ < 0 to correspond to contractile forces.The equations are solved using a hybrid lattice Boltzmann algorithm [40].The MDCK simulations are performed in a 200×200 box with periodic boundary conditions over 120000 lattice-Boltzmann time-steps, and data is collected every 300 time-steps.The measurements are performed in steady state when the mean number of defects and the fraction of the extensile area do not, change over time.The initial orientation of both n and m is random, and the magnitude of the order is S n = S m = 1.We choose values of parameters that lead to an active fluid in a low Reynolds number regime: We set up the LP-9 simulations such that nematic order only forms inside a circular region with radius R = 80.The free energy in this region is again given by Eq. ( 2), but the bulk free energy outside the circle is which sets the magnitude of the order to be zero.We impose a defect in the shape director n at the centre of the inner region by setting the director angle with the x-axis to be equal to ϕ/2, where ϕ is the polar angle in the co-ordinate system centered at the defect core.We do not allow n to vary in time; m relaxes towards the minimum of the free energy.We use the same parameter values as for the MDCK cells except K m = 0.02, K n = 0.01, ζ = 0, ∆ m = ∆ n = 0, C ′ = 0.003.We note that the quantitative fits achieved using this model must be viewed with some caution as there are several adjustable parameters.However, the qualitative behaviour is insensitive to the numerical values of the parameters.

Movie captions:
Movie 1: Tissue dynamics with the cell orientation field n shown as black lines on top of a colour map distinguishing contractile (red) and extensile (blue) regions.Topological defects in the cell orientation are indicated by yellow (+1/2) and green (−1/2) symbols.

FIG. 1 .
FIG. 1.(a) Nematic ordering of cell shape in a confluent cell monolayer.Red lines indicate the nematic directors n, which lie along the long axis of a cell, and which tend to align parallel.(b) Cell shape orientations around a +1/2 topological defect shown in yellow.(c) Cell shape orientations around a −1/2 topological defect shown in green.(d) Geometry of the MDCK monolayers used in the experiments.(e) Example of a typical experiment measuring the long axis of cell shape orientation, n, (red lines) in an MDCK monolayer.(f) A typical set of results from experiments measuring the orientation of maximum contractile stress m (white lines) overlaid on the orientations of cell shape, n (red lines).In (e)-(f) the lengths of the red (white) lines are proportional to the deviation of the cell aspect ratio from unity (magnitude of the contractile stress).

FIG. 2 .
FIG. 2. (a) Definition of the misalignment angle θ between the shape orientation axis n, and the principal axis of contractile stress m.(b) Distribution of the misalignment angle θ.Red/blue colouring denotes contractile (θ < 45 o )/extensile (θ > 45 o ) values.The distribution contains data from 11 independent MDCK islands and all time points.The black dotted line shows the result from simulation.(c) left: Tissue snapshot with the cell orientation field n shown as black lines on top of a colour map distinguishing contractile (red) and extensile (blue) regions.Topological defects in the cell orientation are indicated by yellow (+1/2) and green (−1/2) symbols.Right: Snapshots of a region of the same tissue taken 15 minutes apart showing the evolution of extensile clusters.The time axis is from left to right.(d) Experimental time average of the area fraction of contractile (red) and extensile (blue) cells for 11 different cell islands.The crosses show the results from simulations.(e) Defects are preferentially found in the vicinity of boundaries between extensile and contractile regions.(f) Snapshot from simulations with the cell orientation field n shown as black lines on top of a colour map distinguishing contractile (red) and extensile (blue) regions.Topological defects in the cell orientation are indicated by yellow arrow (+1/2) and green trefoil (−1/2) symbols.The time axis is from left to right.

FIG. 3 .
FIG. 3. (a) Decay of the spatial correlation functions, C x (r) = ⟨cos 2[ψx(r + r 0 , t 0 ) − ψx(r 0 , t 0 )]⟩t 0 ,r 0 , where ψx represents shape director angle C n , stress angle C m , or the mismatch angle C θ .Comparison of the simulation and the experiment assumes 10LB spatial units ∼ 44µm.(b) Decay of the time correlation functions, C x (t) = ⟨cos 2[ψx(r 0 , t + t 0 ) − ψx(r 0 , t 0 )]⟩t 0 ,r 0 , where ψ represents shape director angle C n or the stress angle C m .Comparison of the simulation and the experiment assumes 100LB time units ∼ 10 min.(c) Magnitude of the stress σ as a function of the misalignment angle θ.The magnitude of the stress is re-scaled with its maximum value in each experiment before averaging.In all plots the error bars show the standard deviation over 11 experiments.

FIG. 4 .
FIG. 4. (a) Snapshot of an LP-9 island.(b) left: Tissue snapshot with the cell shape orientation field n shown as black lines on top of a colour map distinguishing contractile (red) and extensile (blue) regions.Topological defects in the cell shape orientation are indicated by yellow (+1/2) and green (−1/2) symbols.Middle: Snapshots of a region of the same tissue taken 40 minutes apart showing that the extensile cluster remains fixed.The time axis is from top to bottom.Right: Close-up view of the same tissue.Cell shape orientation is shown in red and stress orientation in white.The colormap shows the extensile (blue) and contractile (orange) regions (smoothed).The defect in the cell shape (stress) orientation is shown in yellow (magenta).(c) Distribution of the misalignment angle θ.Red/blue colouring denotes contractile (θ < 45 o )/extensile (θ > 45 o ) values.The LP-9 cells predominantly form contractile regions (compare with Fig. 2b for MDCK cells).(d) Left: Snapshot from simulations with the cell shape orientation n shown as black lines on top of a colour map distinguishing contractile (red) and extensile (blue) regions.The +1/2 defect in the cell shape is indicated by a yellow arrow.Right: Close-up view of the same tissue.Cell shape orientation is shown in red and stress orientation is white.The colormap shows the extensible (blue) and contractile (orange) regions (smoothed).The defect in the cell shape (stress) orientation is shown in yellow (magenta).

FIG. 5 .
FIG. 5. Construction of the director fields for cell orientation and contractile stress: a) Cell position and orientation, as shown in Fig. 1. b) Cell position and orientation after mapping to a lattice.The lattice unit is equal to dl = L/√ N , where L is system size and N is the number of cells in the snapshot.c) Linear interpolation is used to find the director field on a lattice with a smaller mesh.This step is required to find the position and orientation of defects correctly.d) Using the director on the new lattice, we then use the defect-identifying algorithm introduced in Ref.[35] to find the direction and orientation of ±1/2 defects.

Movie 2 :
Dynamics in the simulations, with the cell orientation field n shown as black lines on top of a colour map distinguishing contractile (red) and extensile (blue) regions.Topological defects in the cell orientation are indicated by yellow (+1/2) and green (−1/2) symbols.Contributions: MRN, LR, GZ and JMY formulated the model and interpreted the results.MRN, MM, JZ, and JN analysed the experimental data.MRN performed the simulations.JN designed and supervised the experiments.MM and JZ performed the experiments.LR, MRN, and JMY drafted the manuscript.All authors commented on the manuscript.