Morphological instabilities of stratified epithelia: a mechanical instability in tumour formation

Thomas Risler; Markus Basan

doi:10.1088/1367-2630/15/6/065011

1. Introduction

Most human cancers are carcinomas, tumours that originate in epithelial tissues [1]. From there, they may invade through the basement membrane into the supporting connective tissue—the stroma—and eventually lead to the formation of metastases [2, 3]. The process of invasion and the modes of motility associated with it are a subject of intense study. Depending on cancer type, cells may break away from the primary tumour and migrate to the blood vessels as single cells, collectively as detached clusters or as multicellular, three-dimensional invasive strands [4, 5]. Undoubtedly, the acquisition of mesenchymal traits by the invasive cancer cells, together with cell motility and the expression of proteases play central roles in many of these cases and particularly in malignant tumours [6–9]. However, undulations and protrusions of the epithelium into the stroma are commonly found in benign tumours and even in healthy epithelia, where the basement membrane separating the two tissues remains intact. As the tissue progresses to higher grades of malignancy, the number of dividing layers within the epithelium increases and such protrusions typically grow in size [1, 10]. This is the case for example in cervical intraepithelial neoplasia [10, 11] and in the epithelial dysplasia of the oral mucosa [12, 13]. In this work, we investigate a potential mechanism underlying these commonly observed, yet striking morphological features of stratified epithelia and carcinomas, which relies on proliferation-induced mechanical stresses without invoking proper cell motility. The mechanics of this process is of particular interest, considering that these protrusions arise in epithelia whose apical surfaces are free of stress, such that it is non-trivial how an increased proliferation rate alone like that observed in neoplastic tissues can lead to invasive protrusions into the stroma.

In previous work [14], we have demonstrated the existence of a mechanical instability based entirely on the flow caused by cell renewal in the epithelium in combination with its viscosity due to cell–cell adhesion. We showed that when the epithelium–stroma interface is displaced from an originally flat configuration, the excess of cell division above a nascent protrusion creates a differential flow of cells, and the resulting viscous shear stress drives the protrusion further. A finite-wavelength instability develops from the combination of surface-stabilizing factors such as interfacial tension and mechanical resistance of the stroma on the one hand, and destabilizing factors enhancing cell-division-driven flows and viscous shear stresses on the other hand. In our previous work, we studied the two limits of a purely elastic or purely viscous rheology of the stroma. This paper serves as a more detailed presentation of this instability as well as proposes two important extensions of our earlier work. Firstly, we study the instability in the general case of a viscoelastic stroma and investigate the transition between the two regimes presented previously, which we recover as limits. Secondly, and most importantly, we include a dependence of the cell-division rate in the epithelium on the local concentration of a substance necessary for cell division—representing nutrients, growth factors or oxygen—which diffuses through the epithelium from the stroma where it is delivered by blood vessels [15]. Compared with our previous study where the rate of cell renewal in the epithelium was pre-defined as a function of distance from the stroma, this coupling leads to a significant enhancement of the instability by a diffusion-limitation mechanism, which is similar to mechanisms leading to instabilities in other contexts [16–19]. We comment further in the discussion on the comparison with classical instabilities known from non-equilibrium physics as well as on other dynamical instabilities observed in living systems such as bacterial colonies.

In tissues, mechanical instabilities may play a role in the morphogenesis of certain shapes and patterns exhibited by growing cell populations. For example, differential growth has been proposed as a mechanism underlying the large-scale looping morphology of the gut [20], and a buckling instability of a monolayered epithelium has been suggested for the formation of villi and crypts in the colon [21, 22]. In the case of multilayered epithelia, a buckling instability of the basal layer of the fetal epidermis has been proposed to be at the origin of the formation of epidermal ridges [23]. The complex network of finger-like protrusions at the dermal–epidermal junction of human skin has been proposed to result from incompatible growth of elastic tissue layers [24]. Similar instabilities have been investigated in different geometries, such as that of a growing elastic shell with differential growth [25–29] or that of a tube [30]. These studies are part of the broader field related to the role of mechanics for growing tissues and of their microenvironment [31–34].

In the specific case of tumour growth, pattern formation has been reported and studied theoretically using reaction–diffusion descriptions of the supply of nutrients, oxygen or growth factors [35]. Such a coupling affects tumour-growth dynamics [36], as well as tumour shape and interfacial structure [37]. Fingering instabilities may develop depending on the degree and spatial structure of vascularization [38]. Taking into account the mechanical properties of the growing tissue and that of its microenvironment, residual stresses [39] and fingering instabilities may or may not develop [40]. In particular, weakened cell–cell adhesions and cell–matrix adhesion have been proposed to favour such instabilities [41, 42]. Instabilities may also result from a combination of pushing forces due to cell proliferation and pulling forces of the tissue on either a substrate or the extra-cellular matrix in the case of a three-dimensional, multi-cellular spheroid [43]. However, the coupling of growth dynamics to the diffusion field of nutrients has never been investigated within the framework of the instability proposed here, which differs from the aforementioned studies because it is based on a purely hydrodynamic description of the tissue and because it arises in the specific geometry of an epithelium with a free apical surface rather than that of a tumour spheroid or a two-dimensional tissue spreading on a substrate. Because this additional diffusion-dependent mechanism can lead to an unstable epithelium–stroma interface over a broader range of parameters than the instability reported previously [14], we expect this mechanism to potentially play an important role in the generation of undulated patterns as they occur in real epithelial tissues.

