Active Fingering Instability in Tissue Spreading

TL;DR

This paper demonstrates that epithelial tissue spreading can exhibit an active fingering instability driven by cellular forces and hydrodynamic interactions, leading to fingerlike protrusions during tissue expansion.

Contribution

It introduces an active polar fluid model predicting a generic instability in tissue front dynamics, highlighting the role of hydrodynamics and cellular forces in tissue morphology.

Findings

Active cellular traction induces a velocity gradient causing instability.

Contractile stresses stabilize short wavelengths.

Hydrodynamic interactions determine fingering pattern wavelength.

Abstract

During the spreading of epithelial tissues, the advancing tissue front often develops fingerlike protrusions. Their resemblance to traditional viscous fingering patterns in driven fluids suggests that epithelial fingers could arise from an interfacial instability. However, the existence and physical mechanism of such a putative instability remain unclear. Here, based on an active polar fluid model for epithelial spreading, we analytically predict a generic instability of the tissue front. On the one hand, active cellular traction forces impose a velocity gradient that leads to an accelerated front, which is, thus, unstable to long-wavelength perturbations. On the other hand, contractile intercellular stresses typically dominate over surface tension in stabilizing short-wavelength perturbations. Finally, the finite range of hydrodynamic interactions in the tissue selects a wavelength for the fingering pattern, which is, thus, given by the smallest between the tissue size and the hydrodynamic screening length. Overall, we show that spreading epithelia experience an active fingering instability based on a simple kinematic mechanism. Moreover, our results underscore the crucial role of long-range hydrodynamic interactions in the dynamics of tissue morphology.

Tables

Table 1. Table 1: Estimates of model parameters.

Symbol	Description	Estimate
$L_0$	monolayer half-width	$200$ $\mu$m
$h$	monolayer height	$5$ $\mu$m Trepat et al. (2009); Pérez-González et al. (2019)
$L_c$	nematic length	$25$ $\mu$m Blanch-Mercader et al. (2017); Pérez-González et al. (2019)
$T_0$	maximal traction	$0.5$ kPa Blanch-Mercader et al. (2017); Pérez-González et al. (2019)
$-\zeta$	intercellular contractility	$20$ kPa Pérez-González et al. (2019)
$\xi$	friction coefficient	$100$ Pa$\cdot$s/$\mu$m2 Cochet-Escartin et al. (2014)
$\eta$	monolayer viscosity	$25$ MPa$\cdot$s Blanch-Mercader et al. (2017); Pérez-González et al. (2019)
$\lambda$	hydrodynamic screening length	$0.5$ mm ($\sqrt{\eta /\xi }$)

Equations

$$F=\int \left[\frac{a}{2}{p}_{\alpha }{p}_{\alpha }+\frac{K}{2}\left({\partial }_{\alpha }{p}_{\beta }\right)\left({\partial }_{\alpha }{p}_{\beta }\right)\right]{\text{differential}}^3\vec{r}.$$

$$L_c^2{\nabla }^2{p}_{\alpha }=p_{\alpha },$$

$${\partial }_{\beta }{\sigma }_{\alpha \beta }+f_{\alpha }=0,$$

$${\sigma }_{\alpha \beta }$$

$$f_{\alpha }$$

$$V_0\approx \frac{T_0{L}_c}{2\eta h}{L}_0\equiv \frac{L_0}{\tau },$$

$$\omega \left(q\right)\approx \frac{1}{\tau }+\left[\frac{1}{3\tau }-\frac{\gamma }{2\eta L_0}\right]{\left(q{L}_0\right)}^2+\mathcal{O}\left({\left(q{L}_0\right)}^4\right),$$

$$V_0\approx \frac{T_0{L}_c}{2\eta h}{L}_0+\frac{L_c}{2\eta }\left[\frac{\zeta }{2}-\frac{T_0{L}_c}{h}\right].$$

$$\omega \left(q\right)\approx \frac{1}{\tau }+\left[\frac{1}{3\tau }+\frac{\zeta L_c^2}{8\eta L_0^2}\right]{\left(q{L}_0\right)}^2+\mathcal{O}\left({\left(q{L}_0\right)}^4\right).$$

$$\frac{\partial p_{\alpha }}{\partial t}=-\frac{1}{{\gamma }_1}\frac{\delta F}{\delta p_{\alpha }},$$

$${\partial }_t{p}_{\alpha }=\frac{1}{{\gamma }_1}\left(-a{p}_{\alpha }+K{\nabla }^2{p}_{\alpha }\right).$$

$$L_c^2{\nabla }^2{p}_{\alpha }=p_{\alpha },$$

$$0={\partial }_{\beta }{\sigma }_{\alpha \beta }+f_{\alpha },{\sigma }_{\alpha \beta }={\sigma }_{\alpha \beta }^s+{\sigma }_{\alpha \beta }^a+{\sigma }_{\alpha \beta }^{E,s}$$

$${\sigma }_{\alpha \beta }^E=-P{\delta }_{\alpha \beta }-\frac{\partial f}{\partial \left({\partial }_{\beta }{p}_{\gamma }\right)}{\partial }_{\alpha }{p}_{\gamma },$$

$${\sigma }_{\alpha \beta }={\sigma }_{\alpha \beta }^s+{\sigma }_{\alpha \beta }^a-P{\delta }_{\alpha \beta }.$$

$${\sigma }_{\alpha \beta }={\sigma }_{\alpha \beta }^s+{\sigma }_{\alpha \beta }^a.$$

$${\sigma }_{\alpha \beta }={\sigma }_{\alpha \beta }^s.$$

$$\begin{gathered} {\sigma }_{\alpha \beta }^s=2\eta {\tilde{v}}_{\alpha \beta }+\frac{{\nu }_1}{2}\left(p_{\alpha }{h}_{\beta }+h_{\alpha }{p}_{\beta }-\frac{2}{d}{p}_{\gamma }{h}_{\gamma }{\delta }_{\alpha \beta }\right)-\zeta q_{\alpha \beta } \\ +\left(\bar{\eta }d{v}_{\gamma \gamma }+{\bar{\nu }}_{1}d{p}_{\gamma }{h}_{\gamma }-\bar{\zeta }-{\zeta }^{\prime }{p}_{\gamma }{p}_{\gamma }\right){\delta }_{\alpha \beta }, \end{gathered}$$

