A population dynamics model of cell-cell adhesion incorporating population pressure and density saturation

TL;DR

This paper develops an improved continuum model for cell-cell adhesion that captures sharp invasion fronts and intermingling of different cell types, aligning well with experimental observations and applicable to tissue growth.

Contribution

It introduces a novel model incorporating localized repulsion and nonlocal attraction, enhancing the ability to simulate sharp fronts and cell intermingling in tissue growth.

Findings

Model captures sharp invasion fronts and cell intermingling.

Quantitative agreement with experimental data.

Explores effects of cell-cell repulsion on tissue dynamics.

Abstract

We discuss several continuum cell-cell adhesion models based on the underlying microscopic assumptions. We propose an improvement on these models leading to sharp fronts and intermingling invasion fronts between different cell type populations. The model is based on basic principles of localized repulsion and nonlocal attraction due to adhesion forces at the microscopic level. The new model is able to capture both qualitatively and quantitatively experiments by Katsunuma et al. (2016) [J. Cell Biol. 212(5), pp. 561--575]. We also review some of the applications of these models in other areas of tissue growth in developmental biology. We finally explore the resulting qualitative behavior due to cell-cell repulsion.

Tables

Table 1. Table 1: Reference variables used in the numerical simulations.

	Description	Dimensional Value	Dimensionless Value
$t$	Time	$1h$	0.179520
$x$	Spatial position	$1\text{µ}m$	0.01
$u(x,t),v(x,t)$	Population densities at position $x$ and time $t$	$\text{µ}{m}^{-1}$
$c_p$	Dispersivity	$1.2\times {10}^4\text{µ}{m}^3/h$
$\varphi$	A constant of proportionality related to viscosity	$1\text{µ}{m}^2/\left(h\mathrm{RU}\right)$
$R$	Sensing radius	$100\text{µ}m$
$m_1=m_2$	Crowding capacity	$2\times 0.0748\text{µ}{m}^{-1}$
$b_1=b_2$	Proliferation rate	$1/12h^{-1}$	0.464201
$k_1=k_2$	Carrying capacity	$0.0748\text{µ}{m}^{-1}$	0.5
$L$	Length of the computational domain $(-L,L)$	$1200\text{µ}m$	$12$
	Distance between the initial colonies	$500\text{µ}m$	$5$
$a_{ij}(i,j=1,2)$	Adhesion strengths
	n1N-n1N	$30\mathrm{RU}$	0.25
	n1N-n3N	$40\mathrm{RU}$	0.333333
	n3N-n3N	$8\mathrm{RU}$	0.066667
	n3NE-n3NE	$200\mathrm{RU}$	1.666667
	n3NE-n3N	$40\mathrm{RU}$	0.333333
	n3NE-n1N	$160\mathrm{RU}$	1.333333

Equations

$${\dot{x}}_i$$

$${\dot{y}}_i$$

$$\rho_{1}(x,t)\simeq\frac{m_{1}}{N}\sum_{i=1}^{N}\delta_{x_{i}(t)}\qquad\mbox{and}\qquad\rho_{2}(x,t)\simeq\frac{m_{2}}{N}\sum_{i=1}^{N}\delta_{y_{i}(t)}\qquad\mbox{ as $N\to\infty$.}$$