In this paper, we present six different versions of a linear stability analysis of the epithelium–stroma interface for different rheologies of the stroma and cell-renewal functions in the epithelium. In section 2, we start by presenting the common properties of the six models that will be studied in the paper. The associated detailed equations can be found in the supplementary information (available from stacks.iop.org/NJP/15/065011/mmedia). Models 1 and 2, where the cell-production function in the epithelium is pre-defined and the stroma is either purely elastic (model 1) or purely viscous (model 2), have already been presented in [14]. The results associated with these models are therefore presented in the supplementary information including some novel aspects as compared with our previous work [14]. In section 3, we generalize the two previous cases to that of a viscoelastic rheology of the stroma (model 3), and compare the mode structures associated with these three different models. Section 4 is devoted to the description of the coupling of cell division in the epithelium to the local concentration of available nutrients, which diffuse from the stroma and are consumed by the epithelial cells. In sections 5–7, we present the resulting stability analysis successively in the case of an elastic (model 4), a viscous (model 5) and a viscoelastic (model 6) stroma. Domains of validity, similarities and differences between the various models are discussed throughout the presentation. We finally discuss the similarities and differences between the results presented in this study and other instabilities known from the literature.

2. Description of models 1–3

We start this section by presenting the geometry of the system, which underlies the six different models presented in this paper. We then present the structure of models 1–3, where the cell-production function in the epithelium is pre-imposed as a function of distance from the epithelium–stroma interface. The structure of the models is just described here, and the full set of equations can be found in the supplementary information. The results associated with models 1 and 2, already present for the main part in our previous report [14], can be found in the supplementary information. We present here in more detail the results obtained with model 3, which constitutes a first extension of our previous work.

2.1. Description of the geometry

We consider an epithelium of thickness H adjacent to a stroma of thickness L, infinitely extended in directions x and y. In the z-axis, we assume a rigid wall at z = 0 to which the stroma is anchored, followed by the stroma, then the epithelium starting at z = L, and finally the apical surface of the epithelium directly adjacent to a lumen at z = L + H (figure 1). These coordinates correspond to the flat configuration and are strictly valid only in the unperturbed, steady state. We study the stability of these interfaces under a given arbitrary but infinitesimally small perturbation. In this case, the system of equations can be linearized, and it is sufficient to study a deformation δh of the epithelium–stroma interface translationally invariant in the y-direction with the complex form $\delta h(x,t) = \delta h_0 \exp (\omega t + {\mathrm { i}} q x)$ . Similarly, a perturbation of the apical surface of the epithelium can be written as $\delta H(x,t) = \delta H_0 \exp (\omega t + {\mathrm { i}} q x)$ . The stability of the system is then given by the dispersion relations ω(q) for each of the relaxation modes.

**Figure 1.** Schematic representation of the geometry considered in this study. An epithelium of thickness H is adjacent to a stroma of thickness L, itself anchored to a static, rigid structure at its opposite side. Perturbations of the epithelium–stroma interface and of the apical surface of the epithelium have infinitesimally small amplitudes, labelled δh and δH, respectively.
Download figure:
Standard image High-resolution image

2.2. Description of the models

In general, for length scales and timescales large compared with cell size and characteristic times of individual cellular processes, biological tissues can be described as continuous media with a potentially complex rheology, intermediate between those of liquids and solids [44–47]. In our previous work [14], as in other theoretical studies of tissue growth [48, 49], a viscous rheology was chosen to describe the epithelium, based on the assumption that cell–cell junction rearrangements lead to the relaxation of static stresses on long timescales [50, 51] and to effective cell flows [52]. This choice is motivated by the observation that relaxation times are of the order of tens of minutes to several hours for embryonic tissues in compression-plate experiments [44, 46] as well as in pipette-aspiration experiments [53]. Such results have also been proven accurate in the specific case of carcinomas, which show an almost complete stress relaxation under imposed sequential strains [54]. In addition, this choice is supported by the experimentally observed presence of surface tension at tissue boundaries [47, 53, 55, 56], which may drive cell sorting and certain morphogenetic movements [55, 57, 58]. On timescales much longer than the characteristic time of cell turnover in the tissue, repeated cycles of cell division and apoptosis can lead to the same result even for tissues that behave like elastic media or yield-stress fluids on short timescales [59, 60]. We therefore model the epithelium as a viscous fluid of shear viscosity η, which we take to be incompressible for the sake of simplicity, but with a rate of material production k_p corresponding to cellular duplication minus cellular death. In models 1–3, this rate is pre-imposed and takes the form of a single exponential function with a characteristic decay length l, whose amplitude is determined by a rate parameter k (detailed equations can be found in the supplementary information, available from stacks.iop.org/NJP/15/065011/mmedia, section 1).

In contrast to the epithelium where cells are confluent, the stroma is made up of a network of proteoglycans, collagen and elastin fibres, within which cells—mostly fibroblasts—are sparse [15, 61]. Fibroblasts constantly remodel the stromal fibres, which can therefore elastically resist deformation only on short and intermediate timescales but eventually reorganize and follow imposed deformations on long timescales. In addition, matrix metalloproteinases expressed either by the advancing tumour cells or directly by tumour-associated stromal cells can enhance this process by digesting the filaments of the stroma [1, 3, 7, 9]. The stroma should therefore be thought of more as a viscoelastic material than a purely elastic or viscous medium, with a characteristic timescale above which it reorganizes and flows. Here, three different descriptions are envisaged, namely those of an elastic solid with shear modulus μ (model 1), a Newtonian viscous fluid with shear viscosity η_s (model 2) or a viscoelastic material described by a Maxwell model with a single relaxation time τ = η_s/μ (model 3). In each case, the medium is supposed to the incompressible for the sake of simplicity.

The boundary conditions are as follows. At the apical surface of the epithelium, the total stress in contact with the lumen vanishes and the cell velocity is equal to the time derivative of the apical-surface location. Taking into account the epithelium apical surface tension γ_a, the normal component of the stress is given by the Laplace pressure, and its tangential component vanishes. At the epithelium–stroma interface, the discontinuity of the normal component of the stress is given by Laplace's law with interfacial tension γ_i, and the tangential components of the stress are continuous and equal to a finite surface-friction term with coefficient ξ. The normal component of the velocity is continuous and the stroma displacement in the z-direction is equal to δh. At the bottom of the stroma, the displacement vanishes. All the corresponding equations are detailed in section 1 of the supplementary information.

3. Results for model 3

