Mechanically driven interface propagation in biological tissues

Animal tissues are made of many different cell types, and interfaces between different cell populations are ubiquitous in the developing and adult organism [1]. An example is the interface that forms between a growing tumor and the healthy cells of the surrounding host tissue [2]. In this article, we study how such interfaces between two distinct cell populations evolve in the presence of cell division and cell death, focussing on the effects of mechanical coupling between the net cell division rates and pressure. It has been hypothesized that the rates of cell division and apoptosis (cell death) are not only dependent on the biochemical conditions present in the tissue but depend on mechanical forces exerted by the surrounding cells [3, 4]. This suggestion was recently underpinned by experimental findings that demonstrate the variation of the net cell division rate with an externally applied pressure [5]. Here, we develop a continuum theory of two cell populations separated by an interface. Our theory takes into account rates of cell division and death which can be influenced by local stress. Resulting cell flows and interface movements follow from the dynamics of cells coupled to force balance in the tissue, and we find moving front solutions that describe the invasion of one tissue into a second one.

To investigate the dynamics of tissues with two distinct cell populations, we develop a continuum description that accounts for cell flow fields and stress distributions on large scales. We assume that the tissue comprises cells of type A and B with number densities ${{n}_{A,B}}$ and average cell volumes ${{\Omega }_{A,B}}$, respectively, such that ${{n}_{A}}{{\Omega }_{A}}+{{n}_{B}}{{\Omega }_{B}}=1$. Cell number balance is expressed by

$$\begin{eqnarray}&&{{}_{t}}{{n}_{i}}+{{}_{\alpha }}\left( {{n}_{i}}v_{\alpha }^{i} \right)={{n}_{i}}{{k}_{i}},\end{eqnarray} \tag{ 1 }$$

where $v_{\alpha }^{i}$ is the velocity of the cell flow field of type i ($i=A,B$), and ${{k}_{i}}=k_{\text{d}}^{i}-k_{\text{a}}^{i}$ is the respective net cell division rate, i.e., the difference between the rates of cell division and apoptosis. The rates ${{k}_{i}}$ in general depend on pressure. Their variation is characterized by the homeostatic pressures $P_{i}^{\text{h}}$: the homeostatic pressure $P_{i}^{\text{h}}$ of a tissue is the value of the pressure P for which the tissue is in a homeostatic state in which cell division and apoptosis balance and no net growth occurs. Near the homeostatic state, one can expand the net cell division rate to linear order in pressure differences as ${{k}_{i}}\simeq {{\kappa }_{i}}(P_{i}^{\text{h}}-P)$. Here, ${{\kappa }_{i}}$ is a coefficient that in principle is different for the two cell types A and B. For simplicity, we choose the same value ${{\kappa }_{i}}=\kappa$ for the two tissues in the following and investigate the effects of a finite difference of the homeostatic pressures between the two tissues $\Delta {{P}^{\text{h}}}P_{A}^{\text{h}}-P_{B}^{\text{h}}>0$.

Some tissues are separated from each other by physical barriers such as the basal membrane [6]. However, tissue interfaces may also be smooth with possible mixing of cells of different types, see figure 1. On large scales, the shape of the interface is then described by the volume fraction $\phi {{n}_{A}}{{\Omega }_{A}}$ as a function of position. The average cell velocity ${{v}_{\alpha }}$ is given by

$$\begin{eqnarray}&&{{v}_{\alpha }}=\phi v_{\alpha }^{A}+(1-\phi )v_{\alpha }^{B},\end{eqnarray} \tag{ 2 }$$

and the relative flux ${{J}_{\alpha }}$ obeys

$$\begin{eqnarray}&&\phi v_{\alpha }^{A}=\phi {{v}_{\alpha }}+{{J}_{\alpha }},\quad (1-\phi )v_{\alpha }^{B}=(1-\phi ){{v}_{\alpha }}-{{J}_{\alpha }}.\end{eqnarray} \tag{ 3 }$$

Under the assumption that the cells are incompressible, i.e., ${{\Omega }_{A}}$ and ${{\Omega }_{B}}$ constant, the cell number balance equations for each cell type then combine to a condition on the average velocity ${{v}_{\alpha }}$ and a dynamic equation for ϕ,

