1. Introduction

The mathematical formalisation of cell motion is a fascinating topic which has attracted the attention of physicists and mathematicians over the past fifty years. In more recent times, the development of sophisticated technologies capable of capturing high-resolution videos of moving cells has renewed the interest of the physical and mathematical communities. This has promoted the formulation of several models, which rely on different mathematical approaches, to reproduce qualitative behaviours emerging from cell motion. We refer the interested reader to[1, 2, 5, 7, 8, 9, 10, 11, 14, 17, 18, 22, 31, 36, 37] and references therein.

Accumulating evidence indicates that the interaction between epithelial and mesenchymal cells plays a pivotal role in cancer development and metastasis formation [20, 39, 40]. Here we propose a model for the co-culture dynamics of epithelial-like and mesenchymal-like cells. In our model, mesenchymal-like and epithelial-like cells in motion are grouped into two distinct populations. Following a kinetics approach, the microscopic state of each cell is identified by its position and velocity (i.e., a point in the phase space), and the two populations are characterised through the related phase space densities. Macroscopic quantities (e.g., local cell densities) are then expressed in terms of microscopic averages [3, 4]. We make use of the mesoscopic formalism developed in[8, 22] to include microscopic aspects of cell motion which cannot be captured by macroscopic models. Moreover, we adapt the strategies presented in[15] to model proliferation, competition for nutrients and adhesive interactions between epithelial-like cells.

A calibration of the model based on experimental datasets is beyond our present scope. In this work, our aim, rather, is to present a sample of numerical results which display the emergence of different spatial patterns, depending on parameters and initial data. The remainder of the paper is organised as follows. The main features of the biological system and the mathematical model are introduced in Section 2. In Section 3, we discuss the main numerical results. Finally, some research perspectives are summarised in Section 4.

---

2. The model

In this section, we introduce the key features of the biological system and the mathematical model. In Subsection 2.1, we state the biological problem and the assumptions we need in view of the mathematical formalisation. In Subsection 2.2, we describe the modelling strategies we designed to translate the biological phenomena into mathematical terms, and we present the model.

---

2.1. Brief phenomenological overview and main assumptions

The term Epithelial-Mesenchymal Transition (EMT) refers to the temporary and reversible switching between epithelial and mesenchymal phenotypes. During EMT, non-motile epithelial cells, collectively embedded via cell-cell junctions (i.e., homotypic adhesion), convert into motile mesenchymal cells [20, 39]. This strongly reduces cell adhesion and enhances cell motility.

EMT is commonly observed in various non-pathological conditions, namely during embryonic development and tissue repair in the adult organism. However, this phenotypic reprogramming has been linked to cancer progression [39,40]. For instance, EMT may favour the seeding of secondary tumours at distant sites and the creation of metastases [20]. In particular, EMT has been implicated in the formation of the fibrous capsules observed in hepatic tumours, which seem to be mainly composed of cells that express a mesenchymal-like phenotype [25].

We consider a monolayer of epithelial-like and mesenchymal-like cells in co-culture on a regular plastic surface (i.e., a petri dish). Cellular movement is seen as the superposition of persistent spontaneous motion and chemotactic response. The former is due to the tendency of cells to randomly orient themselves, whilst the latter is guided by chemotactic cytokines. In this framework, we focus on the following biological phenomena:

(ⅰ) secretion from cells of chemotactic cytokines;

(ⅱ) diffusion and consumption of chemotactic cytokines;

(ⅲ) random motion of cells and chemotaxis;

(ⅳ) adhesive interactions between epithelial cells;

(ⅴ) cell proliferation and competition for nutrients.

To reduce biological complexity to its essence, we make the prima facie assumption that the diffusion of chemotactic cytokines is isotropic. Moreover, we model cell motion in two dimensions only, since we focus on a monolayer co-culture (i.e., we let cells grow side by side and not one on top of the other), and we assume that the status of motion of a cell is left unaltered by interactions which do not lead to homotypic adhesion. Finally, we stress that this paper does not deal with several evolutionary aspects. For instance, we consider a sample where both epithelial-like and mesenchymal-like cells are present from the beginning, and we do not let EMT occur. Moreover, we do not include the effects of cellular heterogeneity and therapeutic actions, which we explored previously in[12,15, 16, 28, 29].

---

2.2. Mathematical formalisation