We present here the first generalization of our previous study [14], namely that of a viscoelastic stroma with the cell-production rate function in the epithelium pre-defined as a function of distance from the epithelium–stroma interface. This model 3 encompasses both of the previously exposed models 1 and 2 as limit cases, where the stroma was either purely elastic or purely viscous [14] (see also the supplementary information, sections 2 and 3). We shall see the consequences of this more general model and how we recover the two cases studied before in particular regimes.

3.1. Mode structure and comparison with models 1 and 2

We obtain four different relaxation modes. This can be understood from the structure of the equations describing the boundary conditions, where time derivatives appear four times: once in the continuity condition for the velocity at the apical surface of the epithelium, and three times in the boundary conditions at the epithelium–stroma interface (supplementary information, equations (14) and (17)). We obtain four relaxation modes rather than three as in the case of an elastic stroma because of the additional relaxation timescale of the viscoelastic rheology here.

We present in figure 2 the structure of the relaxation modes as a function of the wavenumber q, and compare the present results with those of the purely elastic and purely viscous models 1 and 2, using the same set of parameters. We choose the two tissue viscosities to be equal to 10 MPa s for models 2 and 3, and fix the viscoelastic relaxation rate τ⁻¹ to 0.86 per day for model 3. This corresponds to an elastic modulus μ of 100 Pa for models 1 and 3. The structure of the figure is as follows. Plots corresponding to a first set of parameters are shown in figures 2(A)–(C), and for a second set in figures 2(D)–(F). For each of these parameter sets, the first column corresponds to the viscoelastic model 3, the second column to the elastic model 1 and the third one to the viscous model 2. In addition, for each of the mode structures computed, the corresponding plots are displayed in a range of wavevector values spanning 0–45 mm⁻¹ in the upper subpanels (indices 1), and 0–180 mm⁻¹ in the lower subpanels (indices 2). The range of omega values in the vertical axes is adapted to each individual case to show different relevant parts of the mode structure.

In the first set of graphs, we can see from the curves associated with model 3 in figure 2(A) that one of the modes relaxes towards the inverse viscoelastic relaxation time τ⁻¹ of 0.86 per day. In figure 2(B), the associated plots for the elastic model 1 almost identically reproduce the two lower modes, which are associated with timescales shorter than τ and therefore correspond to the elastic limit. We also recognize features of the two upper modes of the full viscoelastic model in both plots associated with models 1 and 2 (figures 2(B) and (C), respectively), although here these features are shared between the two different models depending on the range of wavenumbers. In the second set of graphs (figures 2(D)–(F)), we can appreciate the similarities between the viscoelastic and elastic relaxation modes whenever these are associated with short timescales, and between the viscoelastic and viscous relaxation modes in the other limit.

3.2. Asymptotic behaviours of the modes

It is instructive to look at the analytic expansions of the different modes in the respective limits of large and small wavenumbers q. In the limit of large wavenumbers, the modes associated with the epithelium–stroma interface decouple from those associated with the apical surface of the epithelium, since their characteristic decay lengths are of order q⁻¹, which is much smaller than H in this limit. Their asymptotic expressions up to constant order in a development in powers of q⁻¹ read

$\begin{eqnarray} \omega_1 &\simeq& -\frac{\gamma_{\rm i}}{2\eta}q-\frac{\mu}{\eta}+k-k_0,\nonumber\\ \omega_2 &\simeq& -\frac{\gamma_{\rm a}}{2\eta}q+k{\rm e}^{-H/l}-k_0,\nonumber\\\\ \omega_3 &\simeq& -2 \frac{\mu}{\xi}q-\frac{(\eta+ \eta_{\rm s})\mu}{\eta\eta_{\rm s}},\nonumber\\ \omega_4 &\simeq& 0-\frac{\mu}{\eta_{\rm s}}.\nonumber \end{eqnarray} \tag{ 1 }$

Among these, the first three expressions correspond to short characteristic times. We therefore expect their behaviours to be similar to those obtained in the case of an elastic stroma of modulus μ as already discussed in [14] and summarized in the supplementary information (available from stacks.iop.org/NJP/15/065011/mmedia), see equation (18). This is indeed the case for all of these three modes to leading order, and for the two first ones even up to constant order. In the third mode, the term of constant order mixes the epithelium and stroma viscosities, and therefore departs from the simpler μ/η term present in the elastic case. The additional fourth mode converges towards the inverse relaxation time τ⁻¹ = μ/η_s of the viscoelastic stroma, as illustrated in figure 2(A2).

In the limit of small wavenumbers, systematic expansions to leading order read

$\begin{eqnarray} \omega_1 &\simeq& -\frac{36\mu}{\eta}\frac{1}{H^3L^3}\frac{1}{q^6},\nonumber\\ \omega_2 &\simeq& -\frac{\mu}{4\eta}\frac{1}{HL}\frac{1}{q^2},\nonumber\\\\ \omega_3 &\simeq& k\,{\rm e}^{-H/l}-k_0,\nonumber\\ \omega_4 &\simeq& -\frac{L^2\left[3H\gamma_{\rm a}+2L(\gamma_{\rm a}+\gamma_{\rm i})\right]}{6\eta_{\rm s}}q^4.\nonumber \end{eqnarray} \tag{ 2 }$

Here, the expressions mix contributions coming from the entire system. The asymptotic expressions of the two first modes are identical to those already obtained in the case of an elastic stroma, and the one of the fourth mode corresponds to that obtained in the case of a viscous stroma (equations (19) and (21), supplementary information). This can be easily understood from the observation that the two first limits correspond to short-time dynamics and the fourth one to long-time dynamics. As for the third asymptotic expression, it is common to the three models. This is because this limit comes from pure mass conservation in the epithelium and is therefore independent of the stroma rheology. One can indeed recover this limit by integrating the continuity equation at q = 0 over the height of the epithelium to leading order in the perturbations. Expressions to next-to-leading order can be found in the supplementary information, section 5.