$$f_{\alpha }=-\xi v_{\alpha }+{\nu }_i{\dot{p}}_{\alpha }+{\zeta }_i{p}_{\alpha },$$

$${\sigma }_{\alpha \beta }=\eta \left({\partial }_{\alpha }{v}_{\beta }+{\partial }_{\beta }{v}_{\alpha }\right)-\zeta p_{\alpha }{p}_{\beta },$$

$$f_{\alpha }=-\xi v_{\alpha }+{\zeta }_i{p}_{\alpha },$$

$$L_c^2\left({\partial }_x^2+{\partial }_y^2\right)p_x$$

$$L_c^2\left({\partial }_x^2+{\partial }_y^2\right)p_y$$

$${\partial }_x{\sigma }_{xx}+{\partial }_y{\sigma }_{xy}$$

$${\partial }_x{\sigma }_{yx}+{\partial }_y{\sigma }_{yy}$$

$${\sigma }_{xx}$$

$${\sigma }_{xy}$$

$${\sigma }_{yy}$$

$$p_x^0\left(x\right)=\frac{\sinh \left(x/L_c\right)}{\sinh \left(L_0/L_c\right)},$$

$$\begin{gathered} v_x^0\left(x\right)=\frac{\bar{\lambda }}{2\eta }[\zeta +\frac{T_0{L}_c{\bar{\lambda }}^2/h}{{\bar{\lambda }}^2-L_c^2}\coth \left(L_0/L_c\right) \\ -\frac{2\zeta {\bar{\lambda }}^2}{4{\bar{\lambda }}^2-L_c^2}\left[2+{\text{csch}}^2\left(L_0/L_c\right)\right]]\frac{\sinh \left(x/\bar{\lambda }\right)}{\cosh \left(L_0/\bar{\lambda }\right)} \\ +\frac{L_c}{\xi \sinh \left(L_0/L_c\right)} \\ \left[\frac{\zeta }{4{\bar{\lambda }}^2-L_c^2}\frac{\sinh \left(2x/L_c\right)}{\sinh \left(L_0/L_c\right)}-\frac{T_0{L}_c/h}{{\bar{\lambda }}^2-L_c^2}\sinh \left(x/L_c\right)\right], \end{gathered}$$

$$L\left(y\right)=L_0+\delta L\left(y\right).$$

$$p_x\left(x,y\right)=p_x^0\left(x\right)+\delta p_x\left(x,y\right),p_y\left(x,y\right)=\delta p_y\left(x,y\right),$$

$$v_x\left(x,y\right)=v_x^0\left(x\right)+\delta v_x\left(x,y\right),v_y\left(x,y\right)=\delta v_y\left(x,y\right).$$

Full text

The spreading of epithelial monolayers by collective cell migration is crucial for tissue morphogenesis, wound healing, and tumor progression. Both in vivo and in vitro, multicellular protrusions called epithelial fingers often appear at the front of spreading tissues (Fig. 1) Vedula et al. (2013); Saw et al. (2014); Hakim and Silberzan (2017); Ladoux and Mège (2017); Omelchenko et al. (2003); Poujade et al. (2007); Petitjean et al. (2010); Reffay et al. (2011); Klarlund (2012); Reffay et al. (2014); Vishwakarma et al. (2018). This epithelial fingering resembles the viscous fingering that occurs via the Saffman-Taylor instability when a viscous fluid displaces a more viscous one Saffman and Taylor (1958); Casademunt (2004). However, the mechanisms of these two phenomena must be different because epithelial monolayers are more viscous than the fluid that they displace. Hence, several models of epithelial fingering have been proposed Hakim and Silberzan (2017).

To induce finger formation, some models directly implement leader cells with a distinct behavior, either via special particles Sepúlveda et al. (2013) or via a dependence of the magnitude of cell motility forces on the curvature of the tissue front Mark et al. (2010); Tarle et al. (2015). Other models recapitulate epithelial fingering by introducing alignment between cell motility forces and the tissue velocity field Basan et al. (2013). These models predict a moving front to be stable and a non-moving front to exhibit an instability with an unbounded growth rate for a number of finite wavelengths Zimmermann et al. (2014); Nesbitt et al. (2017). Fingers were also observed in the numerical solution of other continuum models of spreading epithelia, either treated as active polar fluids Lee and Wolgemuth (2011) or as active nematics with cell proliferation Doostmohammadi et al. (2015). Recently, fingering was also found in a parameter range of an active vertex model Barton et al. (2017). Finally, interface undulations can emerge from the coupling of chemotactic fields to the mechanics of epithelial spreading Ouaknin and Bar-Yoseph (2009); Salm and Pismen (2012); Köpf and Pismen (2013); Ben Amar and Bianca (2016).

Despite the many efforts, the physical mechanism of the fingering instability in epithelia remains a matter of debate. Here, we address this problem by means of a continuum active polar fluid model for epithelial spreading. The model includes hydrodynamic interactions through the tissue, and it implements neither leader-cell behavior nor alignment between cellular traction forces and the flow field. Yet, we analytically predict a long-wavelength instability of the moving front that explains epithelial fingering. The instability is based on a generic kinematic mechanism, namely the front acceleration associated to a fixed velocity gradient. In spreading epithelia, the velocity gradient is imposed by active traction forces at the edge of the viscous cell monolayer. The fastest-growing mode has a finite wavelength, typically a few hundreds of micrometers, consistent with the measured finger spacing Vishwakarma et al. (2018). This characteristic wavelength is selected by the long-range hydrodynamic interactions in the tissue, which are either limited by the tissue size or screened by cell-substrate friction forces. The model also shows that intercellular contractility stabilizes short-wavelength perturbations of the tissue boundary. The stabilizing effect of contractility is typically stronger than that of tissue surface tension. Globally, our analysis shows how, as a result of the flows induced by the traction force field, a morphological instability may naturally take place in a spreading cell monolayer. Leader cells could then appear upon the onset of the instability, influencing finger development.

Model.—