We divide epithelial-like and mesenchymal-like cells in motion into two populations labeled, respectively, by the index $i=1$ and $i=2$. We identify the microscopic state of each moving cell by its instantaneous position and velocity, which are described by the continuous variables $\mathbf{x} = (x,y) \in X \subset \mathbb{R}^2$ and $\mathbf{v} = (v_x,v_y) \in V \subset \mathbb{R}^2$. We let the set $V$ be compact and spherically symmetric [8]. Furthermore, we group epithelial-like cells at rest due to homotypic adhesion and chemotactic cytokines into two additional populations, labeled by the index $i=3$ and $i=4$, and we model their microscopic states by means of the $\mathbf{x}$ variable only.

At any instant of time $t$, these populations are characterised, respectively, by the phase space densities

$ f_1 = f_1 (t,\mathbf{x},\mathbf{v}): \mathbb{R}_+ \times X \times V \to \mathbb{R}_+, f_2 = f_2 (t,\mathbf{x},\mathbf{v}): \mathbb{R}_+ \times X \times V \to \mathbb{R}_+ \nonumber $

and the local densities

$ n_3 = n_3 (t,\mathbf{x}): \mathbb{R}_+ \times X \to \mathbb{R}_+, n_4 = n_4 (t,\mathbf{x}): \mathbb{R}_+ \times X \to \mathbb{R}_+. \nonumber $

The local densities of cells in motion and the local cell density can be computed as

$ n_{1,2}(t,\mathbf{x}) = \int_{V} f_{1,2}(t,\mathbf{x},\mathbf{v}) d\mathbf{v}, \varrho(t,\mathbf{x}) = \sum\limits_{i=1}^3 n_i(t,\mathbf{x}), $

(1)

while the total cell density and the average local cell density are

$ N(t) = \int_{X} \varrho(t,\mathbf{x}) d\mathbf{x}, \bar{\varrho}(t) = \frac{N(t)}{|X|}, $

(2)

where $|X|$ denotes the measure of the set $X$.

The following notations and assumptions are used to model the biological phenomena of interest (vid. Table 1 for a summary of the model parameters):

$\bullet$ Secretion from cells of chemotactic cytokines. Cells in population $i=1,2,3$ secrete cytokines responsible for chemotaxis at an average rate $\nu \in \mathbb{R}_+$.

$\bullet$ Diffusion and consumption of chemotactic cytokines. Chemotactic cytokines diffuse across the culture system with unitary diffusion constant. Moreover, the concentration of cytokines decreases over time due to the consumption by cells, which occurs at a rate described by the functional $\mathcal{A}(\varrho(t,\mathbf{x});\alpha) \geq 0$. The parameter $\alpha \in \mathbb{R}_+$ stands for the average consumption rate and, as a first attempt to take into account the nonlinear nature of cell-cytokine interactions [6], we define the functional $\mathcal{A}$ as

$ \mathcal{A}(\varrho;\alpha) := \alpha \; \varrho(t,\mathbf{x})^2. $

(3)

$\bullet$ Random motion of cells and chemotaxis. Using a velocity jump formalism [8, 22], we assume that a cell in the state $(\mathbf{x},\mathbf{v}_*)$ of the population $i=1,2$ can jump into the state $(\mathbf{x},\mathbf{v})$ of the same population at a rate described by the functional

$ \mathcal{B}(\mathbf{v}, \nabla_{\mathbf{x}} n_4; \beta, |V|) := \frac{1}{|V|} \; + \; b(|\mathbf{v}|, |\nabla_{\mathbf{x}} n_4|; \beta, |V|) \left(\mathbf{v} \cdot \nabla_{\mathbf{x}} n_4(t,\mathbf{x}) \right) $

(4)

with

$ b(|\mathbf{v}|, |\nabla_{\mathbf{x}} n_4|; \beta, |V|) : = \frac{|V|^{-1}}{\beta + |\mathbf{v}| |\nabla_{\mathbf{x}} n_4(t,\mathbf{x})|}. $

(5)

The functional $\mathcal{B}$ satisfies

$ \int_V \mathcal{B}(\mathbf{v}, \cdot; \beta, |V|) d\mathbf{v} = 1, \forall \beta, |V| \in \mathbb{R}_+. $

(6)

