Probing epithelial tissue rheology with a Stokes experiment

We perform a bidimensional Stokes experiment in an active cellular material: we force a monolayer of Madin-Darby Canine Kidney (MDCK) epithelial cells to migrate around a circular obstacle within a long and narrow channel, involving cell shape changes and neighbour rearrangements. Based on image analysis of tissue flow velocity and cell shape, we determine the tissue strain rate, cell deformation and cell rearrangement rate fields, which are spatially heterogeneous. We find that the cell deformation and cell rearrangement rate fields correlate strongly, which is compatible with a viscoelastic liquid behaviour with a relaxation time $\tau = 70\pm 15$ min. In our experiment, the monolayer behaves as a flowing material with a Weissenberg number close to one which shows that both elastic and viscous effects can have comparable contributions in the process of collective cell migration.

Introduction Epithelial tissues are active cellular materials made of constitutive objects, the cells, that can not only deform and exchange neighbors, but also grow, divide, have a polarity and exert active stresses [1][2][3][4]. Tissue mechanical properties have a crucial importance in biological morphogenetic processes [5]. In particular, cell neighbour rearrangements contribute to shaping tissues, as during the Drosophila embryo germ band extension [6]. Rearrangements are often described as either passive, caused by an external stress emerging at the tissue scale; or active, for instance triggered and directed by an anisotropic distribution of molecules (such as myosin or cadherins) at the cell-cell contacts [7]. Active cell contour fluctuations combined to stress emerging at the tissue scale are key ingredients to fluidize the tissue both in vitro [8] and in vivo [9].
Experiments performed on embryonic tissues [10], multicellular spheroids [11,12], or cell monolayers with or without substrate [13,14] have suggested describing tissues as a viscoelastic liquid. However, there is still a debate around the microscopic origin and value of the viscoelastic relaxation time τ and on whether MDCK monolayers behave predominantly as liquids or solids [13,15]. MDCK cells can actively migrate on a substrate while sustaining tissue cohesivity. Collective cell migration is due to the cells displaying simultaneously cryptic lamellipodia on their basal side and adherent junctions on their apical side [16,17]. Studying their bidimensional flow in long and narrow adhesive strips facilitates experiments, simulations, theory and their mutual comparisons [14,[18][19][20][21].
The flow of a liquid relative to a circular obstacle, whether in two or three dimensions, is a classical experiment introduced by Stokes [22] to measure a viscosity, and/or to probe the mechanical behaviour of viscoelastic or viscoplastic materials [23]. More complex materials which are together viscous, elastic and plastic such as liquid foams (bubbles are able to sustain elastic deformations due to their surface tension, until a yield point where they rearrange and exchange neighbours [24,25]) have been well characterised using the Stokes experiment [26]. The Stokes flow geometry has even been applied in contexts as different as granular materials [27], or active materials : to study clogging in animal crowds [28], or cell-substrate interaction in MDCK cell monolayers, where interaction forces were found directed towards the obstacle [29].
The total deformation rate is the strain rate, also called velocity gradient,ε tot . It is the sum of the cell shape deformation rate, which is the time derivative of the elastic cell shape deformation ε e , and of the intercellular topological change rate, which in absence of cell divisions and apoptoses reduces to the cell cell rearrangement rateε r , so thaṫ ε tot =ε e +ε r [30]. In the Stokes flow, these deformation fields are heterogeneous: it displays simultaneously a variety of combinations of these fields, depending on the position with respect to the obstacle; for instance there are places where the strain rate and shape deformation are parallel, other places where they are orthogonal. In that respect, to discriminate between different models, the Stokes flow geometry is more efficient [26] than a homogeneous flow such as that between parallel plates (Couette flow) [31].
Here, we observe over days the Stokes flow of MDCK cell monolayers (Fig. 1a,b) involving mostly cell rearrangements (Fig. 1c,d) and only few cell divisions. By using image analysis methods which do not require cell contour segmentation, we perform quantitative kinematic measurements within the frame of continuum mechanics: velocity, velocity gradient, and cell shape fields (Fig. 2). We deduce the cell rearrangement rate field and quantify the fields correlations to probe the tissue rheology (Fig. 3). Methods The procedure for micropattern printing, cell culture, imaging and velocity measurement is described in details in Ref. [15]. Briefly, a strip is 4 mm long; it is adhesive, while its four boundaries and the circular obstacle are not. To check reproducibility, each experimental batch is composed of two identical strips. To test the effect of experiment dimensions, three batches are used: a first one with obstacle diameter 150 µm and strip width 750 µm; a second one with obstacle diameter 200 µm and strip width 1000 µm; a third one with obstacle diameter 300 µm and strip width 1000 µm.
We inhibit cell divisions using mitomycin; this prevents cell density increase and jamming that slows down migration. MDCK monolayers can then migrate during days (a timescale much longer than the cell cycle duration) and over millimeters (a lengthscale much larger than the cell size). The 2 mm long region upstream of the obstacle serves as a cell reservoir; its density is initially high and decreases, at constant cell number. The typical time for cells to migrate over a distance equivalent to the obstacle diameter of 200 µm is 3 h. For obstacle diameters of 400 µm or larger, the transit time of a cell around the obstacle would be too large to establish a steady flow. For obstacle diameters of 100 µm or smaller, cells would make suspended bridges over the non-adhesive obstacle.
We divide the field of view in square boxes of 128 pixel side, containing typically 10 cells. We measure the twodimensional velocity field v(x, y, t) using a custom-made Matlab optic flow code based on the Kanade Lucas Tomasi (KLT) algorithm [32] with a level 2 pyramid. This coarse-grained velocity field can be averaged in time (see Fig.  2a), yielding components v i (x j ) where i, j = 1 or 2. Using finite differences we obtain in each box the symmetric part of the velocity gradient ∇ v sym , with components (∂ i v j + ∂ j v i )/2. This symmetric tensor can be diagonalized; its deviator part, corresponding to the anisotropic component of the tensor, is diagonalized too. We graphically represent in each box (Fig. 2b) its anisotropic part, with bars oriented along the direction of the eigenvector corresponding to the positive eigenvalue (local tissue expansion direction) and of length proportional to the eigenvalue.
To extract the coarse-grained cell anisotropy, we use Fourier Transform (FT, see Fig. S1); details and validations are provided in the companion paper [33]. This method provides an efficient measurement of the cell deformation field without having to recognize and segment each individual cell contour which is challenging on phase contrast images.
We use the same grid and box size as for the velocity. In each box, we multiply the image by a windowing function to avoid singularities in the FT due to box boundaries. We then compute the FT using the Fast Fourier Transform algorithm implemented in Matlab (fft2.m) and keep only its amplitude (not the phase), which can be time averaged to increase signal to noise-ratio.
We binarize the resulting Fourier space pattern keeping 5% of the brightest pixels (for justification of this percentile value, see the companion paper [33]). We compute the inertia matrix of the binarized FT pattern and diagonalize it, which yields two eigenvalues λ 2 max and λ 2 min in the directions of the pattern main axes. They determine the ellipse back in the real space, with eigenvalues L max = 2m λmin and L min = 2m λmax in the directions of the same axes with m the size of the FT image in pixels.
We define as follows the cell shape deformation tensor ε e with respect to a rest state which we assume to be isotropic, with two equal eigenvalues L 0 = √ L max L min : ε e has the same eigenvectors as the Fourier pattern, and two eigenvalues 1 2 ( . The deformation tensor ε e is equivalent to other definitions of the strain (e.g. true strain [33,34]) within a linear approximation; in addition, ε e has the advantage to have well established transport equations [35].
We represent the anisotropic part of ε e , namely its deviator ε dev e = ε e −Tr(ε e )/2, which only depends on the ratio Lmax Lmin . This represents the average anisotropy depicted as a bar in Fig. 2c, and Fig. S1. We deduce the time averaged cell rearrangement rate ε r by measuring the difference ε tot −ε e as in [30] as divisions are inhibited. From the measurement of ε e , we estimateε e by taking into account the cell deformation advection in the flow, aṡ ε e = Dεe Dt = ∂εe ∂t + v · ∇ε e : this significantly improves the results presented below (we still neglect rotation terms [34,35], which do not change the present results, see Fig. S2). We correlate the deviator of the cell deformation ε e with the deviator of the rearrangement rate ε r = ∇v sym − ∂εe ∂t − v · ∇ε e where both fields are time-averaged over the experiment duration (at least 10 hours).
Results The anisotropic parts of cell deformation field and cell rearrangement rate field display a strong linear correlation (Fig. 3), correlation value R = 0.67. The orientations of these fields are identical, excepted in a region where both fields are small (bottom right of Fig. 3c, coded in blue). Their proportionality coefficient is τ = 70±15 min: this value is independent on the experiment dimension and results from the average of six measurements (Fig. 3a, Supp. Fig. 3), which result from fits which correlations values range from R = 0.53 to 0.77. Note that the numerous data corresponding to small cell deformation (deviator of ε e with an amplitude smaller than 0.05) are not included in the fits as they would add noise but no information to the fit. Each individual experiment provides a self-sufficient data set, which compares well with the line of slope τ (Fig.  3b, Fig. S2) thanks to the Stokes flow heterogeneity; the data which deviate most from it are the cyan points which correspond to points downstream and in direct contact with the obstacle.
Discussion Figure 3 is compatible with a simple model where an intracellular elastic deformation ε e is in series with the intercellular deformation due to rearrangement rate, resulting in Eε e = η rεr , where E is an intracellular Young modulus and η r an intercellular viscosity due to rearrangements. Such Maxwell viscoelastic liquid picture yields ε e =ε r τ as observed in Fig. 3, with τ = η r /E. There is growing evidence that cell shape and mechanical stress strongly correlate [33,36]. The time τ we find is associated with cell shape relaxation. To verify that the cell shape is a good readout of the mechanical stress in the monolayer, it would be interesting to perform direct stress measurement experiments.
Here τ is much longer than the viscoelastic time associated with internal stress dissipation in the cell due to cytoskeleton viscosity (≈ 15 min [37]) and much shorter than the viscoelastic time associated with cell division (several hours) [38]. As in this experiment it is related with cell rearrangements dynamics (see Fig. 2d), it would be In black,εrτ : deviatoric part of the cell rearrangement rate tensor, multiplied by the value 70 min found in (a); in color, deviatoric part of the cell deformation tensoṙ εe. The color codes for the alignement coefficient c align of these tensors, namely the square cosine of their relative angle, with green corresponding to tensors being aligned (c align = 1) and blue to tensors being orthogonal c align = 0.
interesting to investigate how new cell-cell contact formation dynamics affects the viscoelastic time [7].
We had previously observed velocity waves during the migration of MDCK cells in channels without obstacle; our phenomenological model had indicated that a relaxation time of one to a few hours would be required to explain the onset of velocity and strain waves [15]. The viscoelastic time we measure here, 70±15 min, is consistent with these requirements : it means that cells in the monolayer can sustain strains on the hour timescale.
How stresses generated by molecular motors can generate global tissue flows has been investigated recently [39]. How cellular shapes are affected by these flows depends not only on the viscoelastic time of the tissue but also on the amplitude of the strain rate generated by these flows. The Weissenberg number W e = τε tot is a dimensionless number characterizing the elastic nature of the flow: in tissue regions where W e 1 the flow is quasistatic and cell shapes remain close to their rest shape; conversely, wherever W e is comparable or larger than 1, the material flows but the cells deform due to the tissue strain rate; their shape can strongly vary depending on the tissue flow history. We find a Weissenberg number that is the largest near the obstacle, and reaches W e ≈ 0.5, which is in the transition between these regimes. Hence tuning biologically W e by regulating the cell rearrangements dynamics could be a way to adjust the individual cell geometry memory to the global tissue flow.
A visual example of this memory is the extension of the cell deformation wake downstream of the obstacle (see Fig.  2c); here the source of cell deformation is the strain rate localized just after the obstacle. Cell deformation is advected downstream at a velocity of order of 1 µm/min (see Fig. 2a) and its principal source of decay is cell rearrangements. This competition between deformation transport and shape relaxation yields a typical length scale τ v, which has the same order of magnitude as the correlation length scale of 100 µm usually observed in MDCK cells monolayers [40].
We showed in [15] that in an obstacle free band, the effective one-dimensional migration velocity only depends on local cell density. Here, the two-dimensional velocity field results from the interplay between collective migration in a band and boundary conditions imposed by the obstacle. It has a marked and steady x → −x asymmetry due to density gradient, and strong fluctuations in time resulting in transient y → −y asymmetries. All these properties differ from that of passive materials like foams in similar Stokes flow [26]. However, upstream cells are elongated tangentially to the obstacle and downstream cells are elongated radially, in a way which is visually similar to bubbles in a flowing foam. Our analysis is local in space and integrates the Lagrangian transport term along velocity field lines: thus the results we obtain do not require to have a theoretical prediction of the velocity field.
Conclusion and perspectives The Stokes flow geometry does not only help discriminating between models; it also favors time and ensemble averaging. Measured fields are spatially smooth, their signal to noise ratio sufficient to calculate time and space derivatives. Despite the actual velocity fluctuations, it is possible to compute the transport terms of elastic deformation along a cell trajectory. This enables us to determine all fields involved in a continuum mechanics description.
Since it is based on continuum mechanics descriptions, our analysis facilitates the comparison between in vitro, in vivo and possibly in silico measurements. An experimental technique to introduce an obstacle in vivo in a Drosophila embryo has been recently developed: laser-induced tissue cauterization burns a group of cells, attaching it to the vitelline membrane surrounding the embryo [41]. Such mechanical perturbation can help unveil the underlying cause of the morphogenetic tissue flow. In the same spirit, magnetic fluid drops could be introduced as obstacles in 3D tissues such as in Zebrafish embryo [9,10]. Finally, the present analysis method could be used to analyse natural motions in vivo where a tissue flows around a small organ embedded in it [42].