$$\begin{eqnarray}&&{{}_{\alpha }}{{v}_{\alpha }}=\phi {{k}_{A}}+(1-\phi ){{k}_{B}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{}_{t}}\phi +{{v}_{\alpha }}{{}_{\alpha }}\phi =-{{}_{\alpha }}{{J}_{\alpha }}+\phi (1-\phi )({{k}_{A}}-{{k}_{B}}).\end{eqnarray} \tag{ 5 }$$

Equation (4) generalizes the the classical incompressibility condition on the divergence of the velocity in the presence of cell division.

Zoom In Zoom Out Reset image size

In a previous work, we have shown that cell division and apoptosis give rise to a diffusive motion of cells even in the absence of spontaneous cell motility [7]. The existence of such an effective diffusion requires the presence of noise in the tissue, for example due to cell shape fluctuations and to the stochasticity of cell division and apoptosis. Such processes can give rise to mutual diffusion between the two cell types that can be described by an effective diffusion coefficient D and generates a relative flux ${{J}_{\alpha }}=-D{{}_{\alpha }}\phi$. Note that the effective diffusion constant is not caused by thermal fluctuations.

Taking the dependence of cell proliferation on pressure into account, the dynamic equation for ϕ then becomes

$$\begin{eqnarray}&&{{}_{t}}\phi +{{v}_{\alpha }}{{}_{\alpha }}\phi =D_{\alpha }^{2}\phi +\phi (1-\phi )\kappa \Delta {{P}^{\text{h}}}.\end{eqnarray} \tag{ 6 }$$

Equation (6) is a generalized version of the Fisher–Kolmogorov equation which was initially put forward to describe the expansion dynamics of advantageous alleles in a population with shared genome [8, 9]. However, the dynamics of the cell volume fraction described by equation (6) contains an additional term that describes advection with velocity ${{v}_{\alpha }}$. This velocity is generated in the tissue by active processes and is related to the proliferation of the two cell populations via equation (4). This is a situation that is very different from cases where flows are created externally such as in studies of bacterial colony growth as e.g. in [10]. Since the net cell division rates depend on pressure, the Fisher population dynamics is eventually coupled to the tissue mechanics. This coupling is generic for tissues, where cells do not only divide and die but also interact mechanically. To our knowledge, diffusive growth has so far been discussed without mechanical interactions [11].

To solve for the velocity field ${{v}_{\alpha }}$, we need to consider the stress distributions and force balance in the tissue. In the absence of external forces, force balance is expressed by ${{}_{\beta }}{{\sigma }_{\alpha \beta }}=0$, where ${{\sigma }_{\alpha \beta }}$ denotes the stress tensor in the tissue. The stress tensor can be split into an isotropic and a traceless part as ${{\sigma }_{\alpha \beta }}=-P{{\delta }_{\alpha \beta }}+{{\overset{}{\mathop{\sigma }}\,}_{\alpha \beta }}$, where P is the total pressure and ${{\overset{}{\mathop{\sigma }}\,}_{\alpha \alpha }}=0$. For an incompressible tissue, P becomes a Lagrange multiplier to impose the generalized incompressibility constraint and from equation (4) we find

$$\begin{eqnarray}&&P=P_{B}^{\text{h}}+\phi \Delta {{P}^{\text{h}}}-\left( {{}_{\alpha }}{{v}_{\alpha }} \right)/\kappa .\end{eqnarray} \tag{ 7 }$$

In the hydrodynamic limit, the deviatoric stress ${{\overset{}{\mathop{\sigma }}\,}_{\alpha \beta }}$ is given by the constitutive equation

$$\begin{eqnarray}&&{{\overset{}{\mathop{\sigma }}\,}_{\alpha \beta }}=2\eta {{\overset{}{\mathop{v}}\,}_{\alpha \beta }}-B\left[ ({{}_{\alpha }}\phi )({{}_{\beta }}\phi )-\frac{1}{3}({{}_{\gamma }}\phi )({{}_{\gamma }}\phi ){{\delta }_{\alpha \beta }} \right],\end{eqnarray} \tag{ 8 }$$