4. Nutrient-diffusion dynamics

As discussed in the introduction, the first three models exposed above and in the supplementary information (available from stacks.iop.org/NJP/15/065011/mmedia) sections 1–3 rely on a specific pre-defined function for the rate of cell division minus apoptosis in the epithelium, characterized by an exponential decay with a single characteristic length l (equations (2) and (3), supplementary information). This description is motivated by the fact that the supply of nutrients and growth factors in the epithelium comes from diffusion from the stroma. However, the function used above was set a priori, and to be more accurate, nutrient diffusion needs to be solved explicitly. Here, we solve the entire system of equations, where nutrients are produced at a given location in the stroma, diffuse in both tissues with tissue-specific diffusion constants, and can leak from the epithelium at its apical surface if they are in excess. They are consumed by the epithelial cells and influence their division rate.

4.1. Bulk equations

The equations describing the epithelium are still given by equation (1) of the supplementary information, but we now couple the cell-production rate k_p to the nutrient density ρ. In the absence of further detailed knowledge about this coupling, we assume a linear dependence of k_p on ρ:

$\begin{equation} k_{\rm p} = \kappa_1\rho-\kappa_0, \end{equation} \tag{ 3 }$

where κ₁ and κ₀ are effective phenomenological coefficients. These coefficients take positive values, since we expect the cell-division rate to increase as a function of ρ, as well as the overall cell population to starve and progressively die in the absence of nutrients. The nutrient-density function ρ is determined by the following diffusion-consumption equations. In the epithelium, nutrients diffuse with a coefficient D and are consumed by the cells with a rate c:

$\begin{equation} \partial_t\rho=D\Delta\rho-c\rho. \end{equation} \tag{ 4 }$

In the stroma, nutrients have the density ρ_s and diffuse with a coefficient D_s, without being consumed:

$\begin{equation} \partial_t\rho_{\rm s}=D_{\rm s}\Delta\rho_{\rm s}. \end{equation} \tag{ 5 }$

For the stroma, we again consider three different versions of the rheology as investigated above: an elastic rheology for model 4, a viscous one for model 5 and a viscoelastic one for model 6. The associated equations are presented in the supplementary information, section 1.

4.2. Boundary conditions

The mechanical boundary conditions are identical to those presented in the previous sections, but we need to specify the boundary conditions for the nutrient fields ρ and ρ_s. We assume a fixed concentration of nutrients $\bar {\rho }_0$ a distance d away from the epithelium–stroma interface in the stroma compartment:

$\begin{equation} {\rho_{\rm s}}_{|z=L-d}= \bar{\rho}_0. \end{equation} \tag{ 6 }$

At the epithelium–stroma interface, the density as well as the flux of nutrients in the direction perpendicular to the local interface are continuous:

$\begin{eqnarray} &&\rho_{|\rm interface} = {\rho_{\rm s}}_{|\rm interface},\nonumber\\\\ &&D\partial_\perp\rho_{|\rm interface} = D_{\rm s} \partial_\perp{\rho_{\rm s}}_{|\rm interface}.\nonumber \end{eqnarray} \tag{ 7 }$

Here '∂_⊥' stands for the partial derivative in the direction perpendicular to the local interface, oriented positively from the stroma towards the epithelium. At the apical surface of the epithelium, we allow for a potential leakage of nutrients from the epithelium into the lumen. Since we expect this leakage to increase with the local nutrient concentration, we write to linear order that the nutrient flux is proportional to this concentration locally:

$\begin{equation} -D\partial_\perp\rho_{|\rm apical} = v_{\rm off}\rho_{|\rm apical}. \end{equation} \tag{ 8 }$

5. Results for an elastic stroma: model 4

We present here the mode structure for model 4, where the stroma is treated as an elastic, incompressible material. The expressions of the steady-state solutions for the nutrient concentration, velocity and stress fields in the epithelium as well as the perturbed equations to linear order can be found in the supplementary information (available from stacks.iop.org/NJP/15/065011/mmedia), section 4. To solve this system of equations, we consider the case where nutrient diffusion is much faster than the characteristic relaxation of the system, both in the stroma (ω ≪ D_sq²) and in the epithelium (ω ≪ Dq² + c). This approximation is valid sufficiently close to the instability threshold and is a standard approximation of the treatment of diffusion-limited interface dynamics [16, 17]. In this regime, we can solve the equations describing nutrient diffusion at steady state first, and then substitute the obtained solutions into the mechanical equations.

5.1. Mode structure and influence of diffusion

The number of modes that we obtain is identical to the one obtained with a pre-imposed function for the production of cells, since the number of modes is prescribed by the structure of the mechanical boundary conditions. In the particular case of an elastic stroma, we obtain three relaxation modes, which can be seen from the expressions of the boundary conditions, where time derivatives appear three times (equations (14) and (15), supplementary information). In figure 3, we study the influence of nutrient diffusion on the structure of these modes. In figure 3(A), we show all three relaxation modes. The main qualitative difference between the present curves and those associated with model 1 (figure 1 in the supplementary information) is the presence of a clear maximum at low values of the wavenumber q (here around 5 mm⁻¹), that is at long wavelengths, in addition to the other broader maximum at smaller wavelength already present in the previous model (here around q = 30 mm⁻¹). To test whether this new maximum is indeed controlled by nutrient diffusion, we vary in figures 3(B) and (C) the values of the diffusion coefficients of nutrients, respectively, in the stroma and in the epithelium. In doing so, we change the value of the nutrient concentration at its production location ( $\bar {\rho }_0$ ) in order to keep a constant concentration at the epithelium–stroma interface in the unperturbed steady state. This allows us to decouple the influence of diffusion from that of the overall amount of nutrient supply. We see in figure 3(B) that increasing nutrient diffusion in the stroma makes the maximum at low q disappear, and that decreasing it instead makes the curve saturate towards a maximum limit curve with a pronounced peak. This behaviour is in qualitative agreement with that of an instability controlled by diffusion limitation [16, 17, 19]. The situation is different in terms of the variation of D, the diffusion coefficient in the epithelium, as illustrated in figure 3(C). In this case indeed, increasing D favours the instability and decreasing it makes the system stable. This stems from the fact that slowing down diffusion in the epithelium makes the layer of dividing cells very thin. This in turns lowers the driving force for the instability, which originates from differential cell flows in the epithelium due to inhomogeneous cell divisions.