We base our analysis on a continuum active polar fluid model of epithelial spreading, which is thus described in terms of a polarity field $\vec{p}\left(\vec{r},t\right)$ and a velocity field $\vec{v}\left(\vec{r},t\right)$ Blanch-Mercader et al. (2017); Pérez-González et al. (2019); Alert and Casademunt (2018). We neglect cell proliferation and the bulk elasticity of the monolayer, which eventually limit the spreading process Poujade et al. (2007); Serra-Picamal et al. (2012); Basan et al. (2013); Recho et al. (2016); Yabunaka and Marcq (2017). Tissue spreading is primarily driven by the traction forces exerted by cells close to the monolayer edge, which polarize perpendicularly to the edge by extending lamellipodia towards free space. In contrast, the inner region of the monolayer remains essentially unpolarized, featuring much weaker and transient traction forces Blanch-Mercader et al. (2017); Pérez-González et al. (2019). Hence, we take a free energy for the polarity field that favors the unpolarized state $p=0$ in the bulk, with a restoring coefficient $a>0$, and we impose a normal and maximal polarity as a boundary condition at the tissue edge. In addition, the polar free energy includes a cost for polarity gradients, with $K$ the Frank constant of nematic elasticity in the one-constant approximation de Gennes and Prost (1993). Altogether,

[TABLE]

We assume that the polarity field is set by flow-independent mechanisms, so that it follows a purely relaxational dynamics, and that it equilibrates fast compared to the spreading dynamics SM . Hence, $\delta F/\delta p_{\alpha}=0$, which yields

[TABLE]

where $L_{c}=\sqrt{K/a}$ is the characteristic length with which the polarity modulus decays from $p=1$ at the monolayer edge to $p=0$ at the center.

Then, force balance imposes

[TABLE]

where $\sigma_{\alpha\beta}$ is the stress tensor of the monolayer, and $f_{\alpha}$ is the external force density acting on it. Since tissue spreading occurs over time scales of several hours (Fig. 1), we neglect the elastic response of the tissue SM . Thus, we relate tissue forces to the polarity and velocity fields via the following constitutive equations for an active polar fluid Oriola et al. (2017) (see discussion and justification in Ref. SM ):

[TABLE]

Here, $\eta$ is the effective monolayer viscosity, and $\xi$ is the cell-substrate friction coefficient. Respectively, $\zeta<0$ is the active stress coefficient accounting for the contractility of polarized cells, and $\zeta_{i}>0$ is the contact active force coefficient accounting for the maximal traction stress exerted by polarized cells on the substrate, $T_{0}=\zeta_{i}h$, with $h$ the monolayer height.

Stability of the tissue front.—

To study the stability of the advancing front, we consider a rectangular monolayer typical of in vitro experiments (Fig. 1). Thus, the reference state is the flat front solution with $\vec{p}=p_{x}^{0}(x)\,\hat{x}$ and $\vec{v}=v_{x}^{0}(x)\,\hat{x}$ SM (dashed lines in Fig. 2a). In addition to a maximal normal polarity at the edges, we impose stress-free boundary conditions. For an interface of arbitrary shape, $\vec{p}\left(x=\pm L\right)=\hat{n}_{\pm}$, and $\left.\bm{\sigma}\cdot\hat{n}_{\pm}\right|_{x=\pm L}=\vec{0}$, respectively, where $\hat{n}_{\pm}$ is the normal unit vector of the top and bottom interfaces. The tissue width $L$ changes according to $\differential L/\differential t=\left.\vec{v}\cdot\hat{n}\right|_{x=L}$. Then, motivated by experimental observations (Fig. 1), we introduce peristaltic small-amplitude perturbations of the flat interface, namely those that modify the monolayer width (Fig. 2a): $L(y)=L_{0}+\delta L(y)$. From a linear stability analysis SM , we obtain that the growth rate $\omega(q)$ of such perturbations is always real, so that no oscillatory behavior is expected. However, the growth rate evidences a long-wavelength instability of the monolayer front. Moreover, the fastest-growing perturbation has a finite wavelength (Fig. 2b). In the following, we analyze the contribution of the different forces to the instability, which allows us to single out its physical mechanism.

Traction forces.—

We first consider a limit case with neither intercellular contractility nor cell-substrate friction, $\zeta,\xi\rightarrow 0$ Blanch-Mercader et al. (2017). In addition, we also consider that the width of the polarized boundary layer of cells is much smaller than the total tissue width, $L_{c}\ll L_{0}$, which is generally the case in experiments Blanch-Mercader et al. (2017); Pérez-González et al. (2019). In this limit, since active forces are concentrated at the narrow boundary layer, most of the tissue behaves as a passive viscous fluid, for which $\partial_{x}\sigma_{xx}\approx 0$ and $\sigma_{xx}\approx 2\eta\,\differential v_{x}/\differential x$. Therefore, the stress is uniform throughout most of the tissue, with a value given by the stress accumulated across the boundary layer, namely $\sigma_{xx}\approx T_{0}L_{c}/h$. Consequently, the velocity gradient is also fixed and uniform, and hence the spreading velocity $v_{x}(L_{0})=V_{0}=\differential L_{0}/\differential t$ reads

[TABLE]

where we have used that $v_{x}(0)=0$. This result means that, due to the sole action of a constant traction force, the flat front accelerates, consistent with measurements Rosen and Misfeldt (1980); Poujade et al. (2007) and with the size dependence of tissue wetting Pérez-González et al. (2019). Consequently, the $q=0$ perturbation mode is unstable, $\omega\left(q=0\right)=\tau^{-1}>0$, since any uniform displacement of the advancing front makes it depart from its original velocity. Thus, the instability mechanism is kinematic in nature: The advanced regions of the front move faster than the trailing regions (Fig. 2a). Thereby, traction forces contribute to destabilize perturbations of all wavelengths (Fig. 3a). Therefore, since other forces such as surface tension stabilize short-wavelength perturbations, the tissue front experiences a long-wavelength instability.

Viscous stresses.—