where we assumed that the tissue behaves as a viscous fluid on long times [7, 12]. Here, η is the shear viscosity and ${{\overset{}{\mathop{v}}\,}_{\alpha \beta }}=\frac{1}{2}({{}_{\alpha }}{{v}_{\beta }}+{{}_{\beta }}{{v}_{\alpha }}-\frac{2}{3}{{}_{\gamma }}{{v}_{\gamma }}{{\delta }_{\alpha \beta }})$ is the traceless part of the velocity gradient tensor. The second term on the right-hand side is due to the coupling between stress and composition gradients and its form is dictated by symmetry. It corresponds to the Ericksen stress in systems near thermodynamic equilibrium, and it can be derived explicitly from a Landau–Ginzburg free energy depending on a scalar parameter ϕ [13]. If the width of the interface is of order l, the surface tension between A and B cells is of order $B/l$. It describes interfacial tensions as they arise for example from differential adhesion between the two cell types [12].

We are interested in the propagation of an interface between two cell populations that are at their respective homeostatic states and at rest at $x=\pm \infty$. We consider a thin, effectively two-dimensional tissue of fixed height h subject to friction with the underlying substrate [14], which breaks Galilean invariance. Note that for bulk, i.e. three-dimensional, tissues, the role of substrate friction could be played by friction with the interstitial fluid. This can be discussed in descriptions of tissues with permeation [15].

Force balance in the film reads

$$\begin{eqnarray}&&{{}_{\beta }}{{\bar{\sigma }}_{\alpha \beta }}=\gamma {{v}_{\alpha }}\end{eqnarray} \tag{ 9 }$$

($\alpha ,\beta =x,y$), where γ is a coefficient describing friction with the substrate and ${{\bar{\sigma }}_{\alpha \beta }}$ is the stress in the film averaged in the z-direction

$$\begin{eqnarray}&&{{\bar{\sigma }}_{\alpha \beta }}=\frac{1}{h}\mathop{\int }^{}_{0}^{h}{{\sigma }_{\alpha \beta }}\text{d}z.\end{eqnarray} \tag{ 10 }$$

For simplicity, we consider tissues for which ϕ is homogeneous along y and fields vary with x only. From the force balance equation (9), using the expressions for P and ${{\overset{}{\mathop{\sigma }}\,}_{\alpha \beta }}$, we obtain an equation for the velocity profile ${{v}_{x}}$,

$$\begin{eqnarray}&&\left( {{\kappa }^{-1}}+\frac{4}{3}\eta \right)_{x}^{2}{{v}_{x}}-\Delta {{P}^{\text{h}}}{{}_{x}}\phi -\frac{2}{3}B_{x}^{2}\phi {{}_{x}}\phi =\gamma {{v}_{x}}.\end{eqnarray} \tag{ 11 }$$

Equations (6) and (11) define the dynamics of interface propagation, together with the boundary conditions $\phi (-\infty )=1$, $\phi (+\infty )=0$, and ${{v}_{x}}(\pm \infty )=0$. These equations can be written in a dimensionless form,

$$\begin{eqnarray}&&{{}_{T}}\phi +V{{}_{X}}\phi =_{X}^{2}\phi +\phi (1-\phi ),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{{\Lambda }^{2}}_{X}^{2}\;V-V=2\Lambda {{V}_{0}}{{}_{X}}\phi \left( 1+\beta _{X}^{2}\phi \right).\end{eqnarray} \tag{ 13 }$$

Here, we have used the characteristic interface width ${{l}_{0}}$ and timescale ${{t}_{0}}$ with

$$\begin{eqnarray}&&{{l}_{0}}={{(D/(\Delta {{P}^{\text{h}}}\kappa ))}^{1/2}},\quad \quad \quad {{t}_{0}}=1/(\Delta {{P}^{\text{h}}}\kappa ),\end{eqnarray} \tag{ 14 }$$

to normalize lengths and times: $X=x/{{l}_{0}}$, $T=t/{{t}_{0}}$, and $V={{v}_{x}}{{t}_{0}}/{{l}_{0}}$. We have furthermore introduced the dimensionless parameters

$$\begin{eqnarray*}&&{{\Lambda }^{2}}=\left( 1+\frac{4}{3}\kappa \eta \right)\Delta {{P}^{\text{h}}}/(D\gamma ),\quad \quad \quad {{V}_{0}}=\Delta {{P}^{\text{h}}}/(2D\gamma \Lambda ),\quad \quad \quad \beta =\frac{2}{3}B\kappa /D.\end{eqnarray*}$$