**Figure 3.** Relaxation modes for model 4 as a function of nutrient diffusion. (A) Parameters are as follows: η = 30 MPa s, μ = 50 Pa, γ_i = 10 mN m⁻¹, γ_a = 1 mN m⁻¹, ξ = 10 GPa s m⁻¹, H = 300 μm, L = 600 μm, d = 75 μm, κ₁ = 100 × 10⁻⁶ m³ s⁻¹, D = 200 × 10⁻¹² m² s⁻¹, D_s = 100 × 10⁻¹² m² s⁻¹, c = 0.1 s⁻¹, $\bar {\rho }_0=1\,{\mathrm { m}}^{-3}$ and v_off = 0 m s⁻¹. All three relaxation modes are shown, of which only one becomes unstable and presents a clear maximum at low-q values. (B) The most unstable mode is shown for the same parameter set as in panel (A) for the default curve (1), and D_s is varied in the other curves. It takes the following values: (1) 100 × 10⁻¹² m² s⁻¹, (2) 200 × 10⁻¹² m² s⁻¹, (3) 1000 × 10⁻¹² m² s⁻¹, (4) 20 × 10⁻¹² m² s⁻¹, (5) 2 × 10⁻¹² m² s⁻¹. (C) A similar analysis is presented while varying D, which takes the following values: (1) 200 × 10⁻¹² m² s⁻¹, (2) 500 × 10⁻¹² m² s⁻¹, (3) 100 × 10⁻¹² m² s⁻¹, (4) 10 × 10⁻¹² m² s⁻¹. In both panels, $\bar {\rho }_0$ is changed as D_s or D are varied to keep the concentration of nutrients unchanged at the epithelium–stroma interface. Plots are coded in both colour and linestyles as indicated on the right of the figure.
Download figure:
Standard image High-resolution image

**Figure 3.** Relaxation modes for model 4 as a function of nutrient diffusion. (A) Parameters are as follows: η = 30 MPa s, μ = 50 Pa, γ_i = 10 mN m⁻¹, γ_a = 1 mN m⁻¹, ξ = 10 GPa s m⁻¹, H = 300 μm, L = 600 μm, d = 75 μm, κ₁ = 100 × 10⁻⁶ m³ s⁻¹, D = 200 × 10⁻¹² m² s⁻¹, D_s = 100 × 10⁻¹² m² s⁻¹, c = 0.1 s⁻¹, $\bar {\rho }_0=1\,{\mathrm { m}}^{-3}$ and v_off = 0 m s⁻¹. All three relaxation modes are shown, of which only one becomes unstable and presents a clear maximum at low-q values. (B) The most unstable mode is shown for the same parameter set as in panel (A) for the default curve (1), and D_s is varied in the other curves. It takes the following values: (1) 100 × 10⁻¹² m² s⁻¹, (2) 200 × 10⁻¹² m² s⁻¹, (3) 1000 × 10⁻¹² m² s⁻¹, (4) 20 × 10⁻¹² m² s⁻¹, (5) 2 × 10⁻¹² m² s⁻¹. (C) A similar analysis is presented while varying D, which takes the following values: (1) 200 × 10⁻¹² m² s⁻¹, (2) 500 × 10⁻¹² m² s⁻¹, (3) 100 × 10⁻¹² m² s⁻¹, (4) 10 × 10⁻¹² m² s⁻¹. In both panels, $\bar {\rho }_0$ is changed as D_s or D are varied to keep the concentration of nutrients unchanged at the epithelium–stroma interface. Plots are coded in both colour and linestyles as indicated on the right of the figure.
Download figure:
Standard image High-resolution image

We now compare the results of the two models more directly, namely those obtained with a pre-defined cell-production function (model 1) versus those obtained while coupling the cell-production function to nutrient diffusion (model 4). In order to do so, we derive from model 4 a cell-production rate as a function of z at steady state. We then fit the cell-production rate of model 1 to that function to generate the associated relaxation modes. As illustrated in figure 4 of the supplementary information, the fitting procedure gives a nearly perfect agreement between the two rate functions. We can therefore assume that the cell-production profiles are identical in the unperturbed situation for both models, and that all the differences that are observed are due to the effect of being coupled to nutrient diffusion to first order in the case of model 4 versus keeping the cell production unchanged in the case of model 1. We can see in figure 4 that coupling to nutrient diffusion overall increases the instability and introduces an additional maximum for the most unstable mode at a low q as commented on above. The remaining mode structure is unchanged and the stable modes are nearly identical in both models. In particular, in both the large- and low-q limits, we have the same asymptotic behaviours for the three modes with both models, as we now discuss in section 5.2.

**Figure 4.** Comparison of the modes obtained with model 4 (linestyle (1)) and with model 1 (linestyle (2)), following the procedure described in the text. The parameters are those of figure 3(A), together with the fitting parameters k = 2 divisions per day and l = 44 μm for model 1. Panels (A) and (B) present the same modes but over two different wavevector ranges.
Download figure:
Standard image High-resolution image

5.2. Asymptotic behaviours of the modes

To leading order, the large-q expansions of the three modes are identical to those obtained for model 1 and given by equation (18) of the supplementary information. This is also the case in the small-q limit for the two diverging modes of model 1, given by equation (19) of the supplementary information. This stems from the fact that these regimes correspond to a fast dynamics, where pure mechanics is at play without feeling the influence of the relatively slow cell-division events. However, the constant-order terms of these expansions do depend on the particular form of the cell-production function. The exact expressions have been derived but are very long and do not yield any particular physical insight in their full extent. They are therefore not presented in this paper. In figure 4 of the previous paragraph, the fitting procedure ensures that the limit of the upper mode is the same for both versions of the model.