In definition (4), $|V|$ denotes the measure of the set $V$, and the parameter $\beta \in \mathbb{R}_+$ provides an average measure of the relative contribution of random motion versus chemotactic reorientation on cell movement. The first term of definition (4) translates into mathematical terms the idea that random motion does not favour any precise velocity (i.e., when a moving cell reorients itself and jump from the state $(\mathbf{x},\mathbf{v}_*)$ to the state $(\mathbf{x},\mathbf{v})$, all values of $\mathbf{v}$ are equiprobable). In the second term, the scalar product models the tendency of cells to follow the cytokine gradient independently from their current velocity (i.e., the chemotactic response leads a moving cell in the state $(\mathbf{x},\mathbf{v}_*)$ to jump to the state $(\mathbf{x},\mathbf{v})$, with $\mathbf{v}$ pointing in the direction of $\nabla_{\mathbf{x}} n_4$ independently from the direction of the vector $\mathbf{v}_*$). Definition (5) preserves the non-negativity of the functional $\mathcal{B}$, and ensures that smaller values of the parameter $\beta$ will correspond to a stronger influence of chemotactic reorientation on cell motion.

$\bullet$ Adhesive interactions between epithelial cells. Encounters between a cell in the state $(\mathbf{x},\mathbf{v}^*)$ of population $j=1$ (or the state $\mathbf{x}$ of population $j=3$) and a cell in the state $(\mathbf{x},\mathbf{v}_*)$ of population $i=1$ can lead the latter to reach either the state $(\mathbf{x},\mathbf{v})$ of the same population (i.e., the adhesive interaction fails) or the state $\mathbf{x}$ of population $h=3$ (i.e., the adhesive interaction leads to homotypic adhesion). For the sake of simplicity, we assume that encounters between epithelial-like cells occur at unitary rate, and we define the rate of adhesive interactions as

$ \mathcal{G}^h_{ij}(\mathbf{v}; \gamma, |V|) := $$\begin{cases} \displaystyle{\frac{\gamma}{|V|}},&\mbox{if }\; i = 1, \; j=1,3 \; \mbox{and } h = 3 \\ \\ \displaystyle{\frac{1 - \gamma}{|V|}},&\mbox{if }\; i = 1, \; j=1,3 \; \mbox{and } h = 1 \\ \\ 0,& \mbox{otherwise}, \end{cases}$$ \forall \; \mathbf{v} \in V, $

(7)

$ \sum\limits_{h=1}^{3} \int_V \mathcal{G}^h_{ij}(\mathbf{v}; \gamma, |V|) d\mathbf{v} = 1, \mbox{for} \; \gamma \in [0,1], \; |V| \in \mathbb{R}_+, \; i=1 \mbox{ and } j=1,3. $

(8)

In definition (7), the parameter $\gamma \in [0,1]$ stands for the average rate of homotypic adhesion, while the factor $1/{|V|}$ accounts for the fact that the interactions under consideration are independent from the velocities of the interacting cells.

$\bullet$ Cell proliferation and competition for nutrients. The parameter $\kappa \in \mathbb{R}_+$ stands for the average cell proliferation rate. Furthermore, since the proliferation of cells is limited by the competition for nutrients, we introduce a death term $\mathcal{M}(N;\mu)$. The parameter $\mu \in \mathbb{R}_+$ stands for the average rate of death. Among all possible definitions of the functional $\mathcal{M}$, we make use of the following one

$ \mathcal{M}(N;\mu) := \mu \; N(t), $

(9)

which relies on the natural assumption that a higher total density of cells corresponds to a lower concentration of available nutrients, and thus to a higher chance of cell death.

The evolution of the functions $f_1(t,\mathbf{x},\mathbf{v})$, $f_2(t,\mathbf{x},\mathbf{v})$, $n_3(t,\mathbf{x})$ and $n_4(t,\mathbf{x})$ is governed by the equations given hereafter, which describe the net inlet of moving epithelial-like and mesenchymal-like cells through the volume element $d\mathbf{x} \; d\mathbf{v}$ centered at $(\mathbf{x},\mathbf{v})$, as well as the net flux of epithelial-like cells at rest and chemotactic cytokines through the volume element $d\mathbf{x}$ centered at $\mathbf{x}$:

$ $$\begin{array}{l} \partial_t f_1(t,\mathbf{x},\mathbf{v}) + \mathbf{v} \cdot \nabla_\mathbf{x} f_1(t,\mathbf{x},\mathbf{v})\nonumber\\= \underbrace{\int_V \mathcal{B}(\mathbf{v}, \nabla_{\mathbf{x}} n_4; \beta, |V|) f_1(t,\mathbf{x},\mathbf{v}_*) d\mathbf{v}_* - f_1(t,\mathbf{x},\mathbf{v})}_{\mbox{ inflow & outflow due to random motion and chemotactic reorientation}} \nonumber\\ +\underbrace{\int_V \int_V \mathcal{G}^1_{11}(\mathbf{v}; \gamma, |V|) f_1(t,\mathbf{x},\mathbf{v}_*) f_1(t,\mathbf{x},\mathbf{v}^*) d\mathbf{v}_* d\mathbf{v}^*}_{\mbox{ inflow due to changes of velocity caused by adhesive interactions}} \nonumber\\ +\underbrace{n_3(t,\mathbf{x}) \int_V \mathcal{G}^1_{13}(\mathbf{v}; \gamma, |V|) f_1(t,\mathbf{x},\mathbf{v}_*) d\mathbf{v}_*}_{\mbox{ inflow due to changes of velocity caused by adhesive interactions}} \nonumber\\ -\underbrace{f_1(t,\mathbf{x},\mathbf{v}) \int_V \int_V \Big(\mathcal{G}^1_{11}(\mathbf{v_*}; \gamma, |V|) + \mathcal{G}^3_{11}(\mathbf{v_*}; \gamma, |V|) \Big) f_1(t,\mathbf{x},\mathbf{v}^*) d\mathbf{v}_* d\mathbf{v}^*}_{\mbox{ outflow due to homotypic adhesion}} \nonumber\\ -\underbrace{f_1(t,\mathbf{x},\mathbf{v}) \; n_3(t,\mathbf{x}) \int_V \Big(\mathcal{G}^1_{13}(\mathbf{v_*}; \gamma, |V|) + \mathcal{G}^3_{13}(\mathbf{v_*}; \gamma, |V|) \Big) d\mathbf{v}_*}_{\mbox{ outflow due to homotypic adhesion}} \nonumber\\ +\underbrace{\kappa \; f_1(t,\mathbf{x},\mathbf{v}) - \mathcal{M}(N;\mu) \; f_1(t,\mathbf{x},\mathbf{v})}_{\mbox{ proliferation and competition}}, \end{array}$$ $

(10)

$ $$\begin{array}{l} \partial_t f_2(t,\mathbf{x},\mathbf{v}) + \mathbf{v} \cdot \nabla_\mathbf{x} f_2(t,\mathbf{x},\mathbf{v}) \!\!\!\!\ &=& \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \underbrace{\int_V \mathcal{B}(\mathbf{v}, \nabla_{\mathbf{x}} n_4; \beta, |V|) f_2(t,\mathbf{x},\mathbf{v}_*) d\mathbf{v}_* - f_2(t,\mathbf{x},\mathbf{v})}_{\mbox{ inflow & outflow due to random motion and chemotactic reorientation}} \nonumber\\ & +&\underbrace{\kappa \; f_2(t,\mathbf{x},\mathbf{v}) - \mathcal{M}(N;\mu) \; f_2(t,\mathbf{x},\mathbf{v})}_{\mbox{ proliferation and competition}}, \end{array}$$ $

(11)

$ $$\begin{array}{l} \partial_t n_3(t,\mathbf{x})&=&\underbrace{\int_V \int_V \int_V \mathcal{G}^3_{11}(\mathbf{v}; \gamma, |V|) f_1(t,\mathbf{x},\mathbf{v}_*) f_1(t,\mathbf{x},\mathbf{v}^*) d\mathbf{v}_* d\mathbf{v}^* d\mathbf{v}}_{\mbox{ inflow due to homotypic adhesion}}\\& +& \underbrace{n_3(t,\mathbf{x}) \int_V \int_V \mathcal{G}^3_{13}(\mathbf{v}; \gamma, |V|) f_1(t,\mathbf{x},\mathbf{v}_*) d\mathbf{v}_* d\mathbf{v}}_{\mbox{ inflow due to homotypic adhesion}} \nonumber\\&+&\underbrace{\kappa \; n_3(t,\mathbf{x}) - \mathcal{M}(N;\mu) \; n_3(t,\mathbf{x})}_{\mbox{ proliferation and competition}}, \end{array}$$ $

(12)

$ \partial_t n_4(t,\mathbf{x}) = \underbrace{\nu \; \varrho(t,\mathbf{x}) + \Delta_\mathbf{x} n_4(t,\mathbf{x})}_{\mbox{ secretion and diffusion}} - \underbrace{\mathcal{A}(\varrho;\alpha) \; n_4(t,\mathbf{x}).}_{\mbox{ consumption}} $

(13)

Plugging definitions (3), (4), (5), (7) and (9) into equations (10)-(13), we obtain