The kinematic mechanism explains why long-wavelength modes are unstable. However, it does not explain why the most unstable mode occurs at a finite wavelength (Fig. 3a). In fact, the existence of a peak in the growth rate is due to the transmission of viscous stress across the monolayer. In the so-called wet limit $\lambda=\sqrt{\eta/\xi}\gg L_{0}$, corresponding to $\xi\rightarrow 0$, viscous stresses transmit through the entire monolayer. Thus, a given perturbation of the front generates a flow perturbation that penetrates a distance of the order of its wavelength, $\pi/q$, into the monolayer. At the monolayer edge, the stress-free boundary condition imposes $\delta\sigma_{xx}(\pm L)=\mp\partial_{x}\sigma_{xx}^{0}(\pm L_{0})\;\delta L$. Hence, since $\delta\sigma_{xx}=2\eta\,\partial_{x}\delta v_{x}$ in the absence of contractility ($\zeta\rightarrow 0$), the gradient of the velocity perturbation profile is fixed at the boundary, being positive (negative) for advanced (trailing) regions of the tissue front. Therefore, flow perturbations further destabilize the flat front in a wavelength-dependent manner.

For short wavelengths, $\pi/q<L_{0}$, the penetration distance of flow perturbations is shorter than the tissue width. Thus, since the slope of velocity perturbations at the interface is fixed, the longer the wavelength, the larger the interfacial velocity perturbation (Fig. 4a). Hence, the growth rate increases with the wavelength (Fig. 3a). In contrast, for long wavelengths, $\pi/q>L_{0}$, the penetration distance of flow perturbations is longer than the tissue width. Thus, in this case, the decay of flow perturbations becomes nearly linear, with a slope that decreases with increasing wavelength. Consequently, perturbations of longer wavelength feature a smaller interfacial velocity perturbation (Fig. 4a), and hence they are less unstable (Fig. 3a). In conclusion, in the absence of cell-substrate friction forces, the finite width of the monolayer limits the range of hydrodynamic interactions in the tissue, thus giving rise to the peak in the growth rate at $q^{*}\sim\pi/L_{0}$ (Fig. 3b).

Cell-substrate friction forces.—

Cell-substrate friction screens the transmission of viscous stresses over distances larger than $\lambda=\sqrt{\eta/\xi}$. Consequently, the peak of the growth rate occurs at $q^{*}\sim\pi/\lambda$ if $\lambda\lesssim L_{0}$ (Fig. 3c). Thus, for sufficiently strong friction, the fingering wavelength is given by the hydrodynamic screening length $\lambda$ instead of the monolayer width $L_{0}$ (Fig. 4b). We estimate $\lambda\sim 0.5$ mm (Table 1). Therefore, the crossover from a viscosity-dominated to a friction-dominated regime of monolayer spreading and fingering should be observable in usual in vitro experiments.

Surface tension.—

The monolayer edge presents a surface tension $\gamma$, namely the work per unit area required to expand it. For a curved interface, surface tension gives rise to a normal force: $\left.\hat{n}_{\pm}\cdot\bm{\sigma}\cdot\hat{n}_{\pm}\right|_{x=\pm L}=-\gamma\vec{\nabla}\cdot\hat{n}_{\pm}|_{x=\pm L}\approx\pm\gamma\,\differential^{2}\delta L/\differential y^{2}$. Then, the expansion of the growth rate at long wavelengths reads

[TABLE]

with $\tau=2\eta h/(T_{0}L_{c})$. This expression reveals the existence of a critical size $L_{\gamma}\approx 3\gamma h/(T_{0}L_{c})$ above which the growth rate curves upwards at $q\rightarrow 0$ ($\omega^{\prime\prime}(0)>0$), hence exhibiting the aforementioned peak at a finite wavelength. Alternatively, if $\gamma>\gamma^{*}\approx T_{0}L_{c}L_{0}/(3h)$, surface tension prevents the hydrodynamic selection of a finite fingering wavelength, which is then only limited by the length of the tissue front.

The surface tension of the monolayer could be due to actin cables found along its edge, specially along the sides of epithelial fingers Poujade et al. (2007); Klarlund (2012); Reffay et al. (2014). Traction force measurements suggest that the tension of such cables is $\gamma\sim 0.2$ mN/m Reffay et al. (2014), lower than typical surface tensions of cell aggregates, $\gamma\sim 1-10$ mN/m Foty et al. (1994); Forgacs et al. (1998); Guevorkian et al. (2010); Stirbat et al. (2013). Combining these value ranges with $T_{0}=0.2-0.8$ kPa and typical values of $h$ and $L_{c}$ (Table 1), the critical monolayer width for fingering is $L_{\gamma}\sim 0.3-10$ $\mu$m. Therefore, we expect surface tension not to play a major role in the fingering instability in monolayers of typical widths $L_{0}\sim 0.1-1$ mm.

Intercellular contractility.—

Because it decreases the spreading velocity, intercellular contractility has an additional stabilizing effect on the monolayer front (Fig. 3d). We discuss the effects of a uniform contractility in SM . Here, we consider an intercellular contractility $-\zeta$ concentrated at the polarized boundary layer, which has a size-independent contribution to the spreading velocity:

[TABLE]

Consequently, in the limit $L_{c}\ll L_{0}\ll\lambda$, this contractility has no impact on the growth rate of the uniform mode, $\omega\left(0\right)=\tau^{-1}$, but it contributes a stabilizing quadratic term to the long-wavelength expansion of the growth rate:

[TABLE]

Thus, as surface tension, contractility defines a critical size $L_{\zeta}\approx\sqrt{-3\zeta L_{c}h/(4T_{0})}$ above which the growth rate features a finite-wavelength peak. Alternatively, if $-\zeta>-\zeta^{*}\approx 4T_{0}L_{0}^{2}/(3hL_{c})$, contractility supresses the hydrodynamic selection of a finite fingering wavelength.

For typical contractilities, $-\zeta\sim 1-50$ kPa Pérez-González et al. (2019), the critical monolayer width for fingering is $L_{\zeta}\sim 10-100$ $\mu$m, where we have used $T_{0}=0.2-0.8$ kPa and estimates for $h$ and $L_{c}$ (Table 1). Therefore, we do not expect contractility to prevent fingering wavelength selection. However, our estimates give $L_{\gamma}<L_{\zeta}$, indicating that contractility typically dominates over surface tension in stabilizing short-wavelength shape perturbations. Thus, the competition between the destabilizing effect of traction forces and the stabilizing effect of intercellular contractility defines the band of unstable modes.

Uniformly polarized tissues.—