6. Discussion for a viscous stroma: model 5

We now investigate the equivalent model for a viscous stroma.

6.1. Effect of nutrient diffusion on mode structure

As for model 4, we study in figure 5 the influence of the dynamics of nutrient diffusion. Here, we have two relaxation modes, since as compared with the elastic case, we lose the mode that originates from the tangential stress-balance condition at the epithelium–stroma interface (see equation (16), supplementary information). Similarly to the previous case, the main difference with model 2 is the presence of an extra maximum of the unstable mode at long wavelengths (figure 5(A) here and figure 2(A) of the supplementary information). In figure 5(B), we study the behaviour of this maximum as we vary the nutrient diffusion constant in the stroma D_s while keeping the nutrient concentration constant at the epithelium–stroma interface in the unperturbed steady state. We see that, similarly to what happened for model 4, decreasing D_s enhances the instability at long wavelengths which saturates as D_s goes towards zero. Increasing D_s instead tends to eliminate this maximum. In figure 5(C), we study the effect of the diffusion coefficient D in the epithelium compartment, following the same procedure. Here, increasing D first favours the instability by increasing the thickness of the dividing region, but eventually suppresses the instability at short wavelengths. This is because increasing D eventually renders cell division homogeneous over short lengthscales. Decreasing D instead tends to stabilize the system, which stems from the fact that the thickness of the dividing region is decreased as nutrient penetration is decreased.

**Figure 5.** Relaxation modes for model 5 as a function of nutrient diffusion. (A) Parameters are as follows: η = 10 MPa s, η_s = 10 kPa s, γ_i = 1 mN m⁻¹, γ_a = 1 mN m⁻¹, ξ = 10 GPa s m⁻¹, H = L = 300 μm, d = 100 μm, κ₁ = 100 × 10⁻⁶ m³ s⁻¹, D = 100 × 10⁻¹² m² s⁻¹, D_s = 200 × 10⁻¹² m² s⁻¹, c = 0.2 s⁻¹, $\bar {\rho }_0=0.3\,{\mathrm { m}}^{-3}$ and v_off = 0 m s⁻¹. Both relaxation modes are shown, of which only one becomes unstable and presents a clear maximum at low-q values. (B) The most unstable mode is shown for the same parameter set as in panel (A) for curve (1), and D_s is varied in the other curves. As for figure 3, the value of $\bar {\rho }_0$ is set accordingly to keep the concentration of nutrients unchanged at the epithelium–stroma interface. Values of D_s are (1) 200 × 10⁻¹² m² s⁻¹, (2) 20 × 10⁻¹² m² s⁻¹, (3) 2 × 10⁻¹² m² s⁻¹, (4) 400 × 10⁻¹² m² s⁻¹ and (5) 2000 × 10⁻¹² m² s⁻¹. (C) The same procedure is applied to study the variation of D, which takes the values (1) 100 × 10⁻¹² m² s⁻¹, (2) 1000 × 10⁻¹² m² s⁻¹, (3) 100 × 10⁻⁹ m² s⁻¹, (4) 20 × 10⁻¹² m² s⁻¹ and (5) 5 × 10⁻¹² m² s⁻¹.
Download figure:
Standard image High-resolution image

As in the case of an elastic stroma, we compare in figure 6 the relaxation modes obtained with model 5 and those obtained with model 2 where the cell-production function is obtained by the fitting procedure described above. As for figure 4, we see that the instability is globally enhanced, and especially at long wavelengths where a second maximum appears.

6.2. Asymptotic behaviours of the modes

In terms of the asymptotic expressions of the modes in the limit of small and large wavenumbers, we find as for models 1 and 4 that the expressions associated with models 2 and 5 match to leading order in the short-wavelength limit and correspond to those given by equation (20) of the supplementary information (available from stacks.iop.org/NJP/15/065011/mmedia). Beyond that, the expressions are complicated and depend on nutrient coupling. We illustrate in equation (35) of the supplementary information the asymptotic expression of one of these modes in a particular case to see how it differs from its equivalent for model 2 (equation (21), supplementary information).

7. Discussion for a viscoelastic stroma: model 6

We now consider the full model with coupling of the mechanical equations to the reaction–diffusion dynamics of nutrients, and with the full viscoelastic rheology for the stroma.

7.1. Parameter-variation study

In figure 7, we study the influence of some of the total 12 independent parameters on which the model depends. We choose here to present the dependence of the most unstable mode while independently varying six out of these parameters, which gives us an almost complete picture of the mode structure. The dependence in the diffusion constants D and D_s has been studied in detail for models 4 and 5 and is therefore not reproduced here for model 6. In figures 7(A) and (B), we show the mode structure for a default parameter set, which is then used in the other panels as a starting point for a systematic parameter-variation study. Only three out of the four modes are visible on the panels, since one relaxation rate takes very large absolute values and stands out of the plot. One of the four existing modes presents a region of instability. We note the presence of two distinct maxima in the unstable region, one around q = 5 mm⁻¹ and one around q = 20 mm⁻¹, corresponding to the two distinct maxima already discussed in the case of models 4 and 5. The viscoelastic relaxation rate τ⁻¹ is equal to 0.86 per day here. We recognize features of the elastic model for the lower mode, which corresponds to relaxation rates larger than τ⁻¹, and features of the viscous model for the two upper modes, which correspond to relatively long relaxation times. In figure 7(B), we display the same mode structure as in figure 7(A) but over a larger wavenumber domain, in order to visualize the convergence of the upper mode to τ⁻¹ in the short-wavelength limit.