For ${{V}_{0}}=0$, the convective velocity vanishes, V = 0. Starting with a steep initial condition $\phi (X)={{(1+{{e}^{X/\epsilon }})}^{-1}}$, with $\epsilon 1$, the system reaches for long times traveling-wave solutions of the form $\phi (X,T)={\Phi}(U)$, where $U=X-CT$ with wave speed $C={{C}_{0}}$ and interface profiles ${\Phi}(U)={\Phi_{0}}(U)$ where ${{C}_{0}}=2$ [16, 17]. This solution of the classical Fisher equation is a so-called pulled-front solution for which the wave speed is determined by the linearized dynamics in the tail of the profile. Pulled-front solutions cannot have larger wave speeds than C = 2. Here we ask whether traveling-wave solutions persist when ${{V}_{0}}>0$ and tissue mechanics and advection have to be taken into account. A simple limit occurs when $\Lambda 1$. In this case, equations (12) and (13) decouple because $V(U=0)\simeq {{V}_{0}}$ remains approximately constant along the interface where $\phi (X)$ varies. From this argument, we find for $\beta =0$ that $V{{}_{X}}\phi \simeq {{V}_{0}}{{}_{X}}\phi$ in equation (12). Thus, the interface profile is again given by ${\Phi_{0}}$, however with traveling speed $C={{C}_{0}}+{{V}_{0}}$. Finite interfacial tension due to $\beta >0$ does not change this result as long as $\beta \Lambda$. Such solutions are so-called pushed-front solutions for which the wave speed is determined by nonlinearities as discussed below [17].

We obtain numerical solutions to the dynamic equations (12) and (13) for different values of ${{V}_{0}}$, $\Lambda$, and β. Starting from a steep initial interface, the interface evolves to a stationary profile ${\Phi}(U)$ and moves with a speed C for all tested parameter values. Examples of traveling-wave profiles $\phi (U)$ and V(U) are shown in figure 2. For small ${{V}_{0}}$ as well as for $\Lambda 1$, the profiles are well captured by the approximate solution ${\Phi_{0}}$ to the original Fisher equation. Deviations start to be important for larger values of ${{V}_{0}}$ and finite interfacial-tension coefficient β, see figure 2(g): because of the term $\propto \ _{X}^{2}\phi$ in equation (13), the interfacial profile broadens, and wave speed increases. The dimensionless wave speed C is shown in figure 3 for different values of ${{V}_{0}}$, $\Lambda$, and β. Our numerical analysis confirms the simple Fisher-wave limit for $\Lambda 1$, in which the advective velocity ${{V}_{0}}$ adds up to the wave speed ${{C}_{0}}$. It is interesting to note the effect of finite β on the observed wave speeds. For $\beta =0$, the wave speed C is bounded from above by ${{C}_{0}}+{{V}_{0}}$, and C increases with increasing $\Lambda$ towards the predicted value ${{C}_{0}}+{{V}_{0}}$. For large enough β, the reverse can be observed: for finite $\Lambda$, the wave speed C actually exceeds ${{C}_{0}}+{{V}_{0}}$ and decreases with $\Lambda$ for fixed ${{V}_{0}}$. Taken together, these results highlight the role of the interplay between cellular diffusion, cell–substrate friction, and the mechanical coupling of cell proliferation to pressure.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One can distinguish between pulled fronts and pushed fronts [17]. The velocity of pulled fronts is set by the dynamics of the asymptotic tail which pulls the main profile behind. Because the tail decays to zero, its dynamics can be understood by a linearized theory. Pushed fronts are driven by the part of the profile where nonlinearities are essential. The profile pushes the linear tail forward at a speed that is larger than the pulled-front velocity. Both pulled- and pushed-front solutions can be obtained as profiles $\phi (U)$ and V(U) that solve equations (12) and (13). Such solutions exist for all $C>0$. The tails of the corresponding profiles are of the form $\phi (U)\simeq {{d}_{1}}{{e}^{-{{\lambda }_{1}}U}}+{{d}_{2}}{{e}^{-{{\lambda }_{2}}U}}$ which follows from linearizing equations (12) and (13). Here, ${{d}_{1}}$ and ${{d}_{2}}$ are coefficients and the (inverse) decay lengths are given by ${{\lambda }_{1,2}}=C/2\mp {{[{{(C/2)}^{2}}-1]}^{1/2}}$, with ${{\lambda }_{1}}{{\lambda }_{2}}$. For localized initial conditions, the system either reaches a pulled-front solution with $C={{C}_{0}}=2$ and ${{\lambda }_{1,2}}=1$ or a pushed-front solution for which $C>{{C}_{0}}$. In the latter case, the velocity C is set by the requirement that the coefficient ${{d}_{1}}=0$ must vanish such that the asymptotic behavior of the front is given by $\phi (U)\simeq {{d}_{2}}{{e}^{-{{\lambda }_{2}}U}}$, corresponding to the steeper front.