To consider monolayers with bulk polarity Trepat et al. (2009), we analyze the morphological stability of a uniformly polarized monolayer SM . In this case, in the absence of contractility, the spreading velocity is constant: $V_{0}=T_{0}/(\xi h)$. Hence, the $q=0$ mode is marginally stable, $\omega(q=0)=0$. Moreover, contractility stabilizes it, $\omega(0)<0$. Nevertheless, the viscous effects discussed above still give rise to a peak of the growth rate at a finite wavelength. Therefore, even though uniformly polarized tissues do not feature an accelerating front, they still exhibit a fingering instability for sufficiently small contractility SM .

Conclusions.—

Motivated by the observation of finger-like protrusions during the spreading of epithelial monolayers, we studied the stability of the advancing front. Modeling the cell monolayer as an active polar fluid, we showed that active traction forces are responsible for a long-wavelength instability of the monolayer front. Several features distinguish this instability from previous proposals. First, it is generic; it takes place for any value of the active traction force. Second, the wavelength of the fingering pattern is selected by the range of hydrodynamic interactions in the tissue. And third, active intercellular forces stabilize short-wavelength perturbations, typically dominating over surface tension effects.

Our analysis identifies the physical mechanism of the instability. Cellular traction forces at the monolayer edge set the velocity gradient in the spreading monolayer. Hence, under the same traction force, a larger monolayer spreads faster Pérez-González et al. (2019). Consequently, when the monolayer front is perturbed, the protruding regions of the interface advance faster than the trailing regions, thus making the perturbation grow. Therefore, the instability is based on a simple kinematic mechanism, which takes place generically in viscous fluids that sustain a fixed velocity gradient in the direction of spreading. In particular, the same morphological instability should occur in the so-called squeeze flow Engmann et al. (2005), in which an incompressible fluid is forced to spread by decreasing the gap between two plates. In this case, under perfect slip conditions at the plates, the rate of gap reduction sets the fixed velocity gradient.

Regarding spreading epithelia, we conclude that neither leader-cell behavior nor regulation of cell motility by curvature or by chemotactic fields are necessary for the fingering instability. Therefore, our results are consistent with the emergence of leader cells concomitantly with finger growth Poujade et al. (2007); Petitjean et al. (2010); Reffay et al. (2011, 2014); Vishwakarma et al. (2018). However, the viscous rheology of the monolayer is essential for the instability. On the one hand, it underpins the velocity gradient that renders the interface unstable and, on the other, it enables wavelength selection for the fingering pattern. Concomitant with the fingering instability, shear stresses give rise to flows transversal to the spreading direction, which might lead to the swirls observed in experiments Poujade et al. (2007); Petitjean et al. (2010). Finally, in addition to explaining fingering in tissue spreading, our results also account for the morphological instability recently observed during tissue dewetting Pérez-González et al. (2019).

Our predictions, such as the absence of a traction force threshold for the instability, and whether the fingering wavelength is given by either the monolayer width or the screening length $\lambda=\sqrt{\eta/\xi}$, are experimentally testable. Indeed, consistent with our result, recent work has shown that the finger spacing is an intrinsic quantity that coincides with the stress correlation length Vishwakarma et al. (2018). To further test our predictions, future experiments could perturb active cellular forces, cell-cell and cell-substrate adhesion, and vary monolayer width.

Our findings illustrate how hydrodynamic interactions impact tissue morphodynamics. In particular, we propose that epithelial fingering can naturally arise from a generic morphological instability in a fluid film driven by interfacial active forces. Thus, our results showcase the relevance of interfacial instabilities in driven Troian et al. (1989); Melo et al. (1989); Ben Amar and Cummings (2001) and active Callan-Jones et al. (2008); Sankararaman and Ramaswamy (2009); Sarkar and Basu (2012, 2013); Nagilla et al. (2018); Williamson and Salbreux (2018); Bogdan and Savin (2018) fluids for tissue spreading.

Appendix A ACTIVE POLAR FLUID MODEL OF EPITHELIAL SPREADING

Here, we briefly justify the description of epithelial spreading in terms of our continuum active polar fluid model.

A.1 Polarity dynamics

The outwards polarization of cells at the monolayer edge is likely due to contact inhibition of locomotion, a cell-cell interaction whereby cells repolarize in opposite directions upon contact Mayor and Carmona-Fontaine (2010); Stramer and Mayor (2017). In fact, this interaction is mediated by cell-cell adhesion, with front-rear differences in cadherin-based junctions acting as a cue for the repolarization Desai et al. (2009); Khalil and Friedl (2010); Weber et al. (2012); Theveneau and Mayor (2013); Vedula et al. (2013); Ladoux et al. (2016). Although originally proposed for mesenchymal cells, contact inhibition of locomotion is being increasingly recognized to play a key role in orchestrating the collective migration of epithelial monolayers Mayor and Carmona-Fontaine (2010); Theveneau and Mayor (2013); Vedula et al. (2013); Ladoux et al. (2016); Mayor and Etienne-Manneville (2016); Hakim and Silberzan (2017); Zimmermann et al. (2016); Coburn et al. (2016); Smeets et al. (2016). In a cohesive monolayer, this interaction naturally leads to polarization of cells at the edge towards free space, leaving the inner region of the monolayer unpolarized. Such a polarity profile, in turn, explains the localization of traction forces at the edge and the build-up of tension at the center of epithelial monolayers Zimmermann et al. (2016); Coburn et al. (2016). Therefore, we assume the polarity field $\vec{p}\left(\vec{r},t\right)$ to be set by an autonomous cellular mechanism such as contact inhibition of locomotion, which polarizes cells within a time scale $\tau_{\text{CIL}}\sim 10$ min Smeets et al. (2016); Weber et al. (2012). Hence, $\vec{p}\left(\vec{r},t\right)$ should remain essentially independent of flows in the monolayer, which occur over a longer time scale given by the strain rate, at least of order $\tau_{s}\sim 100$ min Blanch-Mercader et al. (2017); Vincent et al. (2015). Consequently, taking a phenomenological approach, we propose the polarity field to follow a purely relaxational dynamics given by

[TABLE]

where $F\left[\vec{p}\,\right]$ is the coarse-grained free energy functional for the polarity field (Eq. 1), and $\gamma_{1}$ is a kinetic coefficient (the rotational viscosity for the angular degree of freedom). With respect to the most general dynamics of the polarity field in an active polar fluid, Eq. S1 neglects polarity advection and corotation, as well as flow alignment and active spontaneous polarization effects. Thus, using Eq. 1, the dynamics of the polarity is given by