$ $$\begin{array}{l} \partial_t f_1 + \mathbf{v} \cdot \nabla_\mathbf{x} f_1 &=& \underbrace{\frac{1}{|V|} \left(1+ \frac{\mathbf{v} \cdot \nabla_\mathbf{x} n_4}{\beta + |\mathbf{v}| |\nabla_\mathbf{x} n_4|} \right) n_1 - f_1}_{\mbox{ random motion and chemotactic reorientation}} \nonumber\\ &+& \underbrace{ \frac{(1-\gamma)}{|V|} n_1\left (n_1 + n_3 \right) - (n_1+n_3) f_1}_{\mbox{ homotypic adhesion}} \nonumber\\ &+& \underbrace{\left(\kappa - \mu \; N \right) f_1,}_{\mbox{ proliferation and competition}} \end{array}$$ $

(14)

$ $$\begin{array}{l} \partial_t f_2 + \mathbf{v} \cdot \nabla_\mathbf{x} f_2 &=& \underbrace{\frac{1}{|V|}\left(1 + \frac{\mathbf{v} \cdot \nabla_\mathbf{x} n_4 }{\beta + |\mathbf{v}| |\nabla_\mathbf{x} n_4|} \right) n_2 - f_2}_{\mbox{ random motion and chemotactic reorientation}} \nonumber\\ &+& \underbrace{\left(\kappa - \mu \; N \right) f_2,}_{\mbox{ proliferation and competition}} \end{array}$$ $

(15)

$ $$\begin{array}{l} \partial_t n_3&= & \underbrace{ \gamma \left(n_1 + n_3 \right)n_1}_{\mbox{ homotypic adhesion}} + \underbrace{\left( \kappa - \mu \; N \right) n_3,}_{\mbox{ proliferation and competition}} \end{array}$$ $

(16)

$ $$\begin{array}{l} \partial_t n_4 &=& \underbrace{\nu \; \varrho + \Delta_\mathbf{x} n_4}_{\mbox{ secretion and diffusion}} - \underbrace{ \alpha \; n_4 \varrho^2.}_{\mbox{ consumption}} \end{array}$$ $

(17)

---

3. Main results

In this section, we discuss a sample of numerical results which display the emergence of different spatial patterns, depending on parameters and initial data. In Subsection 3.1, we describe the simulation setup and the method we use for calculating numerical solutions. In Subsection 3.2, we focus on a sample composed of mesenchymal-like cells only, where proliferation and competition phenomena do not take place. The results we present highlight how different initial cell densities can cause the emergence of different spatial patterns, such as spots, stripes and hole structures. In Subsection 3.3, we show that, when there is little competition for nutrients, epithelial-like and mesenchymal-like cells can segregate and create honeycomb structures. Finally, the simulations discussed in Subsection 3.4 provide a possible explanation for how the interplay between epithelial-cell adhesion and mesenchymal-cell spreading paves the way for the formation of ring-like structures.

---

3.1. Numerical method and simulation setup

For simplicity, we follow [33] and approximate cellular velocities in polar coordinates. We also assume that all moving cells are characterised by the same modulus of the velocity, which is normalised to unity. This can be justified by noting that the main differences between the adhesive behaviours of epithelial-like and mesenchymal-like cells are already captured by the modelling strategies described in the previous section. As a result, we set

$ v_x = \cos\theta, v_y = \sin\theta, \theta \in [0,2\pi), $

(18)

so that the phase space distributions can be computed by tracking velocity angles and space only.

To perform numerical simulations, we set $X = [-L,L] \times [-L,L]$ and we discretise the computational domain with a uniform mesh as

$ $$\begin{array}{l} \Delta x &=& \frac{L}{2 M}, x_i =i \Delta x, i=-M, \ldots, 0, \ldots, M, \nonumber\\ \Delta y &=& \frac{L}{2 M}, y_j=j \Delta y, j=-M, \ldots, 0, \ldots, M, \nonumber\\ \Delta \theta &=& \frac{2 \pi}{m + 1}, \theta_k=k \Delta \theta, k= 0, \ldots, m. \end{array}$$ $

The phase space densities and the local densities are approximated as

$ f_{1,2}(t,\mathbf{x},\mathbf{v}) \approx f_{1,2}(h \Delta t,x_i,y_j,\theta_k) = f^h_{1,2}(x_i,y_j,\theta_k) $

and

$ n_{3,4}(t,\mathbf{x}) \approx n_{3,4}(h \Delta t, x_i,y_j) = n^h_{3,4}(x_i,y_j), $

where $\Delta t$ is the time step.

We numerically solve the problem defined by the discretised versions of equations (14)-(17), periodic boundary conditions and suitable initial data. Simulations are performed in MATLAB with $L=15$, $M=56$, $m=36$ and $\Delta t=0.014$.