This analysis allows us to distinguish pulled-front solutions from pushed fronts in our numerical study. The inset of figure 3 shows a line separating two regions in which pulled and pushed fronts occur as a function of ${{V}_{0}}$ and $\Lambda$, in the case where $\beta =0$. For small ${{V}_{0}}$ front propagation is dominated by diffusion and the system develops a pulled front which moves at the same dimensionless speed $C={{C}_{0}}$ as a Fisher wave without advection. Beyond a critical value of ${{V}_{0}}$, the front is pushed by nonlinearities and moves at an increased speed C that depends on ${{V}_{0}}$ in this advection-dominated regime. Both regimes are separated by a sharp transition from pulled to pushed fronts shown as a solid line.

The case of a linear geometry discussed above has been chosen for conceptual clarity. For a small clone of A cells surrounded by B cells on a substrate, the expansion is more realistically captured by a growing circular domain of A cells in two dimensions. Using circular symmetry for simplicity, the dynamic equation for the cell volume fraction ϕ reads in polar coordinates

$$\begin{eqnarray}&&{{}_{t}}\phi +{{v}_{r}}{{}_{r}}\phi =D\left( _{r}^{2}\phi +\frac{{{}_{r}}\phi }{r} \right)+\phi (1-\phi )\kappa \Delta {{P}^{\text{h}}},\end{eqnarray} \tag{ 15 }$$

where r is the radial coordinate. The force balance reads

$$\begin{eqnarray}&&\left( {{\kappa }^{-1}}+\frac{4}{3}\eta \right)\left[ _{r}^{2}{{v}_{r}}+\frac{{{}_{r}}{{v}_{r}}}{r}-\frac{{{v}_{r}}}{{{r}^{2}}} \right]-\gamma {{v}_{r}}=\Delta {{P}^{\text{h}}}{{}_{r}}\phi +\frac{2}{3}B\left( _{r}^{2}\phi +\frac{{{}_{r}}\phi }{r} \right){{}_{r}}\phi .\end{eqnarray} \tag{ 16 }$$

The boundary conditions are given by $\phi \prime (0)=0$, $\phi (\infty )=0$, $v_{r}^{\prime }(0)=0$, and ${{v}_{r}}(\infty )=0$. Note that traveling waves of the form $\phi (r,t)={\Phi}(r-ct)$ cannot exist for this system because the dynamic equations are not invariant with respect to changes of r. In the limit of long times, however, the radius of the interface becomes large, ${{r}_{\text{i}}}\max ({{l}_{0}},\Lambda {{l}_{0}})$, and one approaches the one-dimensional case described above. For smaller ${{r}_{\text{i}}}$, the finite curvature of the interface has to be taken into account.

We obtain numerical solutions to the dynamic equations, starting from localized initial conditions. The resulting behavior is shown in figure 4, where an initially localized distribution of A cells first decreases in magnitude and widens before it increases its magnitude and grows out radially. Interestingly, the number of A cells first decreases, because of the increased pressure due to the interfacial tension, which counteracts the proliferative advantage of the stronger A cells (see inset in figure 4). The time dependence of the total number of A cells ${{N}_{A}}=2\pi \Omega _{A}^{-1}\mathop{\int }^{}_{0}^{\infty }\phi r\text{d}r$, where we take ${{\Omega }_{A}}$ to be an effective two-dimensional cell volume, can be discussed using equation (15), which implies