[TABLE]

In the limit of fast polarity dynamics compared to the spreading dynamics, the polarity field is always at equilibrium, $\partial_{t}p_{\alpha}=0$, adiabatically adapting to the shape of the monolayer. Under this approximation, the polarity field is given by

[TABLE]

where we have defined the characteristic length $L_{c}\equiv\sqrt{K/a}$ of the polar order in the monolayer.

A.2 Force balance

Flows in cell monolayers occur at very low Reynolds numbers. Therefore, inertial forces are negligible, and hence momentum conservation reduces to the force balance condition

[TABLE]

where $\sigma_{\alpha\beta}^{s}$ and $\sigma_{\alpha\beta}^{a}$ are the symmetric and antisymmetric parts of the deviatoric stress tensor, and $f_{\alpha}$ is the external force density. Respectively, $\sigma_{\alpha\beta}^{E,s}$ is the symmetric part of the Ericksen tensor. This tensor generalizes the pressure $P$ to include anisotropic elastic stresses associated to the orientational degrees of freedom in liquid crystals de Gennes and Prost (1993):

[TABLE]

where $f$ is the Frank free energy density, namely the integrand of Eq. 1. Thus, the orientational contribution to the Ericksen tensor is of second order in gradients of the polarity field, and hence we neglect it, so that the stress tensor reads

[TABLE]

Then, the pressure is related to the cell number surface density $\rho$ by the equation of state of the monolayer. For the sake of an estimate, we assume the simplest form for an equation of state, $P\left(\rho\right)=B\left(\rho-\rho_{0}\right)/\rho_{0}$, where $B$ is the bulk modulus of the monolayer, and $\rho_{0}$ is a reference density defined by $P\left(\rho_{0}\right)=0$. Taking the pressure origin at the monolayer edge, $\rho_{0}\sim 2\cdot 10^{3}$ cells/mm2 Trepat et al. (2009); Pérez-González et al. (2019). In turn, density differences in the monolayer are, at most, $\rho-\rho_{0}\sim 6\cdot 10^{3}$ cells/mm2 Trepat et al. (2009); Pérez-González et al. (2019). Then, the monolayer is expected to be highly compressible because area changes can in principle be accommodated by changes in height, resisted only by the shear modulus of the tissue. Hence, we estimate the bulk modulus of the monolayer by typical shear moduli of cell aggregates, which are in the range $G\sim 10^{2}-10^{3}$ Pa Forgacs et al. (1998); Marmottant et al. (2009); Guevorkian et al. (2010). Thus, the pressure in the monolayer should be $P\lesssim 30-300$ Pa. In fact, isotropic compressive stresses (pressures) of $\sim 50$ Pa were shown to induce cell extrusion Saw et al. (2017). In conclusion, if tissue spreading is not dominated by cell proliferation Basan et al. (2013); Recho et al. (2016); Yabunaka and Marcq (2017), the magnitude of the pressure in the monolayer is expected to be much smaller than the tensile stress (tension) induced by traction forces, as measured by monolayer stress microscopy, which is of the order of kPa Trepat et al. (2009); Pérez-González et al. (2019). Hence, we neglect pressure:

[TABLE]

Now, for a nematic medium, the antisymmetric part of the stress tensor is given by $\sigma_{\alpha\beta}^{a}=1/2\left(p_{\alpha}h_{\beta}-h_{\alpha}p_{\beta}\right)$, where $h_{\alpha}=-\delta F/\delta p_{\alpha}$ is the molecular field. From Eq. S1, the adiabatic approximation for the polarity dynamics, $\partial_{t}p_{\alpha}=0$, implies $h_{\alpha}=0$. Therefore, the antisymmetric part of the stress tensor vanishes under this approximation, $\sigma_{\alpha\beta}^{a}=0$. Thus, the stress tensor reduces to

[TABLE]

A.3 Constitutive equations

Next, constitutive equations must be given to specify the deviatoric stress tensor $\sigma_{\alpha\beta}^{s}$ and the external force $f_{\alpha}$ in terms of the polarity and velocity fields. The generic constitutive equations of an active liquid crystal are provided by active gel theory Kruse et al. (2005); Jülicher (2011); Marchetti et al. (2013); Prost et al. (2015). Here, based on the previous assumptions for the dynamics of the polarity field, we propose a simplified version of the generic constitutive equations of an active polar gel to describe epithelial spreading.

First, the spreading occurs on timescales of the order of $\tau_{s}\sim 100$ min Blanch-Mercader et al. (2017); Vincent et al. (2015), at which the tissue should have a fluid rheology. This time scale is much slower than the turnover time scales of proteins in the cytoskeleton or in cell-cell junctions, which are of the order of tens of minutes at most Wyatt et al. (2016); Khalilgharibi et al. (2016). Intra- or intercellular processes such as cytoskeletal reorganizations or cell-cell slidings dissipate energy over these time scales, so that elastic energy may only be stored in the tissue at shorter times. In addition, other processes such as cell division, death, and extrusion Ranft et al. (2010); Matoz-Fernandez et al. (2017), as well as cell shape fluctuations Marmottant et al. (2009); Etournay et al. (2015) and topological rearrengements Etournay et al. (2015); Krajnc et al. (2018) also fluidize the tissue. Therefore, to describe the slow spreading dynamics, we will not consider the elastic response of the tissue at short time scales.

Then, in the viscous limit, the constitutive equations for the internal stress and the interfacial force of an active polar medium are:

[TABLE]

[TABLE]

where, $q_{\alpha\beta}=p_{\alpha}p_{\beta}-p_{\gamma}p_{\gamma}/d\;\delta_{\alpha\beta}$ is the traceless symmetric nematic order parameter tensor, with $d$ the system dimensionality, and $v_{\alpha}$ is the velocity of the fluid with respect to the substrate. The coefficients $\eta$ and $\bar{\eta}$ are the shear and bulk viscosities of the medium, $\zeta$ is the anisotropic active stress coefficient, and $\bar{\zeta}$ and $\zeta^{\prime}$ are two isotropic active stress coefficients. Finally, $\xi$, $\nu_{i}$, and $\zeta_{i}$ are the corresponding interfacial versions of the viscosity (viscous friction), flow alignment (polar friction), and active stress (active force) coefficients. The constitutive equation for the internal stress, Eq. S9, is that of an active polar gel with a variable modulus of the polarity Jülicher (2011). In turn, the constitutive equation for the interfacial force, Eq. S10, is less conventional Jülicher and Prost (2009), but it was derived from a mesoscopic model of an active polar gel Oriola et al. (2017).