$$\left\{\begin{aligned} \partial_{t}\rho_{1}&=\nabla\cdot\Big{(}\rho_{1}\nabla\big{(}W_{11}\star\rho_{1}+W_{12}\star\rho_{2}+\epsilon(\rho_{1}+\rho_{2})\big{)}\Big{)},\\ \partial_{t}\rho_{2}&=\nabla\cdot\Big{(}\rho_{2}\nabla\big{(}W_{22}\star\rho_{2}+W_{21}\star\rho_{1}+\epsilon(\rho_{1}+\rho_{2})\big{)}\Big{)}.\end{aligned}\right.$$

$$\left\{\begin{aligned} &\frac{\partial u}{\partial t}=\Delta u-\nabla\cdot\left(u\mbox{\boldmath$K$}_{\!g1}(u,v)\right)+f_{1}(u,v),\\ &\frac{\partial v}{\partial t}=\Delta v-\nabla\cdot\left(v\mbox{\boldmath$K$}_{\!g2}(u,v)\right)+f_{2}(u,v).\end{aligned}\right.$$

$$\frac{\partial u}{\partial t}=-\nabla\cdot(u(\mbox{\boldmath$V$}_{\!\!\mathrm{p}}+\mbox{\boldmath$V$}_{\!\!\mathrm{a}})),$$

$$\mbox{\boldmath$V$}_{\!\!\mathrm{p}}=-\nabla p=-c_{\mathrm{p}}\nabla u.$$

$$\displaystyle\mbox{\boldmath$V$}_{\!\!\mathrm{a}}(\mbox{\boldmath$x$})=\frac{\phi}{R}\mbox{\boldmath$K$}_{\!g}(u)(\mbox{\boldmath$x$})=\frac{\phi}{R}\int_{0}^{R}\!\!\!\int_{S^{d-1}}a\,g(u(\mbox{\boldmath$x$}+r\mbox{\boldmath$\eta$}))\,\omega(r)r^{d-1}\mbox{\boldmath$\eta$}\,d\mbox{\boldmath$\eta$}dr,$$

$$g(u)=\left\{\begin{array}[]{lll}u(1-u/m)&\mathrm{if}\ \ u<m,\\ 0&\mathrm{otherwise},\end{array}\right.$$

$$\displaystyle\mbox{\boldmath$V$}_{\!\!\mathrm{a}}(\mbox{\boldmath$x$})=\frac{\phi}{R}\int_{B(0,R)}\!\!\!a\,g(u(\mbox{\boldmath$x$}+\mbox{\boldmath$y$}))\,\nabla\Omega(\mbox{\boldmath$y$})\,d\mbox{\boldmath$y$}=-\frac{\phi}{R}\int_{B(\mbox{\boldmath$x$},R)}\!\!\!a\,g(u(\mbox{\boldmath$y$}))\,\nabla\Omega(\mbox{\boldmath$x-y$})\,d\mbox{\boldmath$y$},$$

$$\frac{\partial u}{\partial t}=c_{\mathrm{p}}\nabla\cdot(u\nabla u)-\frac{\phi}{R}\nabla\cdot(u\mbox{\boldmath$K$}_{\!g}(u))+bu(1-u/k).$$

$$\begin{array}[]{cc}\displaystyle x^{*}=\frac{x}{R},\quad t^{*}=\frac{c_{\mathrm{p}}m}{R^{2}}t,\quad u^{*}=\frac{u}{m},\quad\displaystyle a^{*}=\frac{R^{d}\phi}{c_{\mathrm{p}}}a,\quad b^{*}=\frac{R^{2}b}{c_{\mathrm{p}}m},\quad k^{*}=\frac{k}{m},\end{array}$$

$$\frac{\partial u}{\partial t}=\nabla\cdot(u\nabla u)-\nabla\cdot(u\mbox{\boldmath$K$}_{\!g}(u))+f(u).$$

$$\mbox{\boldmath$K$}_{\!g}(u)(\mbox{\boldmath$x$})=a\,\int_{0}^{1}\!\!\!\int_{S^{d-1}}g(u(\mbox{\boldmath$x$}+r\mbox{\boldmath$\eta$}))\,\omega(r)r^{d-1}\mbox{\boldmath$\eta$}\,d\mbox{\boldmath$\eta$}dr,$$

$$g(u)=\left\{\begin{array}[]{lll}u(1-u)&\mathrm{if}\ \ u<1,\\ 0&\mathrm{otherwise},\end{array}\right.$$

$$f(u)=bu\left(1-u/k\right),$$

$$\left\{\begin{aligned} &\frac{\partial u}{\partial t}=\nabla\cdot\left(u\nabla(u+v)\right)-\nabla\cdot\left(u\mbox{\boldmath$K$}_{\!g1}(u,v)\right)+f_{1}(u,v),\\ &\frac{\partial v}{\partial t}=\nabla\cdot\left(v\nabla(u+v)\right)-\nabla\cdot\left(v\mbox{\boldmath$K$}_{\!g2}(u,v)\right)+f_{2}(u,v).\end{aligned}\right.$$

$$\displaystyle\mbox{\boldmath$K$}_{\!gi}(u,v)(\mbox{\boldmath$x$})=\int_{0}^{1}\!\!\!\int_{S^{d-1}}$$

$$\displaystyle+a_{i2}g_{i2}(u(\mbox{\boldmath$x$}+r\mbox{\boldmath$\eta$}),v(\mbox{\boldmath$x$}+r\mbox{\boldmath$\eta$}))\Big{]}\omega(r)r^{d-1}\mbox{\boldmath$\eta$}\,d\mbox{\boldmath$\eta$}dr,$$

$$g_{11}(u,v)$$

$$g_{22}(u,v)$$

$$f_1(u,v)=b_{1}u\left(1-(u+v)/k_1\right),$$

$$f_2(u,v)=b_{2}v\left(1-(u+v)/k_2\right),$$

$$\displaystyle\mbox{\boldmath$V$}_{\!\!\mathrm{a}}(\mbox{\boldmath$x$})$$

$$\frac{\partial u}{\partial t}=\nabla\cdot(u\nabla u)-\nabla\cdot(u(1-u)\mbox{\boldmath$K$}(u))+f(u),$$

$$\displaystyle\mbox{\boldmath$K$}(u)(\mbox{\boldmath$x$})=a\int_{0}^{1}\!\!\!\int_{S^{d-1}}u(\mbox{\boldmath$x$}+r\mbox{\boldmath$\eta$})\,\omega(r)r^{d-1}\mbox{\boldmath$\eta$}\,d\mbox{\boldmath$\eta$}dr.$$

$$\left\{\begin{aligned} &\frac{\partial u}{\partial t}=\nabla\cdot\left(u\nabla(u+v)\right)-\nabla\cdot\left(u(1-u-v)\mbox{\boldmath$K$}_{\!1}(u,v)\right)+f_{1}(u,v),\\ &\frac{\partial v}{\partial t}=\nabla\cdot\left(v\nabla(u+v)\right)-\nabla\cdot\left(v(1-u-v)\mbox{\boldmath$K$}_{\!2}(u,v)\right)+f_{2}(u,v).\end{aligned}\right.$$

$$\displaystyle\mbox{\boldmath$K$}_{\!i}(u,v)(\mbox{\boldmath$x$})=\int_{0}^{1}\!\!\!\int_{S^{d-1}}$$

$$\left\{\begin{aligned} &\frac{\partial u}{\partial t}=\nabla\cdot\left(u\nabla\chi_{1}(u,v)\right)-\nabla\cdot\left(ug_{1}(u,v)\mbox{\boldmath$K$}_{\!1}(u,v)\right)+f_{1}(u,v),\\ &\frac{\partial v}{\partial t}=\nabla\cdot\left(v\nabla\chi_{2}(u,v)\right)-\nabla\cdot\left(vg_{2}(u,v)\mbox{\boldmath$K$}_{\!2}(u,v)\right)+f_{2}(u,v).\end{aligned}\right.$$

$$\displaystyle\begin{array}[]{ll}u(x,t)=0,\\ v(x,t)=k_{2}/m_{2}=0.5\end{array}\quad\mathrm{if}\quad x\leq-L,\ t>0,$$

$$\displaystyle\begin{array}[]{ll}u(x,t)=k_{1}/m_{1}=0.5,\\ v(x,t)=0\end{array}\quad\mathrm{if}\quad x\geq L,\ t>0.$$

$$\omega(r)=\left\{\begin{array}[]{lll}a_{1}&\mbox{if}&0\leq r<R_{0},\\ a_{2}&\mbox{if}&R_{0}\leq r\leq R=1,\end{array}\right.$$

Full text

1 Introduction

More than 50 years ago, Steinberg proposed the differential adhesion hypothesis to explain the force directing cell motility during morphogenesis [35, 36, 37, 38]. Selective cell-cell adhesion between different cell types is fundamental for sorting different cell types in morphogenesis. For examples, cells dissociated from various tissues of vertebrate embryos preferentially reaggregate with cells from the same tissue when they are mixed and cultured together. This process is mediated by various cell adhesion molecules that hold cells together by homophilic or heterophilic interactions between adjacent cells. The major cell adhesion molecules in vertebrates are cadherins and nectins. Cadherins are homophilic adhesion molecules involved in various morphogenetic processes. In contrast to cadherins, nectins prefer heterophilic to homophilic adhesion.

This paper considers mathematical models capable of replicating cell sorting phenomena due to differential adhesion focusing on the different adhesive properties of molecules such as cadherins or nectins. A number of single-cell-based models have been developed and used for cell sorting phenomena, e.g., the cellular Potts models [14, 19, 22, 25, 26], lattice-free models [32, 39, 45], the vertex dynamics models [1, 21, 23] and so on. These models can show excellent agreement between numerical experiments and some cell sorting phenomena. However, these models have some drawbacks; it is necessary to set many rules and parameters for numerical calculations; multiple factors having causal relations such as adhesion and motility have to be stipulated separately; there are restrictions on the size and shape of cells, and it is difficult to consider aggregations of cells with complicated shapes such as neurons; it is very difficult to elucidate the essence of the phenomenon from analytical points of view. In order to overcome these drawbacks, cell population dynamics models are considered. Our goal in this paper is to propose a new population model to explain cell sorting due to differential adhesion.

Let us consider that we have two type of cells expressing different levels of protein surface ligands such as cadherins or nectins. Let us consider that $\{x_{i}\}$, $i=1,\dots,N$, and $\{y_{j}\}$, $j=1,\dots,M$ represent the positions of the nuclei of each of the cells forming these two groups and avoid dealing with cell membrane remodelling since we will work at the population level. Our assumption, introduced in [11], is that cells will interact with other cells either by attraction at medium distances through the formation of protusions or filopodia or by strong repulsion if the interparticle cell distance becomes very small due to volume size constraints around the nuclei. For simplicity, we assume that $N=M$ from now on. Let us also model forces exerted by cell type $j$ onto cell type $i$ as radial nonlinear springs in the direction of the centers of the nuclei, and therefore they derived from a radial potential denoted by $W^{N}_{ij}$, $i,j=1,2$. The basic agent based model for these two cell types reads as

[TABLE]

As we deal with a large number of cells, we will work at the population level, where we model cells densities rather that cell positions. With this aim, we now introduce the empirical measures associated to the above agent-based models, and describe the two populations via the mean-field scaling, the factor $1/N$, on each of these potentials with equal number of agents in order to keep a finite mass $m_{i}$, $i=1,2$, for each of the populations. We denote by $\rho_{i}$, $i=1,2$, the populations densities of each cell-type with total masses $m_{i}$, $i=1,2$, that in the limit of large number of agents $N\to\infty$, are given by the limits of the empirical measures, i.e.

[TABLE]

Here, the notation $\delta_{x_{o}}$ refers to the Dirac Delta measure at the point $x_{o}$. Let us remark that agent-based models of this type with a finite number of cells are also interesting in detailed models where differential adhesion is important, see [11] and the references therein. They include rosette formation in the early migration of the zebrafish lateral line primordium [17, 18], whose dynamics are fundamental for the correct embryonic development of the animal, and zebrafish stripe skin patterning [48].

We now put further assumptions on the scaling of these potentials reflecting the attraction for distances less than some cut-off radius $R$ together with the volume size restriction modelled by localized repulsion [9]. Each potential is scaled in $N$ in such a way that $W^{N}_{ij}\simeq\epsilon\delta_{0}+W_{ij}$ as $N\to\infty$ with $W_{ij}$ being radially symmetric, compactly supported on the ball of radius $R$, and strictly attractive inside this ball. The limit $N\to\infty$ leads to the following system for the densities $\rho_{i}$, $i=1,2$,

[TABLE]

Here, star $\star$ denotes the convolution of two functions. The rigorous derivation for one single cell type can be done following the blueprint in [30], see [6, 5] and the references therein too. This basic model for two populations shows very rich dynamical properties and complex set of stationary states and stability [8, 10, 13, 7]. However, it does not establish any upper bound on the maximal density of cells, density saturation, that is sensible for cell population models and it does not include reaction terms to take into account cell apoptosis and cell division/growth. Another well-known model was proposed by Armstrong, Painter and Sherratt [2]. They include interaction and reaction terms to model cell migration and growth/death with nonlinearities depending on the two cell densities $u,v$ while assuming linear diffusion for both, it reads as

[TABLE]

Here, the advection velocities \mbox{\boldmathK}_{\!gi}(u,v), $i=1,2$, are nonlocal terms modelling the cell adhesion with neighbouring cells with a cut-off radius $R$ but proportional not to the density of cells but to a suitable nonlinear saturation of the density to impose density saturation. We will explain this in full details below. The resulting model is able to capture concentrated densities but it does not lead to full segregation and/or sharp boundaries [4, 13] between cell types as in the case of (1) due to the linear diffusion terms and the lack of population pressure, similar drawbacks are shared by classical segregation models in population dynamics [34, 33] and references therein. Similar systems to (2) without cross-diffusion as in (1) appear in modelling cell adhesion in other mathematical biology contexts such as zebrafish patterning and tumour growth models, see [16, 15, 31, 48] for instance.

A variation of the model in (2) taking into account the population pressure coming from volume size instead of linear diffusion, as in (1), was proposed by part of the authors in [29]. They obtained sharp boundaries and segregated densities but several drawbacks due to the modelling of the migration term has been identified since then. We will elaborate on this in section 2. In this work, we propose a variation of the models in [2, 29] taking into account the population pressure coming from volume size instead of linear diffusion as in (1) together with a density saturation mechanism different from (2) for which forces are still computed linearly with respect to the local population but the saturation of the density is taken care on the mobility term. We will present the model and its improvements with respect to the previous literature in Section 2. We here validate that this model captures well the differential adhesion hypothesis [35, 36, 37, 38] explaining patterns of cell sorting and cellular movement during morphogenesis with thermodynamic principles, see a differential cell adhesion scheme in Fig. 4A below.

We will show in Section 3 that this model is able to capture qualitatively the cell sorting mechanism in [24]. Moreover, once the parameters have been carefully identified, the model is able to capture quantitatively finite interpenetration zone lengths between cell types and their invasion fronts for the experiments developed in [24]. The continuous models have some advantages to deal with problems that were difficult to be handled with conventional computable models such as single-cell-based models introduced above. Indeed, the continuous models have been applied to real problems in life sciences recently. Two applications will be reviewed briefly in Section 4. Our model can represent not only cell-cell adhesion but also cell-cell repulsion. In Section 5, we show the advantage by carrying out numerical experiments of the single-component model.

2 Comparison and derivation of cell sorting models

Let us start by introducing the Murakawa–Togashi model [29] for cell-cell adhesion. We first deal with a single population of cells for the sake of simplicity. The population density at position $x$ and time $t$ is denoted by $u(x,t)$. The model is based on the following conservation of mass:

[TABLE]

where the velocity vector is composed of velocities due to pressure \mbox{\boldmathV}_{\!\!\mathrm{p}} and adhesion \mbox{\boldmathV}_{\!\!\mathrm{a}}. Assume that the pressure $p$ is proportional to the population density, namely,

[TABLE]

Here $c_{\mathrm{p}}$ is called the dispersivity. More generally, the pressure can also be chosen as a nonlinear function of density, namely, $p=\chi(u)$. Following the argument by Armstrong, Painter and Sherratt [2], the adhesion velocity vector is given as follows:

[TABLE]

where $\phi$ is a constant of proportionality related to viscosity, $R$ is a positive constant called sensing radius, $a$ is an adhesive strength parameter, $S^{d-1}$ is the $d$-dimensional unit spherical surface. Note that in one dimension $S^{0}$ contains only two points and the integral becomes a sum. The function $g$ represents how the adhesive force depends on the local cell density, which is defined

[TABLE]

where $m$ denotes the crowding capacity of the population. The function $\omega$ describes how the force is dependent on the distance from $x$. In fact, by defining \Omega(\mbox{\boldmathx})=\Omega(|\mbox{\boldmathx}|), then \nabla\Omega(\mbox{\boldmathx})=\omega(r)\frac{\mbox{\boldmathx}}{|\mbox{\boldmathx}|}. Thus, the adhesion velocity can also be written as

[TABLE]

that is similar to the aggregation equation models studied in [28, 46, 3, 12] except for the nonlinearity $g(u)$ in the convolution. We will comment of the choice of the force $\omega(r)$ below. Adding the logistic growth term, the following model is obtained:

[TABLE]

Here, $b$ represents the growth rate and $k$ is the carrying capacity. We now rescale the system to produce a dimensionless formulation, with this goal we choose the following units

[TABLE]

and dropping the stars, the following non-dimensional model for the single population of cells is obtained.

[TABLE]

Here, the functions \mbox{\boldmathK}_{\!g} and $f$ are redefined as follows:

[TABLE]

[TABLE]

[TABLE]

reducing the model to 3 parameters $(a,b,k)$: adhesion strength, reaction strength and dimensionless carrying capacity. This derivation can be extended to a model of two types of cell populations by introducing a total population pressure leading to

[TABLE]

Here the functions \mbox{\boldmathK}_{\!gi} ($i=1,2$) are given by

[TABLE]

where $a_{ij}$ ($i,j=1,2$) denote rescaled adhesive strength parameters between the $i{\rm th}$ and $j{\rm th}$ cell populations, and the functions $g_{ij}$ are defined by

[TABLE]

[TABLE]

The growth terms $f_{1}$ and $f_{2}$ in (9) are defined by

[TABLE]

where $b_{i}$ and $k_{i}$ are the rescaled growth rate and the ratio of the carrying capacity of the $i$th population to the crowding capacity, respectively.

We now introduced a new modified model with respect to [29] to improve several of its drawbacks. Let us consider two typical situations shown in Fig. 1, where there are two cell colonies which intersect with the sensing region of a red cell, one with moderate density and the other density is almost reaching crowding capacity.

As shown in Fig. 1(I–A), according to the adhesive velocity (4), an isolated red cell moves towards the moderate density colony because the advection term \mbox{\boldmathK}_{\!gi}(u,v)(\mbox{\boldmathx}) in (7) penalizes more the high density region. However, the probability of finding other cells is higher in the right colony than in the left colony. Therefore, it is more reasonable that the force is directed towards the higher density region as shown in Fig. 1($\mathrm{II}$–A). Another unreasonable situation is depicted in Fig. 1(I–B). Here, the adhesive force due to (4) is directed towards the moderate density colony even if the red cell is surrounded by highly concentrated green cell density since both small and high densities are penalized due to the form of $g(u)$ in (8). However, it is more reasonable that the red cell does not move in such a situation as shown in Fig. 1($\mathrm{II}$–B).

To realize these ideal situations, we propose to redefine the following adhesion velocity vector instead of (4):

[TABLE]

Here, each cell counts its surrounding linearly increasing in density to determine the direction of movement. The magnitude of the total force decreases as the density locally at the cell position increases. Using the same velocity vector due to pressure population and analogous growth term as above, and applying the same rescaling as in (5), we have the following modified model for single cell population:

[TABLE]

where

[TABLE]

The choice of the radial force $\omega$ can be used as a parameter too. However, the balance between the volume exclusion term and the aggregation term leads to stationary states in the absence of the reaction term $f\equiv 0$. In Fig. 2, we observe these steady states for different choices of the potential $\Omega$.

Since their qualitative behavior is analogous for all radial attractive power-law potentials, being streched in the vertical direction for larger exponents $\lambda$, we decided to fix $\omega\equiv 1$ or equivalently $\lambda=1$ for all numerical simulations below.

All two dimensional numerical simulations in this paper are carried out in a fixed domain $[-5,5)^{2}$ with the periodic boundary condition. The standard explicit upwind finite volume method is employed and the nonlocal terms are calculated by numerical integrations (see [29] for more details).

Let us now compare two dimensional numerical solutions of (2) (Fig. 3 (B)), (6) (Fig. 3 (C)), and of (13) (Fig. 3 (D)) without growth terms. We observe how the new interaction mechanism leads to the formation of separated smooth bumps leading to stationary states of the same shape but with different masses in contrast with the previous models leading to corrugated surfaces depending on the mass of each of the bumps.

Extending this idea to two-component problem leads to the system

[TABLE]

Here the functions $f_{i}$ ($i=1,2$) are defined as in (11), and \mbox{\boldmathK}_{\!i} ($i=1,2$) are the functions \mbox{\boldmathK}_{\!gi} of (10) with $g_{i1}(u,v)=u$ and $g_{i2}(u,v)=v$, namely,

[TABLE]

Now let us look at two dimensional numerical simulations of (2), (9) and (14) for the differential adhesion hypothesis (Fig. 4). Here, we set $f_{i}\equiv 0$ ($i=1,2$). We use translucent green (resp. translucent red) to plot $u$ (resp. $v$) in order to make it easy to understand whether the two types of cells are mixed or separated. The region where $u,v<10^{-3}$ is painted in white.

Model (2) does not lead to segregation and the two types of cells are intermixed in engulfment (Fig. 4(B–$\mathrm{II}$)), partial engulfment (Fig. 4(B–$\mathrm{III}$)) and complete sorting (Fig. 4(B–$\mathrm{IV}$)) cases because of an oversimplification of the diffusion part of the model. Both models (9) (Fig. 4(C)) and (14) (Fig. 4(D)) can replicate different cell sorting patterns explained by the differential adhesion hypothesis and show sharp boundaries between the green and red cells. As seen in one-component numerical simulations, Model (14) leads to smooth bumps, whereas Model (9) leads to corrugated surfaces. For the mixing pattern Fig. 4(I), Model (2) and (9) show that the green and red cells are intermixed about fifty-fifty. In contrast, in Model (9), these cells form marble pattern colonies. Since the green and red cells are expressing the same type of nectin molecule in these cases and these cells are scattered apart in the initial state, the state of being mixed randomly like Fig. 4($\mathrm{I}$–D) is more realistic. Incidentally, the numerical solutions $u$(green)+$v$(red) in Fig. 4 ($\mathrm{I}$–B), ($\mathrm{I}$–C) and ($\mathrm{I}$–D) coincide with those in Fig. 3 (B), (C) and (D), respectively. Thus, the colonies in Fig. 4($\mathrm{I}$–D) are smooth.

In this paper, we suggested specific choices of the population pressure in (3) and the density saturation in (12). Of course, one can consider more general form of the model as follows:

[TABLE]

Here, $\chi_{i}$ and $g_{i}$ ($i=1,2$) are given functions.

3 Numerical simulations for the Togashi et al.’s experiment

Nectins are immunoglobulin-like cell-adhesion molecules that compose a family of four members (nectin-1, -2, -3, and -4) [40]. Nectins prefer heterophilic interactions to homophilic ones, in contrast with the homophilic nature of cadherin interactions. For example, when cells expressing nectin-1 and -3 are mixed, these cells are arranged in a mosaic pattern due to their heterophilic interactions [42, 41].

Top figures of each Fig. 5(A)–(D) show time-lapse microscopy of a mosaic forming assay of HEK293 cells expressing different nectins and cadherins. Depending on the combinations of cell-adhesion molecules, different patterns are observed. Note that HEK293 cells naturally express N-cadherin, but not E-cadherin. Before the time-lapse experiments, these cells were separately cultured to allow the formation of independent colonies. When their colony edges came into contact with one another, their boundaries were examined. The shape variations of the colony edges at the start point are dependent on the experimental manipulation. Time $t$ indicates elapsed time in hours. In (A), HEK293 cells expressing nectin-1, N-cadherin and EGFP (green) and those expressing nectin-1, N-cadherin and mCherry (red) were cultured. In early stage of culture, cells grow proliferously and expand their habitats. Each red cell moves from the left to the right and each green cell moves from the right to the left. And then, they fill in the gap around $t=30$h. They do not overlap at later times in this case. We can clearly see the boundary between green and red colonies. There is no driving force for intermingling with one another because the green and red cells are the same type. After the density reaches the carrying capacity, they stop growing and moving. In (B), HEK293 cells expressing nectin-1, N-cadherin and EGFP (green) and those expressing nectin-3, N-cadherin and mCherry (red) were cultured. In this case, they are intermingling with one another after contact, because each nectin-1+ green cell prefers to adhere to nectin-3+ red cells and vice versa.

In (C), HEK293 cells expressing nectin-3, N-cadherin and E-cadherin-EGFP (green) and those expressing nectin-3, N-cadherin and mCherry (red) were cultured. These cells do not invade the counter colony each other. However, the red cells push back the counter colony of green cells. In (D), HEK293 cells expressing nectin-3, N-cadherin and E-cadherin-EGFP (green) and those expressing nectin-1, N-cadherin and mCherry (red) were cultured. The red cells invade the counter colony of green cells, but the green cells do not. The habitat (front position) of green cells does not change at later times.

We perform numerical simulations for these experiments in one space dimension using (14). The following actual measured values are used for the dimensional equations: $R=100\ \mbox{\textmu}\mathrm{m}$, $b_{1}=b_{2}=1/12\,\mathrm{h}^{-1}$, $k_{1}=k_{2}=0.0748\,\mbox{\textmu}\mathrm{m}^{-1}$. The sensing radius $R$ is determined by the length from the cell body to the tip of leading edge. The distance between the initial colonies is $500\,\mbox{\textmu}\mathrm{m}$. The dispersivity $c_{\mathrm{p}}$ controls the spreading speed of the colony edges. We set $c_{\mathrm{p}}=1.2\times 10^{4}\,\mbox{\textmu}\mathrm{m}^{3}/\mathrm{h}$ such that the two colonies edges come into contact at 30 hours of incubation in Case (A). We could not measure the values $\phi$ and $m_{i}$, and simply chose $\phi=1\,\mbox{\textmu}\mathrm{m}^{2}/(\mathrm{h}\,\mathrm{RU})$. It is thought that $m_{i}$ is greater than $k_{i}$ because the cells are moving due to cell adhesion even after the cells fill up the region. We chose $m_{i}=2k_{i}$. We consider the domain of length $2400\,\mbox{\textmu}\mathrm{m}$, that is, the computational domain is $(-L,L)$ with $L=12$. In order to compare numerical solutions with the experimental results, the visualized domain is set to be $[-4,4]$. The following Dirichlet boundary conditions are imposed:

[TABLE]

We fix all parameters of the above for every simulations (A)–(D). The adhesive strength parameters are the only difference.

We determined the adhesive strengths from the data given by Harrison et al. [20] and Katsunuma et al. [24]. Hereafter, we use ellipsis notation for the adhesive strength, e.g., n3NE-n1N implies the adhesive strength between nectin-3, N-cadherin and E-cadherin expressing cells and nectin-1 and N-cadherin expressing cells. We chose $\mbox{n1N-n1N}=30\,\mathrm{RU}$(Response Unit) from Fig. 1a in [20], $\mbox{n3N-n3N}=8\,\mathrm{RU}$ from Fig. 1c in [20], and $\mbox{n1N-n3N}=40\,\mathrm{RU}$ from Fig. 1c in [20]111We chose $\mbox{n1N-n3N}=230\,\mathrm{RU}$ in [29] using the datum from Fig. 1a in [20]. But the datum ($=40\,\mathrm{RU}$) from Fig. 1c in [20] is consistent with the experiment (B).. Since n3NE-n3NE is about 8 times of n2N-n2N (Fig. 8E in [24]), and $\mbox{n2N-n2N}=25\,\mathrm{RU}$ from Fig. 1b in [20], we chose $\mbox{n3NE-n3NE}=8\times 25\,\mathrm{RU}=200\,\mathrm{RU}$. There is no data available for n3NE-n3N nor n3NE-n1N. Assuming that n3NE-n3N (resp. n3NE-n1N) is similar to n2NE-n2N (resp. n3NE-n2N), and using the data in Fig. 1b in [20] and in Fig. 8E in [24], we chose $\mbox{n3NE-n3N}=40\,\mathrm{RU}$ and $\mbox{n3NE-n1N}=160\,\mathrm{RU}$. All parameters are summarized in Table 1.

The numerical results of (14) are shown in bottom figures of each Fig. 5(A)–(D). We use green line to plot $u$ and red line to plot $v$, and total population density $u+v$ is drawn with a blue line behind green and red lines. Appearance of blue line implies that the two types of cells are mixed. Under the green (resp. red) line is painted in translucent green (resp. pink) so that we can easily see whether the two populations are mixed or not. The step functions are used for the initial data in Cases (A) and (B). But, in the experiments in Cases (C) and (D), these cells are already moving at $t=0$, and the distributions at $t=0$ are similar to those in Cases (A) and (B) at $t=24$ (Watch Videos 3-6 in [24] or Video S1). Furthermore, the spreading speeds of the green colonies are slower than those of the red colonies. Therefore, the initial data of Cases (A) and (B) shifted a little to the left are set at $t=-24$ as the initial data in Cases (C) and (D), and we plot the numerical results from $t=0$. The numerical results demonstrate excellent agreement between the model (14) and the experiments, and illustrate that the model (14) is able to replicate each pattern not only qualitatively but also quantitatively. The interested readers can find Video S1 of the experiments versus the simulations using (14) in the supplementary material.

Katsunuma et al. [24] performed lattice-based numerical simulations in order to understand the asymmetric pattern formation in Fig. 5 (D). They noted “our time-lapse observations revealed a significant difference in the migratory behavior of 293 cells, which was dependent of E-cadherin expression. … The 293 cells expressing E-cadherin migrated more slowly toward the counter colony compared with those that did not express E-cadherin”. Hence, they added cell mobility to their computational model independently of cell adhesion, and concluded that differential mobility, in addition to adhesiveness, was responsible for the asymmetric mingling patterning even though there is a causal relation between cell adhesion and cell mobility as they mentioned. In contrast, Model (14) is capable of replicating the differential cell mobility as a consequence of the differential adhesion. The density of cells expressing E-cadherin near the propagating front is larger than that of cells which did not express E-cadherin. Model (14) can also replicate such phenomenon naturally.

Finally, we turn our attention to the comparison of our modelling approach with respect to individual based modelling for cell adhesion. When single-cell-based models such as the cellular Potts models [14, 19, 22, 25, 26], lattice-free models [32, 39, 45] and the vertex dynamics models [1, 21, 23] are used, rules on adhesion, mobility, density, etc. have to be defined independently even if they have a causal relationship. And then, the number of parameters increases and the simulations become complicated and difficult to track and parameterize. Our model has the advantage of connecting the individual based modelling to macroscopic populations models as explained in the introduction with a minimal set of parameters amenable for easier tuning.

4 Applications

There are tight limitations in size and shape of the cells in single-cell-based models. However, the terminals of neurons show very irregular shapes that change dynamically. When investigating cell aggregations of such dynamic cells, it would be useful to use a continuous model. Here, we give brief summaries of two applications.

4.1 Role of Reelin during layer formation in the mannmalian cerebral neocortex

The mannmalian cerebral neocortex has a highly organized layered structure of neurons. During cortical development, a glycoprotein, Reelin, plays a crucial role in the neuronal migration and neocortical lamination because Reelin-deficient mice show disrupted neuronal organization. However, its precise role in neuronal layer organization remained unclear. It was thought that Reelin promotes neuronal adhesion and induces neuronal cell aggregation. To uncover how Reelin controls the intercellular adhesion among cortical cells, Matsunaga et al. [27] performed Reelin stimulation experiments. Two types of cells, radial glial cells (RGCs) and neurons (Neus), played key roles there. If Reelin promotes N-cadherin-dependent neuronal adhesion directly, the differential adhesivenesses should be in the order of Neu-Neu$>$Neu-RGC$>$RGC-RGC in the Reelin stimulation experiments. In this situation, according to the differential adhesion hypothesis, the neurons form clusters surrounded by the RGCs. However, the actual results were opposite, that is, the RGC clusters were engulfed by the neurons. A mathematical model based on (9) was used to understand the reason for this unexpected clustering pattern. The model predicted that Reelin causes an increase of cell-cell adhesion among neurons transiently but not persistently. Many expriments revealed that the prediction is correct. Finally, it was concluded that transient but not persistent increase in cell-cell adhesion might be necessary for the highly organized layered structure of neurons in the mammalian neocortex [27].

4.2 Role of differential adhesion during columnar unit formation in the Drosophila brain

In the developing visual center of the Drosophila brain, multiple neurons gather to form a columnar structure, a basic morphological and functional unit of the brain. Three types of core columnar neurons, R7, R8 and Mi1, play key roles at the initial step of column formation along a two dimensional layer of the brain during larval stage. A series of biological studies demonstrated that the differential adhesiveness of these neurons in the order of R7$>$R8$>$Mi1 causes their concentric arrangement. As a result, the terminal of R7 occupies the dot-like central region, the R8 terminal enwraps the R7 terminal forming a donut-like region and the Mi1 terminal occupies a grid-like region outside the R8 terminal. Since the neurites of columnar neurons show very irregular shapes that change dynamically, we formulated a three-component model based on (14) by considering the density of neurites and differential adhesion between them. Our model demonstrated that the differential adhesion among R8, R7 and Mi1 is sufficient to reproduce the wild type and mutant patterns of the columns [47].

5 Further numerical experiments of the single-component model

5.1 Patterning under different initial data

Cells form cell masses by adhesion, but the shape and size of the cell mass greatly differs according to the initial datum. Here, we show several different patterns depending on the difference of the initial data. Numerical experiments of Equation (13) are performed with fixed $a=4$, $\omega\equiv 1$ and $f\equiv 0$. The initial datum is set as a constant $V$ perturbed with 1% random noise. Figure 6 shows that bell-shaped, striped and perforated patterns occur under variation of the initial cell density. The colony size and shape are different for each colony, and depend on the initial data.

5.2 Repulsive interaction and volume-dependent pattern formation

Cell sorting is occasionally caused not by cell-cell adhesion but by cell-cell repulsion, see, e.g., [43, 44, 45] and references therein. The cell-cell repulsion is represented by our model just setting the adhesion parameter $a$ negative or setting $\omega$ negative. Here, recall that the function $\omega$ describes how the force is dependent on the distance from the cell. We carry out numerical simulations for Equation (13) with fixed $a=1$ and $f\equiv 0$. The function $\omega$ is set to be

[TABLE]

where $0<R_{0}<1$ and $a_{1}$ and $a_{2}$ can be positive or negative. For example, when $a_{1}$ is positive and $a_{2}$ is negative, it indicates short-range attractive and middle-range repulsive potential. We use the same initial datum as above. Figure 7 shows numerical results with fixed $R_{0}=0.85$, $a_{1}=4$ and $a_{2}=-32$, and different choices of $V$ that controls total mass. When the total mass is small, a spotted pattern arranged in order is observed (Fig. 7A). The distance among the spots seems to be determined by $R_{0}$ and $R$. Here, we recall that the computational domain is a square with side length 10. When the total mass exceeds a certain threshold, a striped pattern appears (Fig. 7B). When the total mass is increased further, a honeycomb pattern appears (Fig. 7C). When the total mass is large enough, the density becomes homogeneous (Fig. 7D). These patterns are robustly generated depending on the total mass, and distinct from the type of patterns in Fig. 6.

Thus, various patterns appear even with a single component. As mentioned in Section 4.2, our model was applied to understanding the formation of columnar units in the Drosophila brain. Multiple columnar units are neatly arranged and form a hexagonal lattice just like Fig. 7A. Our model might be able to explain the formation of such structures. Besides this, similar patterns to Fig. 7 appear in vegetation patterns [49] and patterns in animal skin [50]. Our model has some advantages to elucidate such pattern formations.

6 Concluding remarks

In this paper, we discussed and improved continuum models for cell-cell adhesion. Armstrong, Painter and Sherratt [2] proposed a celebrated model consisting of linear diffusions and nonlocal advection terms. However, it gives biologically unrealistic numerical solutions. In particular, it can not replicate mutually-immiscible phenomena. The underlying cause of the problem was that the model is based on random movement of each individual cell at the microscopic level. Murakawa and Togashi [29] have changed the basic principles of cell movement from “cells move randomly” to “cells move from high to low pressure regions”, and have proposed a modified continuous model for cell-cell adhesion based on the total population pressure. Their model was able to replicate the mutually-immiscible phenomena and the different types of sorting patterns. However, this model presents a drawback since it leads to corrugated surfaces of cell colonies which are not biologically reasonable. In this paper, we proposed a modified model by rethinking of the nonlocal adhesion terms and imposing the saturation response to adhesive forces in a different manner. Numerical solutions of the modified model show smooth surfaces of cell colonies while being able to capture qualitatively the cell sorting mechanism. By a careful parameter choice, we are able to obtain an excellent agreement with the phenomena observed in experimental data by Katsunuma et al. [24]. We also show how the models depends on the potential shape and we have explored not only cell-cell adhesion but also cell-cell repulsion.

Acknowledgments

JAC was partially supported by the EPSRC grant number EP/P031587/1. HM was partially supported by JSPS KAKENHI Grant Numbers 26287025, 15H03635 and 17K05368, and JST CREST Grant No. JPMJCR14D3. HT was supported by JSPS KAKENHI Grant Numbers 18K06219, 18H04764, and a grant from the Takeda Science Foundation.

Supplementary materials