**Figure 7.** Study of the relaxation-mode structure for model 6, as a function of several parameters entering the model. (A) and (B) The whole mode structure is shown for a default set of parameter values over two different wavevector ranges. Parameter values are η = 10 MPa s, η_s = 1 MPa s, μ = 10 Pa, γ_i = 5 mN m⁻¹, γ_a = 1 mN m⁻¹, ξ = 10 GPa s m⁻¹, H = L = 300 μm, d = 75 μm, κ₁ = 100 × 10⁻⁶ m³ s⁻¹, D = 100 × 10⁻¹² m² s⁻¹, D_s = 100 × 10⁻¹² m² s⁻¹, c = 0.1 s⁻¹, $\bar {\rho }_0=0.6\,{\mathrm { m}}^{-3}$ and v_off = 0 m s⁻¹. Only three of the four modes are visible, since one relaxation rate takes very large negative values and stands out of scale. In each of the remaining panels, the most unstable mode is displayed for different parameter values, as one parameter at a time is varied as compared with the default parameter set of panels (A) and (B). The plain-red linestyle, associated with number (1), corresponds to the default parameter set. Other parameter values are as follows. (C) Variation of η: (2) η = 5 MPa s; (3) η = 20 MPa s. (D) Variation of η_s: (2) η_s = 20 MPa s, (3) η_s = 100 kPa s. (E) Variation of H: (2) H = 100 μm; (3) H = 50 μm; (4) H = 900 μm. (F) Variation of d: (2) d = 50 μm; (3) d = 150 μm. (G) Variation of κ₁: (2) κ₁ = 50 × 10⁻⁶ m³ s⁻¹; (3) κ₁ = 150 × 10⁻⁶ m³ s⁻¹. (H) Variation of γ_i: (2) γ_i = 10 mN m⁻¹; (3) γ_i = 2 mN m⁻¹.
Download figure:
Standard image High-resolution image

In the remaining panels, we study the dependence of the most unstable mode on different parameters. We first show the dependence on η, the epithelium viscosity in figure 7(C). As expected and discussed in [14], the instability is increased with increasing values of η, which stems from the fact that the instability relies on viscous shear stresses in the epithelium. We recover this dominant tendency here. However, decreasing η never really stabilizes the system as was the case in [14], and the interface always remains unstable on a small interval at sufficiently small wavenumbers. This can be attributed to the coupling to nutrient diffusion, similarly to what happens in the context of crystal growth, where the system is always unstable at arbitrarily long wavelengths [16, 17]. In figure 7(D), we study the dependence on the stroma viscosity η_s. Increasing η_s to very large values tends to eliminate the low-q maximum, but is not enough to render the system stable. It is interesting to note that, contrary to the case of the Saffman–Taylor instability [62, 63], which occurs when a fluid of lower viscosity displaces a more viscous one in a Hele–Shaw cell, the relative values of the two fluid viscosities here is not a criterion to change the stability of the system. This fact was also observed in the context of model 2 and illustrated in figure 2 of the supplementary information (available from stacks.iop.org/NJP/15/065011/mmedia), section 3. In figure 7(E), we study the dependence in H, the thickness of the epithelium layer. A sufficiently thin epithelium does not display any instability. This is because the underlying mechanism of the instability remains the presence of viscous shear stresses within the epithelium, even when nutrient diffusion is introduced. In figure 7(F), we then study the dependence on d, the distance between the epithelium–stroma interface and the source of nutrients in the stroma. Varying this distance has a pronounced effect on the growth rate of the instability. This is natural since this distance directly influences the distribution of nutrients within the stroma, and so dramatically influences the stability of the system. In figure 7(G), we study the dependence on κ₁, which describes the coupling of cell division to nutrient concentration (equation (3)). Decreasing κ₁ allows only the long-wavelength maximum to remain unstable, and increasing it allows the other maximum to reach larger values. Finally, figure 7(H) shows that the epithelium–stroma interfacial tension γ_i has a stabilizing effect, and primarily at short wavelengths.

Among the parameters whose variations have not been displayed in this figure, D and D_s have similar effects to those seen in figures 3 and 5 in the frameworks of models 4 and 5, and the remaining parameters have a simple, intuitive effect. Increasing the elastic modulus μ of the stroma tends to stabilize the system, and so does the increase of c, the rate of nutrient consumption by the epithelial cells. On the other hand, the stroma thickness L as well as the total amount of nutrient produced $\bar {\rho }_0$ have a destabilizing effect (graphs not shown).

7.2. Mode structure and asymptotic expressions

Similarly as before, the different associated expressions of the fast modes to leading order are unchanged as compared with the case of model 3. In the limit of small wavelengths, we retrieve the expressions given by equation (1) to linear order, which is explained for the three first modes by the fact that these correspond to fast relaxations of the system, and for the fourth mode because its limit corresponds to the inverse relaxation time of the viscoelastic material constituting the stroma, which is independent of cell division. In the long-wavelength regime, the expressions to leading order of the two first modes are identical to those obtained without nutrient diffusion and are given by equation (2). The situation is different for the modes of order constant and proportional to q⁴, for which the specific coefficients do depend on nutrient coupling. The general expressions are very long, but the example of equation (35) of the supplementary information (available from stacks.iop.org/NJP/15/065011/mmedia) discussed above also holds here.

8. Discussion

In this work, we have proposed a mechanism for the formation of undulations at the epithelium–stroma interface that arise from a physical instability intrinsic to the structure of a multilayered epithelium. This instability, presented originally in [14] and discussed here in more detail, is a hydrodynamic instability arising from viscous shear stress originating from flow within the tissue due to cell renewal. Such an instability can develop in a healthy epithelium, depending on its physical characteristics such as the thickness of the cell-division region, its long-term viscosity due to cell–cell adhesion and the mechanical resistance of the stroma. These results are in agreement with the observation that the degree of undulations in vivo typically depends on the grade of the dysplastic epithelium. The grade itself correlates with the number of dividing cell layers [12], and is part of the diagnosis of dysplastic tissues, for example of oral cancer [13, 64]. The mechanical resistance of the stromal tissue can also influence the morphology of the basement membrane. This has been observed for example to a mild extent in vivo in the stratified epithelia of the cornea and of the limbus [65], and in a clearer way in reconstituted epithelia in vitro [66]. There, different degrees of undulations were observed as a function of the properties of the supporting scaffold on which the tissue was grown. More generally, this instability could be present in all sufficiently viscous fluids with internal flow due to material production or destruction. Growth factors, proteases, external feedback and abnormal proliferation are not required for the occurrence of this instability.