The method for calculating numerical solutions is based on a time-splitting scheme. We begin by updating $f_1^h$ and $f_2^h$ by discretising and solving the homogeneous equations associated to the advection-reaction equations (14) and (15).

A flux-limiting scheme (see for instance [27]) is used to treat the advective terms. First, we advect $f_1^h$ and $f_2^h$ in the $y$-direction, and we denote the results obtained by $f_1^{h+1/4}$ and $f_2^{h+1/4}$. Secondly, we advect $f_1^{h+1/4}$ and $f_2^{h+1/4}$ in the $x$-direction, and we denote the results obtained by $f_1^{h+1/2}$ and $f_2^{h+1/2}$. In more detail, making use of the notations $f_{1,2}^{h}(x_i, y_j,\theta_k) = f^h_{1,2,i,j,k}$ and dropping the indexes $1,2$ for the sake of readability, we compute

$ f^{h+1/4}_{i,j,k}=f^{h}_{i,j,k}-\lambda [F_{i,j,k}^h - F^h_{i,j-1,k} ], $

with $\lambda = \Delta t / \Delta x$ and the flux $F_{i,j,k}^h$ being defined as the combination of a lower order flux (L) and a higher order flux (H), that is,

$ F _{i,j,k}^h = H_{i,j,k}^h - ( 1 - \phi^h_{i,j,k} ) [H_{i,j,k}^h - L^h_{i,j,k}]. $

In the above equation,

$ L^h_{i,j,k} = \sin(\theta_k) \; f^{h}_{i,j,k}, $

and, according to the Richmyer two-step Lax-Wendroff method,

$ H^h_{i,j,k} = \sin(\theta_k) \; f^{h+1/4}_{i,j+1/2,k}, $

with

$ f^{h+1/4}_{i,j+1/2,k} = {1 \over 2}(f^h_{i,j,k} - f^h_{i,j+1,k} ) - {\lambda \over 2} [L^h_{i,j+1,k} - L^h_{i,j,k}]. $

The quantity $\phi^h_{i,j,k}$ is the Ultrabee flux limiter, that is,

$ \phi^h_{i,j,k} = \max(0,\min(1, 2 \; r_{i,j,k}^h),\min(r_{i,j,k}^h, 2)), $

with

$ r^h_{i,j,k}= {{{f^h_{i,j,k}} - f^h_{i,j-1,k}} \over {f_{i,j+1,k}^h - f_{i,j,k}^h}}. $

An analogous procedure is used to compute $f_1^{h+1/2}$ and $f_2^{h+1/2}$.

Then, we update $f_1^{h+1/2}$ and $f_2^{h+1/2}$ by discretising and solving equations (14) and (15) without advective terms on the left-hand sides. The system of equations obtained after the discretisation is solved by means of a fourth-order Runge-Kutta method. In this way, we obtain $f_1^{h+1}$ and $f_2^{h+1}$.

On the other hand, $n_3^{h}$ and $n_4^{h}$ are updated, respectively, through the systems of equations resulting from the discretisation of equations (16) and (17), which are solved by using a fourth-order Runge-Kutta method. Thus, we obtain $n_3^{h+1}$ and $n_4^{h+1}$. A classical second-order centred finite difference scheme is used to calculate the numerical solutions of the reaction-diffusion equation (17).

We consider different initial conditions to mimic different biological scenarios. Moreover, we vary the values of the parameters $\kappa$ and $\mu$, while we keep the other parameters equal to suitable non-zero values. In particular, we set $\nu=1$, $\alpha=0.1$, $\beta=0.1$ and $\gamma=0.1$. Additional simulations show that variations of these parameters, within reasonable ranges, leave the qualitative properties of the emerging patterns unaltered.

---

3.2. Emergence of spots, stripes and hole patterns

We focus on a sample composed of mesenchymal-like cells and chemotactic cytokines only. Cells are initially characterised by a uniform space distribution parametrised by $n^0_2 \in \mathbb{R}_+$, and their velocities are homogeneously distributed. The initial distribution of cytokines is a small positive random perturbation of the zero level. We test the possibility of obtaining different emerging patterns by tuning the value of the initial cell density (i.e., the value of the parameter $n^0_2$). We are interested in the case where the total cell density does not vary over time. Therefore, we neglect the effects of non-conservative phenomena (i.e., we set $\kappa=\mu=0$).