Now, the adiabatic approximation for the polarity dynamics implies $\dot{p}_{\alpha}=h_{\alpha}=0$, so that flow alignment terms contribute neither to the stress tensor nor to the interfacial force. Next, we assume that polarized cells generate much larger active stresses than unpolarized cells. Hence, we neglect the active stress coefficient $\bar{\zeta}$ in front of $\zeta$ and $\zeta^{\prime}$. Then, assuming that the two-dimensional fluid layer is compressible, we take $\zeta=\zeta^{\prime}\,d=2\zeta^{\prime}$ and $2\eta=\bar{\eta}\,d=2\bar{\eta}$ for simplicity. Under these simplifications, and using Eq. S8, the constitutive equations reduce to

[TABLE]

[TABLE]

which close the set of equations defining the active polar fluid model of the spreading of an epithelial monolayer.

Appendix B LINEAR STABILITY ANALYSIS

Here, we give the details of the linear stability analysis of the tissue front. First, we explicitly write down the equations of the model in Cartesian coordiates, which are most convenient for the rectangular geometry of the monolayer (Fig. 2a). Thus, the equation for the polarity field, Eq. 2, reads

[TABLE]

Respectively, the force balance equation Eq. 3 reads

[TABLE]

where the components of the stress tensor are given by

[TABLE]

Next, we obtain the flat front solution in rectangular geometry, which is the reference state of the linear stability analysis. The long ($\hat{y}$) axis of the rectangle is much longer than the short ($\hat{x}$) axis. Hence, we assume translational invariance along the long axis of the monolayer Blanch-Mercader et al. (2017). Moreover, traction forces are mainly perpendicular to the monolayer boundary Trepat et al. (2009), so that we take the polarity field along the $\hat{x}$ direction: $\vec{p}=p_{x}^{0}\left(x\right)\hat{x}$, where the superindex indicates the zeroth order in the perturbations of the front. Now, imposing symmetry as well as maximal polarity and stress-free boundary conditions, $p_{x}^{0}\left(L_{0}\right)=1$ and $\sigma_{xx}^{0}\left(L_{0}\right)=0$, one obtains the polarity and velocity profiles:

[TABLE]

[TABLE]

where $\bar{\lambda}=\sqrt{2\eta/\xi}=\sqrt{2}\,\lambda$ is a redefined hydrodynamic screening length, and $L_{0}$ is the semi-width of the monolayer, which changes according to $\differential L_{0}/\differential t=v_{x}^{0}\left(L_{0}\right)$.

Next, we introduce peristaltic small-amplitude perturbations of the flat interface of the monolayer (Fig. 2a):

[TABLE]

Under these perturbations, the polarity and velocity fields take the form

[TABLE]

In turn, boundary conditions must keep imposing a normal and maximal polarity, as well as vanishing normal and shear stresses at the interface, which is now curved. To this end, we define the normal and tangential vectors of each interface,

[TABLE]

where $\theta$ is the angle between the normal directions of the flat and perturbed interfaces, and the $\pm$ index stands for the top and bottom interfaces, respectively (Fig. 2a). Thus, the boundary conditions for the polarity read

[TABLE]

For the $x$-component, the conditions imply $p_{x}\left(\pm L\right)\approx\pm 1$. This expands into

[TABLE]

which yields

[TABLE]

as a boundary condition on the polarity perturbation. For the $y$-component of the polarity perturbation, the boundary condition imposes

[TABLE]

Then, the boundary conditions on the stress read

[TABLE]

which respectively ensure vanishing normal and shear stress at the interfaces. Here, for simplicity, we neglect interfacial tension and bending rigidity, which would contribute stabilizing terms to the growth rate. The condition on the normal stress gives $\sigma_{xx}\left(\pm L\right)=0$ which, after expanding as previously, leads to

[TABLE]

for the stress perturbation. In turn, the condition on the shear stress directly gives

[TABLE]

Next, we decompose all perturbations in their Fourier modes, identified by the wave number $q$:

[TABLE]

In terms of the Fourier modes, the equations for the polarity components read

[TABLE]

In turn, the components of the force balance equation, once the constitutive relation is introduced, read

[TABLE]

The boundary conditions must also be translated into the Fourier domain, reading

[TABLE]

[TABLE]

Then, the four coupled differential equations Eqs. S30 and S31 are analytically solved for $\delta\tilde{p}_{\alpha}\left(x,q\right)$ and $\delta\tilde{v}_{\alpha}\left(x,q\right)$. From the Fourier modes of the velocity field, the perturbed spreading velocity $V$ can be computed as

[TABLE]

so that

[TABLE]

Thus, the growth rate $\omega\left(q\right)$ of the tissue shape perturbations follows from

[TABLE]

Hence,

[TABLE]

The expression of the resulting growth rate is omitted here due to its length. Finally, note that, in our free-boundary problem, the amplitude of the front perturbations does not grow exponentially in time. This is because the growth rate depends on time through the monolayer width $L_{0}\left(t\right)$. Consequently, Eq. S36 yields

[TABLE]

Appendix C NEMATIC ELASTICITY

In this section, we discuss the effects of the nematic elasticity of the polarity field on the growth rate of front shape perturbations. Front perturbations distort the polarity field, generating polarity gradients along the tissue front. The elastic energy cost of these transversal polarity gradients is larger for shorter-wavelength perturbations. Thus, to minimize the polar free energy Eq. 1, polarity perturbations decay more steeply for shorter wavelengths. From Eq. 2, the Fourier components of polarity perturbations obey $L_{c}^{2}\left(\partial_{x}^{2}-q^{2}\right)\delta\tilde{p}_{\alpha}=\delta\tilde{p}_{\alpha}$, so that their decay length is $\ell_{c}(q)=L_{c}[1+(qL_{c})^{2}]^{-1/2}$, which decreases with $q$. Similarly to Eq. 5, the corresponding velocity gradient perturbation is then proportional to $T_{0}\ell_{c}(q)$. Therefore, by reducing the size of the polarized boundary layer, $\ell_{c}(q)\leq L_{c}$, nematic elasticity causes a decrease of the growth rate with decreasing wavelength. However, this effect is only significant at wavelengths shorter than the nematic length $L_{c}$ ($qL_{c}>1$), which is typically smaller than the size of the fingers.