Nevertheless, introducing coupling of cell division to the diffusion of nutrients, growth factors or oxygen from the stroma adds a destabilizing effect, similar to that occurring during diffusion-limited growth [18, 19]. This leads to a significant enhancement of the previously proposed instability and to the appearance of an additional maximum of the unstable mode at long wavelengths. Instabilities originating from similar coupling terms have been identified in living systems, for example in the case of bacterial-colony growths, where patterns similar to those associated with aggregation phenomena and viscous fingering have been observed [67–69]. Such a coupling can lead to fractal branching patterns via the process of diffusion-limited aggregation [70, 71] or other types of branching patterns depending on the bacterial morphotype [72, 73].

The destabilizing effect of the coupling to nutrient diffusion is reminiscent of the Mullins–Sekerka instability in the context of crystal growth [16, 17]. This analogy relies on the fact that nutrients here play a similar role as either latent heat when crystal growth occurs in an undercooled melt of a pure substance or solute molecules when it occurs in an isothermal but binary, supersaturated solution: outward-pointing epithelium protrusions have access to a larger amount of nutrients than retracted regions, and therefore proliferate faster, enhancing the already existing protrusion. The instability therefore develops faster when nutrient diffusion is slower, especially at long wavelengths as illustrated in figures 3 and 5. Nevertheless, differences exist between the diffusion-driven instability proposed here and that of growing crystals. For one, in the latter case, the solid phase grows by addition of new material coming from the environment to the interface. This phase therefore grows only at the surface and, once formed, remains rigid and static in the bulk. Here however, the new material comes from cell division within the epithelial tissue itself, and growth occurs as a bulk phenomenon. As a consequence, the instability proposed here can only develop for sufficiently thick epithelia, as illustrated in figure 7(E), and its driving force still relies on differential cell flows, as proposed originally [14].

The instability discussed in this work may also evoke the Saffman–Taylor instability, which occurs when a fluid of lower viscosity displaces a more viscous one in a Hele–Shaw cell [62, 63]. However, the two mechanisms are very different and should not be confused. In the case of the Saffman–Taylor instability, an external force in the form of a global pressure difference is required to drive the system out of equilibrium. The dynamics is governed by Darcy's law arising from the balance of pressure gradient with friction of the fluid with the background. The Laplace equation for the pressure field together with the constant-pressure boundary conditions determine the instability and the mode structure. On the other hand, in the case studied here, forces are generated within the epithelial tissue by internal processes. In addition, the field responsible for the second, long-wavelength maximum is the nutrient-concentration field, which follows a Laplace equation in the stromal compartment in the regime of fast diffusion, rather than the pressure field. A particular illustration of the difference between the two instabilities can be seen from figure 7 of the main text and figure 2 of the supplementary information (available from stacks.iop.org/NJP/15/065011/mmedia), which show that here, the instability can develop when either of the two fluids—epithelium or stroma—is the more viscous one.

Beyond the mechanisms discussed in this work, which rely entirely on cell-division driven cellular flows and their coupling to passive diffusion, there exist more specific, active, biological mechanisms that can cause significant aberrations of the tissue interface in malignant cases. As one example, neoplastic epithelial cells are known to excrete signalling molecules that recruit and stimulate stromal fibroblasts and induce the differentiation of monocytes into macrophages [1, 3]. This further induces the production and release of epithelial growth factors, matrix metalloproteinases and angiogenetic factors [1, 74, 75]. This effect may enhance the progression of protrusions in a positive feedback loop, since advancing protrusions increase the concentration of signalling molecules in regions of the stroma close to them, further stimulating the expression of growth factors by fibroblasts. As a second example, many malignant epithelial tissues and their stromal-associated cells are known to secrete matrix-degrading enzymes such as metalloproteinases [1, 3, 7, 9]. These enzymes play a central role in the invasion of epithelial protrusions in the stroma by cutting through the dense meshwork of collagen fibres, carving out space for the expanding tumour. Proteases such as matrix metalloproteinases 2 and 9 can be produced by the tumour-associated macrophages [9, 75], a mechanism that contributes to the cross-talk mentioned above.

While neither of these effects is directly included in our model, we expect the driving force behind the advancing protrusions to still originate from the fundamental mechanical picture proposed here. Indeed, these effects may enhance an already existing instability, but do not represent an alternative driving force to the present model. The situation is fundamentally different in the case where cells from the epithelium have undergone the transition to a motile phenotype and invade the stroma [4, 5]. Although clearly only relevant for malignant tissues, this mechanism is fundamentally different from the one discussed in this work and provides a true alternative process that can lead to pronounced irregular pathologies observed in vivo [10], and ultimately drive cancerous invasion. However, one may speculate that the acquisition of invasive motile phenotypes itself comes after an original protrusion has formed rather than be responsible for its initiation. In two-dimensional tissue migration for example, the emergence of fingering protrusions seems to arise spontaneously [76] and is often accompanied by the emergence of a different phenotype at the tip of the protrusion [77]. The appearance of such leader cells may however be the signature of a loss of contact inhibition [78–80], and may therefore be induced by the original perturbation rather than be responsible for its initiation. One may therefore speculate that a physical mechanism such as the instability proposed here could initiate the emergence of invasive phenotypes.

Acknowledgments

We thank Professors J Prost and J-F Joanny for thinking about the original ideas with us as well as for numerous constructive discussions. The authors acknowledge support from the European network MitoSys.

Morphological instabilities of stratified epithelia: a mechanical instability in tumour formation

Article metrics

Author e-mails

Author affiliations

Author notes

Dates

Abstract

1. Introduction