The results of Fig. 1 highlight how increasing values of the parameter $n^0_2 \in (0.005, 0.1)$ lead to the emergence of different spatial patterns for large values of $t$. These patterns are similar to those observed in macroscopic models of chemotaxis [33], and in other transport models [21, 35]. They also exhibit a dependance on the total cell density and on the size of the spatial domain which is analogous to that of classical models for chemotactic movement (see also Remark 1). In more detail:

(ⅰ) if $n^0_2 \in (0.005, 0.02)$, the emergent patterns are highly concentrated spots of aggregation [vid. Fig. 1(A)];

(ⅱ) if $n^0_2 \in (0.02, 0.07)$, the emergent patterns are wider aggregation spots [vid. Fig. 1(B)] or stripes [vid. Fig. 1(C)];

(ⅲ) if $n^0_2 \in (0.07, 0.15)$, the emergent patterns have a hole structure [vid. Fig. 1(D)].

In analogy with classical macroscopic models of chemotaxis, cells tend to aggregate in the local maxima of the chemoattractant (data not shown). However, the quadratic dependence of the consumption rate of cytokines on the local cell density [see definition (3)] seems to prevent finite time blow-up, which is observed in standard macroscopic models. This allows for the formation of bounded aggregation patterns. Additional simulations (data not shown) suggest that relevant patterns cannot emerge when $n^0_2 < 0.005$ (the density is too low for any spatial organisation of cells) and $n^0_2 > 0.15$ (overcrowding occurs).

Remark 1.  The total cell density $N(t)$ and the average local cell density $\bar{\varrho}(t)$ [see equations (2)] are preserved, since we are neglecting the effects of non-conservative phenomena (i.e., we set $\kappa=\mu=0$). In particular, $\bar{\varrho}(t)= n^0_2$ for all $t \geq 0$. As a consequence, higher values of the parameter $n^0_2$ correspond to higher values of the average local cell density. This observation supports the idea that the long-time behaviour of the spatial patterns is governed by the asymptotic values of the average local density of cells. Therefore, as a complementary interpretation of the results presented in this subsection, we note that increasing asymptotic values of the average local cell density may induce the formation of different spatial patterns.

---

3.3. Emergence of spots, stripes and honeycomb patterns

We assume that the sample is initially composed of chemotactic cytokines, mesenchymal-like cells and epithelial-like cells in motion. The distribution of cytokines is a small positive random perturbation of the zero level. Cells are uniformly distributed in space and their distributions are parametrised by $n^0_{1} \in \mathbb{R}_+$ and $n^0_{2} \in \mathbb{R}_+$. Moreover, cellular velocities are homogeneously distributed. We test the possibility of obtaining different emerging patterns by varying the value of the quotient of the average rate of proliferation and the average rate of death due to competition for nutrients (i.e., the value of the ratio $\kappa/\mu$). For simplicity, we set $\kappa=1$ and we tune the value of the parameter $\mu$.

The results in Fig. 2 highlight how decreasing values of the parameter $\mu \in (0.007, 0.2)$ lead, in the limit of large times, to the emergence of different spatial patterns. In fact:

(ⅰ) if $\mu \in (0.05, 0.2)$, the emergent patterns are highly concentrated spots of epithelial-like and mesenchymal-like cells [vid. Fig. 2(A) and Fig. 2(B)];

(ⅱ) if $\mu \in (0.015, 0.05)$, the emergent patterns are stripes of epithelial-like and mesenchymal-like cells [vid. Fig. 2(C) and Fig. 2(D)];

(ⅲ) if $\mu \in (0.007, 0.015)$, the emergent patterns display a honeycomb structure, which results from the segregation between epithelial-like and mesenchymal-like cells [vid. Fig. 2(E) and Fig. 2(F)], an emergent behaviour documented in the biological literature [26, 30].

These results have been obtained with $n^0_{1,2}=0.04$. Additional simulations (data not shown) support the idea that the values of the parameters $n^0_{1}$ and $n^0_{2}$ do not affect the qualitative properties of the patterns shown in Fig. 2. In fact, in agreement with the considerations drawn in Remark 1, the numerical results presented here suggest that the long-time behaviour of the spatial patterns depend on the asymptotic values attained by the average local cell density $\bar{\varrho}(t)$, which can be easily proven to be equal to $\kappa \; (\mu |X|)^{-1}$. In more detail,

(ⅰ) when $\mu=0.1$, $\displaystyle{\lim_{t \to +\infty} \bar{\varrho}(t) \approx 0.01}$ and we observe the emergence of highly concentrated spots [to be compared with case (i) of Subsection 3.2];