Appendix D UNIFORM INTERCELLULAR CONTRACTILITY

In this section, we discuss the effects of a uniform intercellular contractility term, characterized by the coefficient $-\bar{\zeta}$ in Eq. S9, which is neglected in the Main Text. Like the intercellular contractility $-\zeta$, which is localized at the polarized boundary layer of the tissue, a uniform contractility has a stabilizing effect on the tissue front (Fig. S1). However, unlike its polarity-related counterpart $-\zeta$, the uniform contractility $-\bar{\zeta}$ has a size-dependent contribution to the spreading velocity:

[TABLE]

in the limit $L_{c}\ll L_{0}\ll\lambda$. Then, the expansion of the growth rate at long wavelengths reads

[TABLE]

where $\tau=2\eta h/(T_{0}L_{c})$. This expression shows that a uniform contractility does not only affect the growth rate of finite-wavelength perturbations but that it also decreases the growth rate of the uniform mode, $\omega(0)$. If the contractility is sufficiently small to allow the tissue to spread ($V_{0}>0$), $-\bar{\zeta}<T_{0}L_{c}/h$, the tissue front remains unstable to long-wavelength perturbations. However, unlike the polarity-related contractility $-\zeta$, if the uniform contractility $-\bar{\zeta}$ induces the retraction of the monolayer front, it also prevents its fingering instability (Fig. S1).

Appendix E LINEAR STABILITY ANALYSIS FOR UNIFORMLY POLARIZED TISSUES

In this section, we address situations in which the epithelial monolayer is polarized not only in a boundary layer but also in its bulk Trepat et al. (2009). To analyze the effects of bulk polarity, we consider the simplest possible situation, namely a uniformly polarized epithelium. The state of uniform polarization, however, presents two issues. First, it may be unstable in the bulk Blanch-Mercader and Casademunt (2017). Here, we assume that the active contractility and traction force coefficients are small enough not to trigger this bulk instability. Second, because the polarity points outward at both fronts of the monolayer (Fig. 2), a spreading monolayer cannot have a uniform polarity. Following previous works Zimmermann et al. (2014); Nesbitt et al. (2017), we avoid this issue by considering the tissue to be bounded by a comoving wall at $x=0$. In this case, we can assume a fixed modulus of the polarity field, $|\vec{p}\,|=1$. Hence, the polar free energy is left only with Frank elasticity:

[TABLE]

Consequently, the equilibrated polarity field, which fulfills $\delta F/\delta p_{\alpha}=0$, is a solution of

[TABLE]

The rest of the model is unchanged with respect to the case with only a polarized boundary layer analyzed in the Main Text. Thus, the force balance is still given by Eqs. 3-4, and the boundary conditions still impose a normal polarity and a vanishing stress at the monolayer front: $\vec{p}(x=L)=\hat{n}$ and $\left.\bm{\sigma}\cdot\hat{n}\right|_{x=L}=0$, respectively, with $\hat{n}$ being the normal unit vector of the tissue front. Moreover, we impose the following boundary conditions at the comoving wall: $\vec{p}(x=0)=\hat{x}$, and $\vec{v}(x=0)=T_{0}/(\xi h)\,\hat{x}$. Thus, the reference state of the stability analysis is the flat front solution with

[TABLE]

where, as previously defined, $\bar{\lambda}=\sqrt{2\eta/\xi}$.

The linear stability analysis can thus be performed as before. In the present case, instead of Eq. S30, the equations for the Fourier modes of the polarity perturbations read

[TABLE]

Respectively, the components of the force balance equation, Eq. S31, reduce to

[TABLE]

Finally, the boundary conditions, previously given by Eqs. S32 and S33, now become

[TABLE]

[TABLE]

Then, Eqs. S44 and S45 are analytically solved for $\delta\tilde{p}_{\alpha}(x,q)$ and $\delta\tilde{v}_{\alpha}(x,q)$. Hence, the growth rate of front perturbations is computed using Eq. S37. The resulting expression is omitted here due to its length. However, the growth rate is plotted in Fig. S2, which shows how it changes under variation of different parameters. From these results, we conclude that, as for tissues with only a polarized boundary layer, uniformly polarized tissues also display an active fingering instability. However, in contrast to the boundary layer case, in the absence of contractility, uniformly polarized tissues spread at a constant velocity $V_{0}=T_{0}/(\xi h)$. As a consequence, the $q=0$ mode, is now marginally stable. Nonetheless, all other modes are destabilized by the interplay between active traction forces and viscous stresses. As for the case with an active boundary layer, front perturbations give rise to flow perturbations that destabilize the flat front (Fig. S2a). Hence, the active fingering instability is robust to the presence of bulk polarity in the tissue.

Moreover, the role of hydrodynamic interactions is unchanged with respect to the active boundary layer case. As in that case, the finite range of hydrodynamic interactions, given by the smallest between the monolayer width $L_{0}$ and the hydrodynamic screening length $\lambda=\sqrt{\eta/\xi}$ determines the wavelength of the most unstable mode. Thus, for polarized tissues, varying the monolayer width $L_{0}$ and the cell-substrate friction coefficient $\xi$ modifies the growth rate in a similar way as for tissues with only a polarized boundary layer (compare Fig. S2b-c to Fig. 3b-c).

Finally, the effect of the contractility in uniformly polarized tissues is a bit different than in unpolarized tissues. As for tissues with only boundary polarity, contractility has a stabilizing effect on the monolayer front. However, for uniformly polarized tissues, the active contractile stress spans throughout the monolayer. As a consequence, instead of just affecting short-wavelength perturbations, contractility now decreases the growth rate of all perturbation modes. Hence, for small contractility, the longest and shortest-wavelength modes become stable, leaving a band of unstable modes at intermediate wavelengths. Therefore, small contractilities do not abrogate the fingering instability. A sufficiently large contractility, however, is able to stabilize all modes, thus suppressing the fingering instability (Fig. S2d).