$$\begin{eqnarray}&&\frac{\text{d}{{N}_{A}}}{\text{d}t}=\frac{2\pi }{{{\Omega }_{A}}}\mathop{\int }^{}_{0}^{\infty }[-{{v}_{r}}{{}_{r}}\phi +\kappa \Delta {{P}^{\text{h}}}\phi (1-\phi )]r\text{d}r\end{eqnarray} \tag{ 17 }$$

(note that the diffusive term disappears under the integral). With equation (7), one eventually obtains

$$\begin{eqnarray}&&\frac{\text{d}{{N}_{A}}}{\text{d}t}=\frac{2\pi }{{{\Omega }_{A}}}\mathop{\int }^{}_{0}^{\infty }\kappa \phi (P_{A}^{\text{h}}-P)r\text{d}r.\end{eqnarray} \tag{ 18 }$$

Zoom In Zoom Out Reset image size

This equation shows that the number of A cells drops as long as $P>P_{A}^{\text{h}}$. If cell motion can be ignored, the typical pressure inside the region of A cells ${{P}_{A}}\simeq P_{B}^{\text{h}}+\Gamma /R$ has a contribution from the Laplace pressure $\Gamma /R$, where $\Gamma \simeq B/l$ is the surface tension (or line tension in two dimensions), l is the interface width, and R is the radius of the A-cell-rich region. We thus have

$$\begin{eqnarray}&&\frac{\text{d}{{N}_{A}}}{\text{d}t}\simeq {{N}_{A}}\kappa (\Delta {{P}^{\text{h}}}-\Gamma /R).\end{eqnarray} \tag{ 19 }$$

The suppression of A cells occurs if R is below a critical radius ${{R}_{\text{c}}}\approx B/(l\Delta {{P}^{\text{h}}})$.

At long times, however, diffusion smears out the interface profile and l first increases as ${{l}^{2}}\simeq Dt$ until it saturates at the finite value ${{l}_{0}}$ defined in equation (14). Therefore the surface tension drops and the critical radius becomes smaller. The stronger A cells with higher homeostatic pressure ${{P}^{\text{h}}}$ can therefore possibly escape suppression and invade the weaker tissue with lower ${{P}^{\text{h}}}$ even if initially the A-cell region is smaller than the critical radius.

This long-time growth of A cells into a region of B cells occurs only if the total number of A cells always remains larger than one. Note however that our continuum description with continuously varying cell volume fraction does not take the discrete character of cells into account. Therefore, this description misses the absorbing state when the last remaining A cell undergoes apoptosis. In this case, the continuum description breaks down and a more refined discrete description is needed.

If we apply this result to the situation of cancer invasion, these arguments therefore suggest that cell diffusion could increase significantly malignancy: when cancer cells become less cohesive and mix with cells of the host tissue, the tumor cannot be suppressed by interface tension only. If tumor and host tissue cannot freely mix (possibly separated by a physical barrier such as the basal membrane), the tumor can be suppressed if smaller than a critical radius ${{R}_{\text{c}}}$ [4] and thus is statistically less dangerous.

The interface dynamics discussed in this article is generic for tissues in which cells undergo cell division and apoptosis and interact. This is the case for developing tissues, where cell competition has been observed for example in the growing larva of Drosophila [18]; another example is the interface between a growing tumor and its surrounding host tissue. Using a continuum description, we find two regimes of interface propagation: a diffusive regime in which relative fluxes dominate the expansion (pulled front), and a propulsive regime in which the proliferation gives rise to convective flows (pushed front). In general, experiments will distinguish the regime in which a given system operates. Recent experiments focused on the growth dynamics of tumor cell spheroids that consist of one cell type only [5]. Taking a closer look at the combined growth dynamics of two cell populations will hopefully provide more insight into the different modes of homeostatic competition discussed in this article.

Here, we restricted our analysis to deterministic equations and noises are considered indirectly in the diffusion term. As cell division and cell death are stochastic events, it is straightforward to add a noise to the cell number balance. It is well known that in classical models for interface dynamics, noises can change the wave speed and the interface profile [19, 20]. We also discussed here only tissues with reduced dimensionality, assuming translational invariance along y or rotational symmetry, respectively. Additional insight into the dynamics in three dimensions may be gained from cell-based tissue simulations [21, 22].