(ⅱ) when $\mu=0.02$, $\displaystyle{\lim_{t \to +\infty} \bar{\varrho}(t) \approx 0.05}$ and we observe the emergence of stripes [to be compared with case (ⅱ) of Subsection 3.2];

(ⅲ) when $\mu=0.01$, $\displaystyle{\lim_{t \to +\infty} \bar{\varrho}(t) \approx 0.1}$ and we observe the emergence of segregation patterns with a honeycomb structure [to be compared with case (ⅲ) of Subsection 3.2].

Additional simulations (data not shown) support the conclusion that relevant patterns cannot emerge for $\mu > 0.2$ (the asymptotic value of $\bar{\varrho}(t)$ is too low for any spatial organisation of cells) and $\mu < 0.007$ (the asymptotic value of $\bar{\varrho}(t)$ is too high and overcrowding occurs).

---

3.4. Formation of ring-like patterns

We consider a sample where epithelial-like cells at rest and chemotactic cytokines are not present at time $t=0$. The initial cell distributions are radially symmetric in space, and the cell velocities are homogeneously distributed, that is,

$ f_{1,2}(t=0,x,y, \theta)= C_{1,2} \; e^{-\frac{{x}^2+{y}^2}{20}}, \forall \; (x,y, \theta) \in X \times [0, 2 \pi), \;C_{1,2} \in \mathbb{R_+}. $

(19)

We are interested in the case where the total cell density does not vary over time. Therefore, we neglect the effects of non-conservative phenomena (i.e., we set $\kappa=\mu=0$).

The results presented in Fig. 3 show that, while epithelial-like cells rapidly stop moving because of adhesive interactions [vid. Fig. 3(A)], mesenchymal-like cells diffuse throughout the sample [vid. Fig. 3(B)] and follow the chemotactic path created by the diffusing cytokines [vid. Fig. 3(C)]. The resulting pattern is an expanding ring-like structure made of mesenchymal-like cells, which surrounds a cluster of epithelial-like cells kept at rest by homotypic adhesion.

The results presented here have been obtained with $C_{1,2}=10$. Additional simulations (data not shown) suggest that the values of the parameters $C_{1,2}$ do not influence the qualitative properties of the patterns in Fig. 3. In fact, these patterns result from the interplay between the radial symmetry of the initial cell distributions, the absence of chemotactic cytokines inside the system at time $t=0$, the tendency of epithelial-like cells to stop moving because of homotypic adhesion, and the fact that mesenchymal-like cells spread throughout the sample following the gradient of chemotactic cytokines.

Such ring-like structures of mesenchymal-like cells resemble the fibrous capsules observed in hepatic tumours, which seem to be mainly composed of cells expressing a mesenchymal-like phenotype [25]. They also share some striking similarities with patterns arising from different biological contexts, such as tumour growth [13] and evolution of bacterial populations [38], and with the outcomes of chemotaxis models that incorporate some specific effects related to the finite size of individual cells [35].

---

4. Conclusions and perspectives

We have developed an integro-differential model for the co-culture dynamics of epithelial-like and mesenchymal-like cells. The strategy that we have used to model the consumption of the chemoattractant seems to prevent blow-up in finite time, which is usually observed in classical macroscopic models of chemotaxis. This allows for the formation of bounded aggregation patterns, whose long-time behaviour depends on the asymptotic value of the average local density of cells. In fact, higher asymptotic values of the average local cell density lead to the emergence of different spatial patterns, such as spots, stripes and honeycomb structures. Furthermore, our results support the idea that epithelial-like and mesenchymal-like cells can segregate when there is little competition for nutrients. Finally, we have shown how the interplay between epithelial-cell adhesion and mesenchymal-cell spreading can pave the way for the formation of ring-like structures, which resemble the fibrous capsules frequently observed in hepatic tumours.

Future research will be addressed to refine the current modelling strategies. For instance, a natural improvement of the model would be to define the rate of death due to competition for nutrients as a function of the local cell density. Furthermore, since cell motion implies resource reallocation (i.e., redistribution of energetic resources from proliferation-oriented tasks toward development and maintenance of motility), it might be worth considering different proliferation/death rates for moving cells and cells at rest. From the analytical point of view, it would be interesting to study the ring-like patterns discussed in Subsection 3.4. In this respect, the techniques employed in [38] may prove to be useful. Finally, spatial dynamics play a pivotal role in the evolution of many living complex systems, including biological and social systems (see for instance [32]). Therefore, another possible research direction would be to investigate if the modelling approach presented here could be used profitably to model the dynamics of other living systems.

---