Modelling pattern formation through differential repulsion

TL;DR

This paper develops a two-species particle interaction model to understand cell segregation, deriving a macroscopic description from microscopic rules, and analyzing the conditions for segregation with comparisons to simulations and experiments.

Contribution

It introduces a novel two-species interaction model for cell segregation and derives a macroscopic PDE from microscopic particle dynamics.

Findings

Segregation occurs within specific parameter ranges.

Macroscopic model aligns well with particle simulations.

Results suggest potential links to biological experiments.

Abstract

Motivated by experiments on cell segregation, we present a two-species model of interacting particles, aiming at a quantitative description of this phenomenon. Under precise scaling hypothesis, we derive from the microscopic model a macroscopic one and we analyze it. In particular, we determine the range of parameters for which segregation is expected. We compare our analytical results and numerical simulations of the macroscopic model to direct simulations of the particles, and comment on possible links with experiments.

Tables

Table 1. Table 1: Model parameters for the inter- and intra- species forces. The value of parameter s ∗ superscript 𝑠 s^{*} has been computed numerically from the formula ( 43 ).

	value of $s$	comment
Case I: ${\kappa }^{AA}={\kappa }^{BB}=2$, ${\tilde{\kappa }}^{AB}={\tilde{\kappa }}^{BA}=2$.
IA	0.5	Stable regime ()
IB	4	Unstable regime ($s>s^{\ast }\approx 1.01$)
Case II: ${\kappa }^{AA}={\kappa }^{BB}=2$, .
IIA	0.5	Stable regime ()
IIB	4	Unstable regime ($s>s^{\ast }\approx 1.43$)
Case III: $2={\kappa }^{AA}>{\kappa }^{BB}=1$, ${\tilde{\kappa }}^{AB}={\tilde{\kappa }}^{BA}=2$.
IIIA	0.5	Stable regime ()
IIIB	4	Unstable regime ($s>s^{\ast }\approx 0.72$)
Case IV: $2={\kappa }^{AA}>{\kappa }^{BB}=1$, .
IVA	0.5	Stable regime ()
IVB	4	Unstable regime ($s>s^{\ast }\approx 1.02$)

Table 2. Table 2: Volume fraction of the green family computed on the simulations of FigS. 5 - 6 at equilibrium for the microscopic model (left column) and for the macroscopic model (right column).

Case number	VF green cells (microscopic model)		VF green cells (macroscopic model)
	Homogeneous IC	Front-like IC	Homogeneous IC	Front-like IC
case (IB)	$48.2\%$	$50.0\%$	$49.7\%$	$50.0\%$
case (IIB)	$40.2\%$	$42.3\%$	$38.5\%$	$42.0\%$
case (IIIB)	$40.6\%$	$42.4\%$	$44.0\%$	$46.0\%$
case (IVB)	$34.8\%$	$35.0\%$	$35.9\%$	$38.0\%$

Equations

$$W^A(X^A,X^B)=$$

$$W^B(X^A,X^B)=$$

$$d{X}_i^A=-\mu {\nabla }_{X_i^A}{W}^A(X^A,X^B)dt+\sqrt{2D^A}d{B}_i,\forall i\in \{1,\dots ,N_A\},$$

$$d{X}_i^B=-\mu {\nabla }_{X_{\ell }^B}{W}^B(X^A,X^B)dt+\sqrt{2D^B}d{B}_{\ell },\forall \ell \in \{1,\dots ,N_B\},$$

$$d{X}_i^A=$$

$$\displaystyle+\sum_{k_{3}=1}^{K_{AB}}\nabla_{x_{1}}\Phi^{AB}(X^{A}_{i(k_{3})},X^{B}_{\ell(k_{3})})\delta_{i(k_{3})}(i)\bigg{]}+\sqrt{2D^{A}}dB_{i}$$

$$d{X}_{\ell }^B=$$

$$\displaystyle+\sum_{k_{3}=1}^{K_{AB}}\nabla_{x_{1}}\Phi^{BA}(X^{B}_{\ell(k_{3})},X^{A}_{i(k_{3})})\delta_{\ell(k_{3})}(\ell)\bigg{]}+\sqrt{2D^{B}}dB_{\ell},$$

$$\begin{array}{ccc}\mathrm{Microscopic}\text{leavevmode}\left\{(X_i^A,X_j^B)\right\}& \stackrel{N,K\to \infty }{\to }& \mathrm{Kinetic}\text{leavevmode}\left\{(f^S,g^{ST},h^{ST}{)}_{S,T=A,B}\right\}\\ \downarrow \epsilon \to 0& & \downarrow \epsilon \to 0\\ \mathrm{Averaged}\text{leavevmode}\mathrm{microscopic}\text{leavevmode}\left\{({\tilde{X}}_i^A,{\tilde{X}}_j^B)\right\}& \stackrel{N,K\to \infty }{\to }& \mathrm{Macroscopic}\text{leavevmode}\left\{(f^A,f^B)\right\}\end{array}$$

$$f_N^A(x,t)=\frac{1}{N_A}\sum _{i=1}^{N_A}{\delta }_{X_i^A(t)}(x),f_N^B(x,t)=\frac{1}{N_B}\sum _{\ell =1}^{N_B}{\delta }_{X_{\ell }^B(t)}(x),$$

$$\begin{array}{rl}{g}_N^{AA}(x_1,x_2,t)=& \frac{1}{2N_A}\sum _{k_1=1}^{K^{AA}}{\delta }_{X_{i(k_1)}^A,X_{j(k_1)}^A}(x_1,x_2)+{\delta }_{X_{j(k_1)}^A,X_{i(k_1)}^A}(x_1,x_2),\\ g_N^{BB}(x_1,x_2,t)=& \frac{1}{2N_B}\sum _{k_2=1}^{K^{BB}}{\delta }_{X_{\ell (k_2)}^B,X_{m(k_2)}^B}(x_1,x_2)+{\delta }_{X_{m(k_2)}^B,X_{\ell (k_2)}^B}(x_1,x_2),\end{array}$$

$$g_N^{AB}(x_1,x_2,t)=\frac{1}{N_A}\sum _{k_3=1}^{K^{AB}}{\delta }_{X_{i(k_3)}^A,X_{\ell (k_3)}^B}(x_1,x_2),$$

$$g_N^{BA}(x_1,x_2,t)=\frac{1}{N_B}\sum _{k_3=1}^{K^{AB}}{\delta }_{X_{\ell (k_3)}^B,X_{i(k_3)}^A}(x_1,x_2)=\frac{N_A}{N_B}{g}_N^{AB}(x_2,x_1,t).$$

$$h_N^{AA}(x_1,x_2,t)$$

$$h_N^{BB}(x_1,x_2,t)$$

$$h_N^{AB}(x_1,x_2,t)$$

$$d{A}_{ij}(t)$$

$$d{B}_{ij}(t)$$

$$d{C}_{ij}(t)$$

$$\begin{array}{rl}{\nu }_{d,\epsilon }^{AA}& ={\nu }_d^{AA}{\epsilon }^{-2},{\nu }_{c,N,\epsilon }^{AA}={\nu }_{c,\epsilon }^{AA}{N}_A^{-1}={\nu }_c^{AA}{N}_A^{-1}{\epsilon }^{-2},\\ {\nu }_{d,\epsilon }^{BB}& ={\nu }_d^{BB}{\epsilon }^{-2},{\nu }_{c,N,\epsilon }^{BB}={\nu }_{c,\epsilon }^{BB}{N}_B^{-1}={\nu }_c^{BB}{N}_B^{-1}{\epsilon }^{-2},\\ {\nu }_{d,\epsilon }^{AB}& ={\nu }_d^{AB}{\epsilon }^{-2},{\nu }_{c,N,\epsilon }^{AB}={\nu }_{c,\epsilon }^{AB}{N}_B^{-1}={\nu }_c^{AB}{N}_B^{-1}{\epsilon }^{-2}.\end{array}$$

$$\mathbb{P}(A_{ij}(t)=1)=\frac{\frac{{\nu }_c^{AA}}{N_A}{\chi }_{\{|X_i^A(t)-X_j^A(t)|\le R\}}}{\frac{{\nu }_c^{AA}}{N_A}+{\nu }_d^{AA}}$$

$$\mathbb{P}(B_{ij}(t)=1)=\frac{\frac{{\nu }_c^{BB}}{N_B}{\chi }_{\{|X_i^B(t)-X_j^B(t)|\le R\}}}{\frac{{\nu }_c^{BB}}{N_B}+{\nu }_d^{BB}}$$

$$\mathbb{P}(C_{ij}(t)=1)=\frac{\frac{{\nu }_c^{AB}}{N_B}{\chi }_{\{|X_i^A(t)-X_j^B(t)|\le R\}}}{\frac{{\nu }_c^{AB}}{N_B}+{\nu }_d^{AB}}$$

$$\begin{array}{rl}d{\tilde{X}}_i^A=& -\mu (\frac{1}{N_A}\frac{{\nu }_c^{AA}}{{\nu }_d^{AA}}\sum _{j=1}^{N_A}{\chi }_{\{|{\tilde{X}}_i^A-{\tilde{X}}_j^A|\le R\}}\nabla {\Phi }^{AA}({\tilde{X}}_i^A-{\tilde{X}}_j^A)\\ & +\frac{1}{N_B}\frac{{\nu }_c^{AB}}{{\nu }_d^{AB}}\sum _{j=1}^{N_B}{\chi }_{\{|{\tilde{X}}_i^A-{\tilde{X}}_j^B|\le R\}}\nabla {\Phi }^{AB}({\tilde{X}}_i^A-{\tilde{X}}_j^B))dt+\sqrt{2D^A}d{B}_i^A\end{array}$$

$$\begin{array}{rl}d{\tilde{X}}_i^B=& -\mu (\frac{N_A}{N_B}\frac{1}{N_A}\frac{{\nu }_c^{AB}}{{\nu }_d^{AB}}\sum _{j=1}^{N_A}{\chi }_{\{|{\tilde{X}}_i^B-{\tilde{X}}_j^A|\le R\}}\nabla {\Phi }^{BA}({\tilde{X}}_i^B-{\tilde{X}}_j^A)\\ & +\frac{1}{N_B}\frac{{\nu }_c^{BB}}{{\nu }_d^{BB}}\sum _{j=1}^{N_B}{\chi }_{\{|{\tilde{X}}_i^B-{\tilde{X}}_j^B|\le R\}}\nabla {\Phi }^{BB}({\tilde{X}}_i^B-{\tilde{X}}_j^B))dt+\sqrt{2D^B}d{B}_i^B.\end{array}$$

$$\begin{split}\partial_{t}f^{A}&=D^{A}\Delta_{x}f^{A}+\nabla_{x}\cdot\bigg{(}f^{A}(x,t)\nabla_{x}\big{(}\tilde{\Phi}^{AA}\ast f^{A}\big{)}(x,t)\bigg{)}+\nabla_{x}\cdot\bigg{(}f^{A}(x,t)\nabla_{x}\big{(}\tilde{\Phi}^{AB}\ast f^{B}\big{)}(x,t)\bigg{)},\\ \partial_{t}f^{B}&=D^{B}\Delta_{x}f^{B}+\nabla_{x}\cdot\bigg{(}f^{B}(x,t)\nabla_{x}\big{(}\tilde{\Phi}^{BB}\ast f^{B}\big{)}(x,t)\bigg{)}+\nabla_{x}\cdot\bigg{(}f^{B}(x,t)\nabla_{x}\big{(}\tilde{\Phi}^{BA}\ast f^{A}\big{)}(x,t)\bigg{)},\end{split}$$

$${\tilde{\Phi }}^{AA}(x)$$

$${\tilde{\Phi }}^{BB}(x)$$

$${\tilde{\Phi }}^{AB}(x)$$

$${\tilde{\Phi }}^{BA}(x)$$

$${\partial }_t{f}^A$$

$${\partial }_t{f}^B$$

$$f^A=f_{\ast }^A+{\tilde{f}}^A,f^B=f_{\ast }^B+{\tilde{f}}^B,$$

$$\begin{split}\partial_{t}f^{A}&=f_{*}^{A}\Delta\big{[}f^{A}\ast\tilde{\Phi}^{AA}+f^{B}\ast\tilde{\Phi}^{AB}\big{]}+D^{A}f^{A}\\ \partial_{t}f^{B}&=f_{*}^{B}\Delta\big{[}f^{B}\ast\tilde{\Phi}^{BB}+f^{A}\ast\tilde{\Phi}^{BA}\big{]}+D^{B}f^{B}.\end{split}$$

$$\hat{F}(y)=\frac{1}{2\pi }{\int }_{{\mathbb{R}}^2}{\exp }^{-ix\cdot y}F(x)dx.$$

$$\begin{split}\partial_{t}\begin{pmatrix}\hat{f}^{A}\\ \hat{f}^{B}\end{pmatrix}(y,t)&=\begin{pmatrix}-f_{*}^{A}|y|^{2}\big{(}2\pi\hat{\tilde{\Phi}}^{AA}(y)+\frac{D^{A}}{f_{*}^{A}}\big{)}&-f_{*}^{A}|y|^{2}2\pi\hat{\tilde{\Phi}}^{AB}(y)\\ -f_{*}^{B}|y|^{2}2\pi\hat{\tilde{\Phi}}^{BA}(y)&-f_{*}^{B}|y|^{2}\big{(}2\pi\hat{\tilde{\Phi}}^{BB}(y)+\frac{D^{B}}{f_{*}^{B}}\big{)}\end{pmatrix}\begin{pmatrix}\hat{f}^{A}\\ \hat{f}^{B}\end{pmatrix}(y,t)\\ &:=M(y)\begin{pmatrix}\hat{f}^{A}\\ \hat{f}^{B}\end{pmatrix}(y,t).\end{split}$$

$$\left(\begin{array}{c}{\hat{f}}^A\\ {\hat{f}}^B\end{array}\right)(y,t)=c_1(y){\exp }^{{\lambda }_1(y)t}{\vec{u}}_1(y)+c_2(y){\exp }^{{\lambda }_2(y)t}{\vec{u}}_2(y),$$

$$\Delta(M)=|y|^{4}f_{*}^{A}f_{*}^{B}\bigg{[}\left(2\pi\hat{\tilde{\Phi}}^{AA}+\frac{D^{A}}{f_{*}^{A}}\right)\left(2\pi\hat{\tilde{\Phi}}^{BB}+\frac{D^{B}}{f_{*}^{B}}\right)-(2\pi)^{2}\hat{\tilde{\Phi}}^{AB}\hat{\tilde{\Phi}}^{BA}\bigg{]},$$

Full text

1. Institut Denis Poisson, Université d’Orléans, CNRS, Université de Tours.

B.P. 6759, 45067 Orléans cedex 2, France.

2. Institut Universitaire de France, Paris, France.

3. Department of Mathematics, Imperial College London,

London SW7 2AZ, United Kingdom.

4. INRIA Team Mamba, INRIA Paris, 2 rue Simone Iff, CS 42112, 75589 Paris, France

Université Pierre et Marie Curie-Paris 6, UMR 7598 LJLL, BC187

4, Place de Jussieu, F-75252 Paris Cedex 5, France.

5. Department of Mathematics, University College London

Gower Street London WC1E 6BT, UK, United Kingdom.

1 Introduction

The organisation of biological tissues during development is accompanied by the formation of sharp borders between distinct cell populations. During the morphogenesis of numerous tissues/organs, cells of the same type regroup into regions, creating niches with specific identities that drive the differentiation of particular cell types. This spatial organization is ensured via cell-cell signalling leading specific cells/tissues to form at the appropriate location. The maintenance of this cell segregation is key in adult tissue homeostatis, and its disruption can lead tumor cells to spread and form metastasis. This segregation is challenged during tissue growth and morphogenesis due to the high mobility of many cells that can lead to intermingling. Therefore, understanding the mechanisms involved in the generation and maintenance of cell segregation is of tremendous importance in tissue morphogenesis, homeostasis, and in the development of various invasive diseases such as tumors.

Numerous experiments have been conducted to identify the mechanisms of cell segregation. Experiments show that mixing cells from different tissues in vitro leads to their segregation, with initially fuzzy borders that sharpen in time [53]. This has been observed for many systems for instance in the development of the wing imaginal disc in Drosophilia [33], in the developing nervous system [28] etc. So far, three types of mechanisms have been identified to have a role in segregation and border formation, namely (i) differential adhesion, (ii) contact cell repulsion and (iii) cortical tension. Indeed, the segregation of cells derived from different tissues in vitro was initially suggested to occur through a combination of directed cell migration and the selective adhesion of cells of the same type [48]. However, other studies have shown that the contact inhibition of cell migration induced by contact repulsion of cells by Eph receptor and ephrin signaling, and finally the induction of cortical tension by actomyosin contraction are important mechanisms that can restrict intermingling between cell populations. However, it remains unclear whether individually each of these mechanisms account for cell segregation or border sharpening, and to what extent the interplay of these different mechanisms is required to achieve cell segregation.

Several modelling efforts have been done to identify the main mechanisms involved in cell segregation. The mathematical models for cell-cell interactions are usually agent-based models, where each cell undergoes a random walk exclusion process and interacts with its neighbours. For instance the cellular Potts model for segregation between two cell populations [32] predicts cell rearrangements in epithelia based on the minimization of a free energy and has been widely used to explore the rate of cell sorting due to differential adhesion, but does not include cell migration as a mechanism. The varying adhesive and repulsive forces between different cell populations, which can result from Eph/ephrin interactions, have also been modelled by representing cells as spheres which can attract or repel each other, giving rise to empirically observed cell sorting pattern [53]. In [1], the authors develop a mathematical model for Eph/ephrin regulated cell-cell segregation and tissue boundary formation, which features independent random cell motion, Eph/ephrin-dependent attraction/repulsion interactions between neighbouring cells, and cell division. The authors in [1] show that the dynamics of Eph/ephrin-mediated cell cluster formation and cell segregation can be captured with these mechanisms. However, this model do not use parameters from measurements of cell behaviour or examine whether cell repulsion is sufficient for border sharpening. Another approach, which simulates cell adhesion, de-adhesion and migration in greater detail [51, 52], was used to model the time course of segregation of cells differing in cadherin expression. The results from this model were acurate for describing cell segregation mediated by differential expression of cadherins but less accurate when simulating the significantly faster rate of Eph-ephrin mediated cell segregation. However when modified to account for repulsive behaviors [53], the model correctly reproduced the experiments and showed that heterotypic repulsion can account for cell segregation and border sharpening, and is more efficient than decreased heterotypic adhesion.

All these results suggest that cell segregation and border sharpening is the result of a complex interplay between homotypic/heterotypic cell adhesion, de-adhesion and repulsion. But how the balance of these phenomena is precisely linked to the existence/size of the segregated zones remains unclear to this day. In this paper, we aim to provide a mathematical framework which enables to quantitatively link the segregation and border sharpening ability of the tissue to these cell-cell interaction phenomena of interest. As agent-based models do not enable precise mathematical analysis of their solutions due to the lack of theoretical results, we turn towards a continuous -macroscopic- model for which the theoretical study gives precise criteria for phase transitions as functions of key model parameters. As a drawback contrary to microscopic models, macroscopic models lose the information at the individual level. In order to overcome this weakness, we aim to derive, as rigorously as possible, the macroscopic model from an agent-based formulation to ensure the good correspondence between the two formulations as it was done in [7, 26, 27].

The starting point is an individual-based model inspired from [7, 8] and which bears similarities with the approach [53]. We consider two families of cells, each cell being modelled as a point particle which interacts with its close neighbors via local cross-links. The links are modeled by springs that are randomly created and destructed. This enables us to model cell-cell attraction and repulsion, with different spring strengths according to the type of link (intraspecies or interspecies). We let the particles move randomly in space to model random motion of cells. In the mean field limit, assuming large numbers of particles and links as well as propagation of chaos, the corresponding kinetic system consists of two equations for the individual particle distribution functions and two equations for the link densities. In the large-scale limit and in the regime where the link creation/destruction frequency is very large, it was shown in [27, 8, 7] that the link density distributions become local functions of the particle distributions. The latter evolve through aggregation diffusion equations. Similar macroscopic model, but with nonlinlinear porous medium type of diffusion has been recently considered analytically and numerically in [19]. Although the derivation of the macroscopic model in the hydrodynamic limit follows closely the steps of [8], the originality of this work lies in the presence of two coupled families of cells which introduce a new level of complexity and make the stability analysis more involved than in previous works. Inspired from the results of [53], we mainly consider repulsive springs and aim to quantify the influence of heterotypic/homotypic repulsion on cell segregation and border sharpening. By addressing the stability of a homogeneous distribution of particles for Hookean repulsive potentials, we obtain a precise condition for the phase transition, which links the system segregation ability to the model parameters and give further insight into the cell segregation processes.

Our study shows that in a system composed of two-species repelling each other, the interspecies forces must be large enough to compensate both for the diffusion and for the intra-species repulsion, which both tend to homogeneize the system. Aggregation will therefore be ensured if and only if interspecies repulsion wins over diffusion and intraspecies repulsion. In the case where attractive interactions are considered, we have noted that a necessary condition for the aggregation of the species (or equivalently instability of the homogeneous steady-state) is that the interspecies forces are of the same sign. To observe aggregates, the two families must therefore either repulse or attract each other, but must have the same effect on each other. On the contrary, if one family is attracted by the other and the other repulses it, we will always observe a homogeneous distribution at equilibrium (intermingling of the two families). A third remark concerns the size of the clusters when aggregation occurs. As the interspecies repulsion force increases, the particles of a given family aggregate more together, leading to a decrease of the size of the local aggregates of the compressed family. These conclusions were confirmed by simulations performed for both the microscopic and macroscopic models and are pictured in Fig. 1. Numerical simulations show that both the micro- and macro- models are in excellent agreement with the predictions of the stability analysis performed on the continuous model. The quantitative agreement obtained between both models show that the macroscopic model is a good approximation of the microscopic model as the number of individuals goes to infinity, provided the interspecies repulsion forces are not too large. Indeed for large interspecies repulsion forces, we find some structural discrepancies between the two models, where the microscopic dynamics seems to favor the formation of rounder aggregates compared to the elongated structures obtained with the macroscopic model. We find that the microscopic dynamics is comprised of two time phases: the first phase consists of a fast segregation between the two families and is followed later by a reorganisation of the clusters which get rounder at large times. The macroscopic dynamics does not seem to contain the second phase (restructuring), suggesting that this phenomenon can be due to finite size effects. We postulate that these discrepancies come from the microscopic noise due to thermal fluctuations, which gives rise to instabilities and allows the agent-based system to reach new states which are not available in the deterministic description, or produce spatial correlations which in turn dominate the macroscopic system behavior. Several works have reported these phenomena [41, 37, 4, 14], particularly at onset for transitions from metastable or unstable phases, in which microscopic noise can be amplified to macroscopic time and length scales. The exploration of these effects will be the subject of future works. Finally, we remark that the segregation process is seen in the numerical simulations to be efficient even close to the instability threshold: as soon as the homogeneous state becomes unstable, the system evolves towards a well segregated configuration. This suggests the presence of a subcritical bifurcation [20, 7, 18], the study of which we also leave for future work.

The paper is organized as follows. In Section 2, we first give the main ingredients of the microscopic model, and then we sketch the main steps of the derivation of the macroscopic dynamics. We present two approaches that involve taking limit of large number of individuals and large scale/fast network remodelling limit of the microscopic model in the different orders; the details of one of them are moved to Appendixes A and B. Section 3 is devoted to the stability analysis around the homogeneous steady-states of the macroscopic model: in Section 3.1, we give the stability results in the whole space, Section 3.2 explores the case of periodic boundary conditions and Section 3.3 is devoted to the case of phase separated initial conditions. Finally, Section 4 presents the numerical results, performed on the microscopic and macroscopic models in different regime of parameters, with a particular focus on the qualitative and quantitative comparison between the two models.

2 Mathematical modelling

2.1 Microscopic model

The model features two families of particles referred to as type $A$ and type $B$. Each particle can link/unlink with neighbors located in a ball of radius $R$ from its center. Each particle can link with a neighbor of its own family as well as with a neighbor of the other family, with no restriction on the number of links per particle and with the same detection radius $R$ no matter the type of link (intra- or inter- species). In order to model tissue plasticity, the links are not permanent but supposed to be created and suppressed via random processes. In this way, the model allows for constant remodelling of the link network. Each link between two particles generates a spring-like interaction potential, which depends on the link type (intra- or inter- species link). Finally, particle positions are subject to random positional noise to model the movements of the tissue.

In this paper, we restrict ourselves to a two-dimensional model. We consider a set of $N_{A}$ particles of type $A$ and $N_{B}$ particles of type $B$ described by their centers $(X^{A}_{i},X^{B}_{\ell})\in\mathbb{R}^{2}\times\mathbb{R}^{2}$, $i\in[1,N_{A}],\ell\in[1,N_{B}]$ respectively. The link creation and suppression are supposed to follow Poisson processes of frequencies $\nu_{c,N,\varepsilon}^{AA},\nu_{c,N,\varepsilon}^{AB},\nu_{c,N,\varepsilon}^{BB}$ and $\nu_{d,\varepsilon}^{AA},\nu_{d,\varepsilon}^{AB},\nu_{d,\varepsilon}^{BB}$, where the subscripts $c$ and $d$ refer to ’creation’ and ’deletion’ respectively, and the superscripts $AA,BB$ and $AB$ denote intraspecies links (AA, BB) and interspecies links (AB); $\varepsilon$ is a scaling parameter, and the subscripts $\varepsilon$, and $N$ for the $\nu_{c}$’s make the dependency of these rates on $N_{A}$, $N_{B}$ and $\varepsilon$ explicit, as will be explained below. We suppose that the intraspecies links generate pairwise symmetric potentials $\Phi^{AA}(X^{A}_{i},X^{A}_{j})$ and $\Phi^{BB}(X^{B}_{\ell},X^{B}_{m})$, not necessarily equal, and that the interspecies links generate non-symmetric potentials $\Phi^{AB}(X^{A}_{i},X^{B}_{\ell})\neq\Phi^{BA}(X^{B}_{\ell},X^{A}_{i})$, modelling the fact that the two particle families act differently on each other. For the moment we do not specify interaction potentials, trying to keep the derivation at maximal level of generality. Note that $\Phi^{AB}$ refers to the action a type B particle exerts on a type $A$ particle while $\Phi^{BA}$ is the action a type $A$ particle exerts on a particle of type $B$. We define the total ”energy” $W^{A}$ of the $A$-particles as the sum over all pairwise link potentials acting on particles $A$, and $W^{B}$ is the sum over all pairwise link potentials acting on particles $B$:

[TABLE]

where $K^{AA},K^{BB},K^{AB}$ denote the (time-dependent) total number of links between particles of type $A$, particles of type $B$, and interspecies links respectively. In the formulas above $(i(k_{1}),j(k_{1}))$ denote the indices of particles of type $A$ connected by the intraspecies link $k_{1}$, $(\ell(k_{2}),m(k_{2}))$ the indices of particles of type $B$ connected by link $k_{2}$. By a slight abuse of notation we denote by $(i(k_{3}),\ell(k_{3}))$ the indices of particles of type A connected to particles of type B by link $k_{3}$.

Particle motion during a time interval between two linking/unlinking events is supposed to occur in the so-called overdamped regime. The resulting equations contain a drift term in the steepest descent direction of the ”energies” $W^{A}$ and $W^{B}$ and a noise term:

[TABLE]

where $\mu>0$ is the mobility coefficient considered to be given, and $B_{i}$ is a 2-dimensional Brownian motion $B_{i}=(B_{i}^{1},B_{i}^{2})$ of intensity $D^{A}>0$ for species $A$ and $D^{B}>0$ for species $B$. Inserting (1)-(2) into (3)-(4), we obtain

[TABLE]

where $\delta_{i}(j)$ stands for the Kronecker delta.

This model bears similarities with the works of [51, 52, 53]. The main difference lies in the fact that cells are modelled as individual spheres here, while in [53] each cell is supposed to be composed of a set of several spheres maintained in a ring of a given radius. Therefore we do not take into account the role of cell deformation in this paper. However, this simpler modelling choice enables us to reduce the complexity of the system and to derive a macroscopic model as performed in the next section. The influence of this modelling simplification on the result will be the subject of future works. In the next section, we present the macroscopic model obtained in the limit of a large number of individuals and in the fast linking-unlinking process.

2.2 Macroscopic model

The derivation of a macroscopic model from the microscopic model defined by Eqs. (5)-(6) requires two limits: (i) limit of large number of individuals and large number of links, denoted $N,K\to\infty$ and (ii) a large scale or fast network remodelling limit, denoted $\varepsilon\to 0$.

Performing the limit (i) first yields, after a mean-field assumption, a kinetic system, from which the macroscopic dynamics can be derived in the $\varepsilon\to 0$ limit. This approach will be referred to as the Approach I.

Performing the limit (ii) first yields, after averaging over the network configurations, an effective dynamics for the particles, from which the same macroscopic dynamics can be derived, in the $N,K\to\infty$ limit. This approach will be referred to as the Approach II.

This structure is summarized in the diagram:

[TABLE]

The meaning of the $f^{S},g^{ST},h^{ST},\tilde{X}^{S}$ is explained below. Our final goal is a macroscopic model describing the evolution in time of the particle distributions $f^{A}(x,t)$ and $f^{B}(x,t)$ of the type-$A$ particles and type-$B$ particles respectively. Let us roughly discuss the two possible approaches mentioned above.

Sketch of Approach I

For finite $N_{A},N_{B}$ we define:

[TABLE]

where $\delta_{X^{S}_{i}(t)}(x)$ denotes the Dirac delta located at $X^{S}_{i}(t)$ for $S$ being either $A$ or $B$. In the large $N_{S}$ limit it gives the probability to find a particle of type $S$ at point $x$ at time $t$.

To write the kinetic model, we need to define the (symmetric) empirical measures $g^{SS}_{N}(x_{1},x_{2},t)$ of the intraspecies links ($S$ being either A or B) by:

[TABLE]

with a similar definition of the Dirac deltas. Such a $g_{N}^{{SS}}(x_{1},x_{2},t)$ gives in the large $N_{S}$ limit the density of links connecting a particle of a given type and located within a volume $dx_{1}$ about $x_{1}$ with a particle of the same type located within a volume $dx_{2}$ about $x_{2}$, normalized by $N_{S}$; note the integral of $g^{SS}_{N}$ is not $1$. As will become clear below, we will be interested in a regime where $K_{SS}$ and $N_{S}$ have the same order of magnitude; hence the chosen normalization ensures that $g^{SS}_{N}$ is of order $1$. Now, we define a non-symmetric empirical measures for the interspecies links $g^{AB}_{N}(x_{1},x_{2},t)$:

[TABLE]

and analogously

[TABLE]

Note that $g^{AB}_{N}(x_{1},x_{2},t)$ gives in the large $N_{A}$ limit the density of interspecies links between a particle of type A located within a volume $dx_{1}$ about $x_{1}$ and a particle of type B located within a volume $dx_{2}$ about $x_{2}$, normalized by $N_{A}$.

Finally, the derivation of the kinetic model also involves the two-particle distribution functions defined by:

[TABLE]

Here, $h^{AA}_{N}$ and $h^{BB}_{N}$ give in the large $N_{A},N_{B}$ limit the probabilities of finding pairs of not necessarily linked particles of the same species around $x_{1}$ and $x_{2}$, while $h^{AB}_{N}(x_{1},x_{2},t)$ gives the probability of finding a particle of type A around $x_{1}$ and a particle of type B around $x_{2}$.

The kinetic system provides evolution equations for $f^{A},f^{B},g^{AA},g^{BB},g^{AB}$, the large $N$-limits of the corresponding empirical densities defined above. Since this derivation follows closely the works of [7, 27] adapted to a two species system, we leave the details in Appendix A. The fast network remodelling limit taken on this kinetic system then formally yields a macroscopic system of evolution equations involving only $f^{A}$ and $f^{B}$, i.e. the macroscopic evolution we are looking for, given by (24), below.

Sketch of Approach II

We denote by $A_{ij}(t),B_{ij}(t),C_{ij}(t)$ the adjacency matrices of particles $A$, $B$, and cross-links $A-B$ respectively. In particular, for $i,j\in\{1,\ldots,N_{A}\}$, $A_{ij}(t)=1$ (resp. $=0$) if particles of type $A$ $i$ and $j$ are connected at time $t$ (resp. not connected). The definition of matrix $B$ is similar. For $i\in\{1,\ldots,N_{A}\}$, $j\in\{1,\ldots,N_{B}\}$, $C_{ij}(t)=1$ (resp. $=0$) if particle $i$ of type $A$ and particle $j$ of type $B$ are connected at time $t$ (resp. not connected). $A$ and $B$ are square symmetric matrices, and $C$ is an $N_{A}\times N_{B}$ rectangular matrix.

The derivation of the reduced microscopic model relies on averaging. The diffusions of particles positions $X_{i}^{A},X_{j}^{B}(t)$ are slow processes, and the links $A_{ij}(t),B_{ij}(t),C_{ij}(t)$ are fast processes: they quickly converge to stationary measures which depend on $X_{i}^{A}(t)$ and $X_{i}^{B}(t)$. We will then compute the evolution of $X_{i}^{A}(t)$ and $X_{i}^{B}(t)$ by averaging the basic dynamical equations (3)-(4) over these stationary measures of the link processes.

The process for the links is written

[TABLE]

where the $N^{AA,d}_{ij},N^{AA,c}_{ij},N^{BB,d}_{ij},N^{BB,c}_{ij},N^{AB,d}_{ij},N^{AB,c}_{ij},$ are independent Poisson processes with rates $\nu_{d,\varepsilon}^{ST}$ for destruction of the link connecting particles of type $S$ and $T$, $S,T\in\{A,B\}$, and $\nu^{ST}_{c,N,\varepsilon}$ for creation of a link between particles of type $S$ and $T$. We will moreover consider the following scaling of these rates

[TABLE]

The subscript $\varepsilon$ (resp. $N_{A},N_{B}$) signals a dependency on $\varepsilon$ (resp. $N_{A},N_{B}$), and the rates without $N_{A},N_{B},\varepsilon$ subscripts are assumed to be independent of $N_{A},N_{B},\varepsilon$. Relations (18) make the scaling with $\varepsilon$ and $N_{A},N_{B}$ needed to perform the limit procedures in (7) explicit. In particular, the scaling with $N_{A},N_{B}$ ensures that the connectivity of any particle remains of order $1$, and the scaling with $\varepsilon$ controls the speed of the linking/unlinking process.

Conditionally on the positions $X^{A}_{i},X^{B}_{j}$, all the processes $A_{ij},B_{ij},C_{ij}$ are independent. The stationary measures of (15)-(16)-(17), for fixed positions $X^{A}_{i},X^{B}_{j}$ are then simply product of Bernoulli measures ($\mathbb{P}$ denotes the probability):

[TABLE]

For $N_{A},N_{B}$ large, the above expressions simplify as the $O(1/N_{A},1/N_{B})$ terms in the denominators are negligible. One can write the equations for the positions, averaged over the stationary measure for the links; calling $\tilde{X}_{i}^{A},\tilde{X}_{j}^{B}$ these new processes, we obtain (neglecting terms of order $1/N_{A},1/N_{B}$)

[TABLE]

[TABLE]

Here we tacitly assumed the translation invariance of the potential $\Phi^{S}T(X_{i},X_{j})=\Phi^{S}T(X_{i}-X_{j})$, which can be relaxed, see the Appendix. From Eqs. (22)-(23), one can deduce the following proposition, describing the the dynamics of the particles density $f^{A}$ and $f^{B}$ in the large $N_{A},N_{B}$ limit:

Proposition 1

Assume Eqs. (22)-(23) and that the potentials are radially symmetric $\Phi^{ST}(X_{i}-X_{j})=\Phi^{ST}(|X_{i}-X_{j}|$; then in the limit $N_{A},N_{B}\to\infty,N_{A}/N_{B}\to r_{AB}>0$, the one particle distribution functions $f^{A}$ and $f^{B}$ are solution of the system

[TABLE]

*where the potentials $\tilde{\Phi}^{ST}$ are given by: *

[TABLE]

and $\ast$ denotes the convolution.

An alternate proof of this proposition based on Approach I is presented in the Appendix. The Hookean potential considered in this paper satisfies the radial symmetry assumption. However, for purposes of derivation of the macroscopic model, this assumption can be relaxed at the expense of more complex formulas (see the Appendix).

The macroscopic model consists of an aggregation-diffusion equations with nonlocal terms, where each particle interacts with its close neighbors of the same family (second term of the right hand side of Eq. (24)) as well as with the ones of the other family (third term of the right hand side of Eq. (24)), and where the diffusive term corresponds to the Brownian motion of individual particles. In the next section, we perform the linear stability analysis to identify phase transitions of the homogeneous steady-state.

3 Analysis of the macroscopic system in the whole space

3.1 Linear stability in the whole space

In this section, we perform a linear stability analysis of the macroscopic model. We recall the macroscopic equations for $f^{A},f^{B}$:

[TABLE]

where the factors $\frac{\nu^{ST}_{c}}{\nu^{ST}_{d}}$ and $r_{AB}$ have been included in the potential functions $\tilde{\Phi}^{ST}$. We first linearize around the homogeneous steady states, i.e $f_{*}^{A}=const.$, $f_{*}^{B}=const.$. Writing

[TABLE]

we have , at the first order:

[TABLE]

For any integrable function $F$ and a vector $y\in\mathbb{R}^{2}$ we recall the definition of the spatial Fourier transform

[TABLE]

Applying it to both sides of (29), we obtain the following system:

[TABLE]

Therefore,

[TABLE]

where $\lambda_{1}(y),\lambda_{2}(y)$ are the eigenvalues of the matrix $M(y)$ and $\vec{u}_{1}(y),\vec{u}_{2}(y)$ are the corresponding eigenvectors.

General case

In the general case, the homogeneous steady state will be unstable if at least one of the eigenvalues of matrix $M(y)$ is positive. By computing the determinant of M:

[TABLE]

we can see that, for general interaction potentials, the constant steady states will be unstable if one of the two conditions is met:

• •

1. $\Delta(M)<0$: $(2\pi\hat{\tilde{\Phi}}^{AA}+\frac{D^{A}}{f_{*}^{A}})(2\pi\hat{\tilde{\Phi}}^{BB}+\frac{D^{B}}{f_{*}^{B}})<4\pi^{2}\hat{\tilde{\Phi}}^{AB}\hat{\tilde{\Phi}}^{BA}$

• •

2. $\Delta(M)\geq 0$ and $Tr(M)>0$: $\begin{cases}(2\pi\hat{\tilde{\Phi}}^{AA}+\frac{D^{A}}{f_{*}^{A}})(2\pi\hat{\tilde{\Phi}}^{BB}+\frac{D^{B}}{f_{*}^{B}})\geq 4\pi^{2}\hat{\tilde{\Phi}}^{AB}\hat{\tilde{\Phi}}^{BA}\\ f_{*}^{A}(2\pi\hat{\tilde{\Phi}}^{AA}+\frac{D^{A}}{f_{*}^{A}})+f_{*}^{B}(2\pi\hat{\tilde{\Phi}}^{BB}+\frac{D^{B}}{f_{*}^{B}})<0.\end{cases}$

Hookean interaction potentials

In order to explicit these conditions as functions of the model parameters, we now suppose that the intra- and inter- species links act as springs of rest length $R$ between the particles. As the detection radius for the interaction is also $R$, this amounts to consider that particles only repulse each other up until distance $R$. To keep enough generality, we consider different interaction intensities between the type A and type B intra- and inter- species springs:

[TABLE]

and not necessarily mutually equal. We have:

[TABLE]

Remembering that the factors $\frac{\nu^{ST}_{c}}{\nu^{ST}_{d}}$ were included in the potential functions $\tilde{\Phi}^{ST}$ and following derivation (70)-(73) from the Appendix B.2, we may introduce:

[TABLE]

which ensures that the equations for $f^{A},\,f^{B}$ are of the form (24) and the linearized system is of the form (29). The form of the new potential is plotted in Fig. 2.

We can now compute the Fourier transform of $\tilde{\Phi}^{ST}$. Using the radial symmetry of $\tilde{\Phi}^{ST}$ and denoting $\Upsilon=|y|$ we can show that

[TABLE]

where $J_{0},J_{1}$ are the Bessel functions of the first kind of order 0 and 1:

[TABLE]

and $H_{0},H_{1}$ are the Struve functions defined by:

[TABLE]

We refer the reader to [7] for the computation of these terms. Replacing the Fourier transforms of the potentials by their expressions in $M(y)$, we can write

[TABLE]

with $c^{ST}=\frac{2\pi\kappa^{ST}f_{*}^{S}\nu^{ST}_{c}R^{2}}{\nu^{ST}_{d}}$ for all $S,T\in\{A,B\}$ and

[TABLE]

Writing $z=|y|R$, the determinant of $M$ can now be written:

[TABLE]

where $c^{\prime ST}=R^{2}c^{ST}=\frac{2\pi\kappa^{ST}f_{*}^{S}\nu^{ST}_{c}R^{4}}{\nu^{ST}_{d}}$ and

[TABLE]

Now, lengthy but straightforward computations show that for $z$ close to the origin we have:

[TABLE]

and so, close to the origin $z=0$, we have

[TABLE]

In order to simplify the analysis, we suppose the following hypothesis:

Hypothesis 1

The intraspecies links generate repulsive potentials, i.e $\kappa^{AA},\kappa^{BB}>0$.

We first note that under Hypothesis 1, the trace of $M$

[TABLE]

is negative for small $z$. Therefore, the homogeneous steady states will be unstable for small $z$ only if $\Delta(M)\leq 0$. Note also that the parameters $c^{\prime AB}$ and $c^{\prime BA}$ should have the same sign to allow the determinant of $M$ to be negative, otherwise the homogeneous steady state will be a stable case. Now, we scale the interspecies link potential intensities with a parameter $s\in\mathbb{R}$ such that $\kappa^{AB}=s\tilde{\kappa}^{AB},\kappa^{BA}=s\tilde{\kappa}^{BA}$. It corresponds to the relevant scaling of the parameters $c^{\prime AB}$ and $c^{\prime BA}$, to simplify the notation we denote the corresponding reference values by the same symbols $c^{\prime AB}$ and $c^{\prime BA}$, then

[TABLE]

and we immediately note that $s$ must be large enough to allow $\Delta(M)$ to be negative. More precisely, the two eigenvalues of $M$ are written:

[TABLE]

where $C_{+}=D^{A}+D^{B}+\frac{c^{\prime AA}+c^{\prime BB}}{24}$ and $C_{-}=D^{A}-D^{B}+\frac{c^{\prime AA}-c^{\prime BB}}{24}$. We can therefore plot their values as functions of $s$ near $z=0$ (see Fig. 3 for $z=0.1$ and parameter values $D^{A}=D^{B}=1,c^{\prime AA}=c^{\prime AB}=c^{\prime BA}=1,c^{\prime BB}=10$).

As shown by Fig. 3, there exists a critical value $s^{*}$ of $s$ such that for $s>s^{*}$ the homogeneous state will be unstable (i.e we will observe cell aggregates). From the definition of $s$, this means that the interspecies repulsion force must be large enough to compensate the intraspecie repulsion and diffusion and enable the two species to separate into clusters. These first results are in accordance with the observations of [53]. Note that, by equating (33) to [math], we can directly compute the value of $s^{*}$ as a function of the model parameters:

[TABLE]

In Fig. 4, we aim to plot the values of $\lambda_{1}(z),\lambda_{2}(z)$ in the unstable regime $s>s^{*}$, $s^{*}$ being determined on Fig.3 (critical value of $s$ such that $\lambda_{1}(z\approx 0)$ becomes positive). We select four values of $s>s^{*}$ and for each of them we compute the functions $\lambda_{1}(z),\lambda_{2}(z)$ using matrix $M(z)$ before Taylor expanding it near 0. The Bessel and Struve functions are approximated numerically.

As one can see in Fig. 4 (I), $\lambda_{1}$ is an increasing function of $s$ at fixed $z$. Furthermore, the most unstable parameter $z$, i.e. the value $z^{*}$ for which $\lambda_{1}(z)$ reaches its maximal value, defined by:

[TABLE]

increases with $s$ (see Fig. 4 (III)). Hence, at small times, the instability should be dominated by Fourier modes with parameter around $z^{\ast}(s)$. Assuming this remains qualitatively true in the time asymptotic regime, one then expects that larger values of $s$ will lead to more aggregated (smaller) patterns.

3.2 Linear stability in the periodic box

For the sake of numerical simulations, we now interpret the above results in the case of a space-periodic domain. In practice, instead of the whole plane we consider the two-dimensional square periodic domain $[-L,L]\times[-L,L]$, however, for $L$ sufficiently large the results change quantitatively but not qualitatively. Introducing the shorthand notation for the Fourier modes

[TABLE]

we may write for $S=A,B$:

[TABLE]

where the Fourier coefficients $\hat{f}^{S}_{k_{1},k_{2}}$ are given by

[TABLE]

The periodic analog of the macroscopic system after Fourier transform is (see (30) for comparison):

[TABLE]

Since $f^{A},f^{B}$ are both probability measures, we take

[TABLE]

which means that the matrix $M$ from (38) has now the form

[TABLE]

This form of the system can be again studied for general potentials. Here we immediately focus on the Hookean case (31), for which the Fourier transform equals

[TABLE]

where we denoted

[TABLE]

As explained in [7], due to specific shape of the potential, we only need to check the stability of the first mode $(k_{1},k_{2})=\pm(1,0)$ or $(k_{1},k_{2})=\pm(0,1)$ in order to find whether the whole system is stable. This is the point where the analysis in the space-periodic domain differs from the whole space case. Note that in (32) we took $z\to 0$, which cannnot be assumed for the discrete values (41). Exactly as before, we compute:

[TABLE]

[TABLE]

Recall however, that assuming Hypothesis 1 we still have that $Tr(M_{1,0})$ is always negative. Based on numerical simulations close to $z=0$ for the whole space case (see Fig. 4 I) we observe that the constant steady state is unstable again only for sufficiently small $z_{1,0}$ and iff $\Delta(M_{1,0})<0$. This leads to two eigenvalues $\lambda_{1}$, $\lambda_{2}$ of different signs, that can be computed using the formulas (34). The value of parameter $s$ that corresponds to the phase transition can be computed by equating (42) to [math] and therefore

[TABLE]

3.3 On phase separated initial conditions

The authors of [53] report segregation experiments between two types of cells. In particular, in some experiments (see Fig.1 i)-k) in [53]) the initial condition is perfectly segregated: at $t=0$ there is a sharp interface between a subdomain with only cells of type $A$ and a subdomain with only cells of type $B$. Depending on various parameters, the experiments show different possible evolutions of the interface (see Fig.1 i)-k) in [53]):

• •

the interface remains sharp and does not move;

• •

the interface remains sharp, does not deform, and moves in one direction;

• •

the interface becomes blurred: the two species start to mix.

Our goal here is to try to relate these different outcomes to different parameter ranges in our model. The heuristic remarks of this paragraph will be backed by numerical simulations in Section 4. We consider here a bounded domain $\Omega\subset\mathbb{R}^{2}$, and make some further hypotheses in order to simplify the system:

• i)

the inter-species interaction is symmetric $\Phi^{AB}=\Phi^{BA}$;

• ii)

the interaction radius $R$ is much smaller than the scale of the experiments;

• iii)

the diffusion is neglected.

Under hypothesis i), the macroscopic system of equation admits a free energy, or Lyapunov functional

[TABLE]

The macroscopic dynamics is a gradient flow of the functional $\mathcal{F}$, with respect to the Wasserstein distance $W_{2}$ [56]. Owing to hypothesis ii), we will replace the potentials $\Phi^{\alpha\beta}$ by Dirac delta functions with weights equal to their integrals, which we call $\gamma^{\alpha\beta}$, respectively. Owing to hypothesis iii), we neglect diffusion, which all together leads to the simplified functional $\mathcal{E}$:

[TABLE]

Under hypotheses i)-iii) we expect that the system evolves in order to minimize (44), with constraints of constant total mass and positivity (we assume that at $t=0$, there is the same amount of $A$ and $B$ cells):

[TABLE]

This problem is easily solved if $\gamma^{AA}<0$ or $\gamma^{BB}<0$. Then the minimum energy is $-\infty$ and it corresponds to one or both species infinitely concentrated. From now on we assume that $\gamma^{AA}>0$ and $\gamma^{BB}>0$, i.e. the intraspecies interactions are repulsive. Another singular case corresponds to the range $\gamma^{AB}<-\sqrt{\gamma^{AA}\gamma^{BB}}$; then there is a pair of constants $(n^{A}>0,n^{B}>0)$ such that

[TABLE]

Now, building a sequence of configurations where $A$ and $B$ cells have concentration $n^{A}/\varepsilon$, $n^{B}/\varepsilon$ on the same domain of measure $\varepsilon$, we see that the associated energy tends to $-\infty$ when $\varepsilon$ tends to [math]. Therefore, we also assume from now on that $\gamma^{AB}>-\sqrt{\gamma^{AA}\gamma^{BB}}$.

Lemma 2

Assume that $\gamma^{AA}>0$, $\gamma^{BB}>0$, and $-\sqrt{\gamma^{AA}\gamma^{BB}}<\gamma^{AB}<\sqrt{\gamma^{AA}\gamma^{BB}}$. Then minimum energy is reached by the homogeneous system with constant densities $f^{A}(x)=f^{B}(x)=1/(2|\Omega|)$.

Proof: Since $\gamma^{AA}\gamma^{BB}-(\gamma^{AB})^{2}>0$, the function

[TABLE]

is convex. Hence, for any functions $f^{A}(x),f^{B}(x)$:

[TABLE]

where $|\Omega|$ is the volume of the domain and $\bar{f}^{S}=\int_{\Omega}f^{S}(x)dx=1/2$, $S=A,B$. This shows immediately that for any function $f^{A}(x),f^{B}(x)$

[TABLE]

The minimum energy is then reached for a homogeneous system, with constant densities, $f^{A}(x)=f^{B}(x)=1/(2|\Omega|)$. $\Box$

We will now assume that $\gamma^{AB}>\sqrt{\gamma^{AA}\gamma^{BB}}$. We want to show that the minimal configuration is perfectly phase separated, ie $\Omega$ divided in two subdomains $\Omega=\Omega^{A}\cup\Omega^{B}$ with $f^{A}=0$ on $\Omega^{B}$ and $f^{B}=0$ on $\Omega^{A}$. Note that the geometry of the subdomains is of no relevance within this simplified model, only their measure matters.

Lemma 3

Among all perfectly phase separated configurations, the following is optimal:

[TABLE]

The associated energy is

[TABLE]

In other words: for any $(f^{A},f^{B})$ perfectly phase separated configuration,

[TABLE]

with equality if and only if (46)-(47) hold.

Proof: First note that the functions

[TABLE]

are convex; hence any perfectly phase separated configuration must be piecewise homogeneous in order to be optimal. Let us now consider a general piecewise homogeneous phase separated configuration:

[TABLE]

with

[TABLE]

Then $n^{A}=1/(2l_{A})$, $n^{B}=1/(2|\Omega|-2l_{A})$, and the configuration is characterized by the parameter $l_{A}$, which we can optimize. The associated energy is

[TABLE]

Minimizing over $l_{A}$, one finds

[TABLE]

with associated energy given by (48). $\Box$

Lemma 4

Assume $\gamma^{AA}>0$ and $\gamma^{BB}>0$, and $\gamma^{AB}>\sqrt{\gamma^{AA}\gamma^{BB}}$. Then the optimal configuration is the perfectly phase separated one described in the previous lemma, ie for any configuration $(f^{A},f^{B})$, we have

[TABLE]

Proof: Take any configuration $f^{A},f^{B}$. Using $\int_{\Omega}f^{A,B}=1/2$ and (48), we have

[TABLE]

where we have used $(\gamma^{AB})>[\gamma^{AA}\gamma^{BB}]^{1/2}$ for the last inequality. Hence the inequality is strict unless $f^{A}f^{B}=0$, ie unless the configuration is phase separated. We have proved that the optimal configuration is phase separated, and given by (47). $\Box$

To summarize, assuming $\gamma^{AA}>0$, $\gamma^{BB}>0$, $\gamma^{AB}>-\sqrt{\gamma^{AA}\gamma^{BB}}$ and putting together Lemmas 2,3,4, one expects the following phenomenology:

• •

If $|\gamma^{AB}|<\gamma_{AA}\gamma_{BB}$, then the homogeneous solution is favoured. The two types of cells should then start to mix and the interface should be blurred.

• •

If $\gamma^{AB}>\sqrt{\gamma_{AA}\gamma_{BB}}$, then the phase separated solution is favoured; in this case the initially sharp interface should stay sharp. There are two subcases:

• –

$\gamma_{AA}\simeq\gamma_{BB}$: then the optimal $l$ is close to $1/2$, which is the initial condition. Hence the interface should not move.

• –

$\gamma_{AA}$ and $\gamma_{BB}$ significantly different: then the optimal $l$ is not close to the initial $1/2$, and one expects the sharp interface to move, as the system tries to approach the energy minimum.

These three scenarios are qualitatively similar to the ones reported in [53] (see their Figure 1), and are seen in numerical simulations, as shown on Fig.6; a quantitative comparison is difficult, since hypothesis i)-iii) are not necessarily satisfied in numerical simulations.

4 Numerical results

The microscopic model described in Section 2.1 is very demanding numerically; it is one reason to introduce the macroscopic model (24), which relies on the double limit $(N_{A},N_{B})\to\infty$, $\varepsilon\to 0$. The averaged microscopic model (22)-(23), obtained in the limit $\varepsilon\to 0$ with $N_{A},N_{B}$ fixed, can be simulated at a reasonable numerical cost, which makes comparisons with the macroscopic model possible.

Numerical simulations for the averaged microscopic model (22)-(23) as well as macroscopic model (24) are performed on a 2D periodic domain $[-7.5,\;7.5]\times[-7.5,\;7.5]$. All simulations are performed with cell detection radii $R_{A}=R_{B}=1$. We fix the diffusion constants $D_{A}=D_{B}=10^{-4}$ and explore different values of the inter- and intra- species repulsion intensities $\kappa^{AA},\kappa^{BB},\kappa^{AB}=s\tilde{\kappa}^{AB},\kappa^{BA}=s\tilde{\kappa}^{BA}$. We explore four different cases:

• •

case 1: For the same intraspecies repulsion $\kappa^{AA}=\kappa^{BB}$ and symmetric inter-species repulsion $\tilde{\kappa}^{AB}=\tilde{\kappa}^{BA}$

• •

case 2: For the same intraspecies repulsion $\kappa^{AA}=\kappa^{BB}$ and non-symmetric inter-species repulsion $\tilde{\kappa}^{AB}<\tilde{\kappa}^{BA}$ ($A$-cells repulse $B$-cells more strongly than the reverse)

• •

case 3: For different intraspecies repulsion $\kappa^{AA}>\kappa^{BB}$ ($A$-cells repulse each other more strongly than $B$-cells), and symmetric inter-species repulsion $\tilde{\kappa}^{AB}=\tilde{\kappa}^{BA}$

• •

case 4: For different intraspecies repulsion $\kappa^{AA}>\kappa^{BB}$ ($A$-cells repulse each other more strongly than $B$-cells) and non-symmetric inter-species repulsion $\tilde{\kappa}^{AB}<\tilde{\kappa}^{BA}$ ($A$-cells repulse $B$-cells more strongly than the reverse)

For each case, we consider two types of initial conditions, (i) when cells are initially randomly distributed, which approaches the homogeneous stationary state, and (ii) when $B$-cells are initially randomly placed on the left-half of the domain and $A$-cells randomly distributed on the right part, which corresponds to the phase separated initial conditions considered in Section 3.3, and we explore two different regimes (stable regime, when $s<s^{*}$ for the interspecies repulsion and unstable regime, when $s>s^{*}$). Table 1 sums up the model parameters used for inter- and intra- species repulsion forces.

In Fig. 5, we show the final states of the simulations for each case and each regime previously described for the microscopic model. $A$-cells are represented as red disks, $B$-cells as green ones and we use $N_{A}=N_{B}=250$ particles for each family of cells. The visualisation of the macroscopic results also uses disks to resemble the microscopic ones, for more details we refer to the Appendix C.

As one can note in the stable regime ($s<s^{*}$) and for initially randomly distributed particles, we observe a homogeneous distribution of particles as expected, with no aggregation. When starting from a front-like intial distribution in the stable regime, $B$- and $A$- cells intermingle at the front. However in the unstable regime (cases $s>s^{*}$), one can observe a segregation of cells by type. In case 1 and 2 (when intraspecies repulsion is the same), we observe the formation of mazes of $B$-cells when starting from an initial homogeneous distribution, and the maintenance of a sharp front when starting from a non-homogeneous initial distribution (Fig. 5 (IB, IIB)). Note that the $B$-cell clusters are smaller in case 2 than in case 1 (i.e when $A$-cells act more strongly on $B$-cells than the reverse compared with the case where inter-species repulsion is symmetric). When starting from a segregated initial condition we observe a slight compression of the green cells by the red ones. This suggests that nonsymmetric interspecies repulsion can favor cell aggregation and domination of a population over the other.

The same observations as in case 2 can be done for case 3 (i.e when intra-species repulsion is stronger in cell $A$ type than in $B$ type and with symmetric interspecies forces, (IIIB)): smaller clusters than in case 1 and slight compression of $B$-family in the case of segregated initial condition. This suggests that inter- and intra- specie repulsion act in the same manner: decreasing the intra-species force for one family has the same impact on the final structures as decreasing the repulsion force of one family onto the other one.

Finally for case 4 (intra-species repulsion is stronger in cell $A$ type than in $B$ type and stronger repulsion of $B$-cells by $A$-cells than the reverse), we observe the emergence of small clusters of $B$-cells in a medium composed of $A$-cells, and a large compression of the $B$-cells with maintenance of sharp borders when starting from a non-homogeneous distribution (IV B).

It is noteworthy that the numerical simulations of the microscopic model in the limit $\epsilon\rightarrow 0$ are in good accordance with the predictions of the macroscopic model. We indeed observe homogeneous or non-homogeneous distributions of particles depending on the model parameters, for values of $s$ in the range predicted by the linear stability analysis.

In Table 2, we show the Volume Fraction (VF) of type B cells (green family) computed on the simulation images of Figs. 5 for the microscopic model and of Fig. 6 for the macroscopic model, starting from homogeneous Initial Conditions (IC) or front-like Initial Conditions. Given a numerical image such as depicted in Fig. 5, the volume fraction corresponds to the number of green pixels over the total amount of pixels in the image. As one can see in Table 2, the volume fraction at equilibrium does not depend on the type of initial conditions, as suggested by the analysis performed in Section 3.3. As expected, we obtain a volume fraction of $50\%$ when the two cell types have the same inter- and intra- species forces, and the volume fraction occupied by specie B decreases as type $B$ cells’ inter- and/or intra- species forces decrease (maintaining the type $A$ cells inter- and intra- species forces constant), due to the compression exerted by the stronger family (type $A$) on the weaker cells (type $B$). These results are in accordance with the theoretical predictions of the macroscopic model, showing that the microscopic and macroscopic model have the same properties. The next section is devoted to deeper numerical comparisons between the two models.

In Figs. 7 (I), we show simulations of the microscopic and macroscopic models for $\kappa^{AA}=4,\kappa^{BB}=\tilde{\kappa}^{AB}=\tilde{\kappa}^{BA}=1$ for which $s^{*}\approx 2.1$ can be computed using the formula (43). Simulations of the microscopic model are performed with $N_{A}=N_{B}=500$ (IA) and $N_{A}=N_{B}=2000$ (IB) particles. Simulations of the macroscopic model correspond to (IC). We consider 7 values of the interspecies repulsion intensity, from left to right: for $2.05=s<s^{*}$, for $s^{*}<s=\{2.15,2.2,2.5,4,6,10\}$. (II) In (IIA-C), we show the values of the quantifiers at time equilibrium as functions of $s$. Fig. (IA) shows the mean elongation of the green clusters, (IIB) shows the number of green clusters and (IIC) shows the Overlapping amount $Q$ described by (74). Black curves are obtained with the microscopic model for $N_{A}=N_{B}=500$ (corresponding to Figures (IA)), yellow curves are for $N_{A}=N_{B}=2000$ and correspond to Figures (IB) and red curves are obtained with the macroscopic model (Figures (IC)).

As one can observe in Fig. 7 (IA-C), we obtain a very good agreement between the microscopic model and the macroscopic simulations. Before the transition from mixed to segregated states (for $s<s^{*}$, first column), the system at equilibrium corresponds to a perfectly mixed state both for the micro- and for the macro- models, while right after the transition (for $s=2.15$), segregation of the two families is observed for both models. As $s$ increases, clusters get more numerous, smaller and rounder. This qualitative observation is supported by the values of the quantifiers plotted in Figs. (IIA-C). For the different values of $s$, we obtain a very good quantitative agreement between the two models, even more so when the number of particles of the micro- model is increased from $N_{A}=N_{B}=500$ to $N_{A}=N_{B}=2000$ (compare black and yellow curves to red ones in Figs. (IIA-C)). These results tend to show that the macroscopic model is a good approximation of the microscopic dynamics as the number of individuals becomes large. However, it is noteworthy that some discrepancy is observed for very large values of $s$. Indeed for $s=10$ (last column of (IA-C)), one can note that the clusters obtained by the macro- model are significantly more elongated than those obtained with the micro- model (compare black/yellow curves to the red one in Fig.7 (IIA)). For the details of image processing used to prepare the figures we refer to the Appendix C.2.

In order to document the discrepancies between the micro- and macro- structures in the case $s=10$, we plot in Fig. 8 the values of the quantifiers computed on the simulation images as functions of time for $\kappa^{AA}=4,\kappa^{BB}=\tilde{\kappa}^{AB}=\tilde{\kappa}^{BA}=1$ for the microscopic model with $N_{A}=N_{B}=500$ (green curves), $N_{A}=N_{B}=2000$ (blue curves), $N_{A}=N_{B}=4000$ (yellow curves) and for the macro model (red curves). Fig. 8 (I) shows the evolution of green cluster elongation, (II) gives the number of cell clusters and (III) shows the overlapping amount $Q$ as function of the simulation time.

Figs. 8 reveal that the structures observed with the microscopic model undergo two transitions as function of the time, enabling us to conclude that the dynamics of the microscopic model exhibits two phases:

• •

Very quickly after initialisation ($t\in[0,100]$), the system segregates the two cell types and reaches an equilibrium value for the number of clusters and the overlapping amount $Q$ (see Figs. 8 (II) and (III). At the time of the segregation, numerous and well-separated elongated clusters are created and then maintained for a long time.

• •

A second transition occurs later in time (around $t\approx 3.10^{3}$), where the clusters change shape to attain a new equilibrium composed of rounder clusters (see the drop in the value of the elongation in Fig. 8 (I)). In this second phase, the number of clusters and their border properties are maintained, but the shape gradually changes to produce very round clusters.

It is noteworthy that this two-phase process is not observed with the macroscopic dynamics (see the red curves of Figs. 8 (I-III). On the contrary for the macroscopic model, the segregation between the two families and the production of clusters appear later than with the microscopic model (around $t\in[100,300]$), and the shape of the clusters (elongated) seems to be at equilibrium. These results tend to show that the macroscopic model fails to capture the second time phase (reorganisation of the clusters) exhibited by the microscopic model.

The good qualitative and quantitative agreement between the micro- and the macro- models has been also confirmed for another set of parameters $\kappa^{AA}=4,\kappa^{BB}=3,\tilde{\kappa}^{AB}=\tilde{\kappa}^{BA}=1$ for which the critical value of $s$ is $s^{*}\approx 3.52$. We do not include the corresponding figures here, as they are very similar to Figs. 8.

5 Conclusions and perspectives

Along the recent biological studies [52] our paper demonstrates that contact cell repulsion on its own can generate pattern formation and cell-sorting in tissues composed of different categories of cells. The present paper provides evidence of this by means of both a microscopic and a consistently derived macroscopic model. The advantage of the macroscopic approach is that it provides a mathematical way to investigate the stability of the equilibria and consequently quantitative criteria for the appearance of these patterns. The validity of this analysis is assessed by numerical comparison between microscopic and macroscopic models. We also show that the model is able to capture the border sharpening observed in the biological study [52]. In the future better quantitative comparison with experiments will allow for more systematic choice of model parameters.

Clearly the macroscopic model captures the behaviour of the mean, and it would be interesting to investigate what fluctuations around these means are induced by the finiteness of the number of cells. Our analysis provides only local information around the threshold for instability, further analysis as in [7] could inform on the type of bifurcation involved, for instance if it is supercritical or subcritical. The consistent derivation of the macroscopic model from the microscopic dynamics is still only formal and it would be desirable to have a rigorous mathematical proof of convergence when the number of particles goes to infinity. Finally, the biological context could be enriched and applied to clinically relevant situations such as cancer for which the type of bifurcation involved could be critically important.

Acknowledgements

The authors wish to thank Anais Khuong for introducing to the biological aspects of the paper. PD acknowledges support from the Engineering and Physical Sciences Research Council (EPSRC) under grant ref. EP/M006883/1 and EP/N014529/1, from the Royal Society and the Wolfson foundation through a Royal Society Wolfson Research Merit Award no. WM130048 and from the National Science Foundation (NSF) under grant RNMS11-07444 (KI-Net). PD is on leave from CNRS, Institut de Mathématiques, Toulouse, France. EZ was supported by the Polish Government MNiSW research grant 2016-2019 ”Iuventus Plus” No. 0888/IP3/2016/74.

Data statement

No new data were generated in the course of this research

Appendix A Derivation of a kinetic model from the microscopic model

In this section, the derivation of a kinetic model from the Individual Based Model of Section 2.1 is performed, following the Approach I. Using the expressions for the individual particle distribution, link distribution and two-particle distribution function defined by Eqs. (8),(9),(10), (12)-(14) and in the limit of a large number of individuals, we have the following formal derivation:

Proposition 5

Assume that in the limit

[TABLE]

the following convergences hold:

[TABLE]

[TABLE]

If we assume the scalings (18) for the rates, and if we assume additional assumptions ((58),(59) and (66),(67) below), then $f^{A},f^{B}$ formally solve:

[TABLE]

where:

[TABLE]

and $g^{AA},g^{BB}$ and $g^{AB}$ formally solve:

[TABLE]

[TABLE]

[TABLE]

*Proof. *

A.1 Evolution equation for the individual particles

For all observable functions $\phi(x)$, we define:

[TABLE]

Similarly, for all two-particle observable functions $\psi(x_{1},x_{2})$, we define:

[TABLE]

where integrals over $x$ are carried over ${\mathbb{R}}^{2}$. Then:

[TABLE]

For the sake of simplicity, the computations for $f^{A}_{N}$ only are developed here. Using (3) and Itô’s formula, we obtain formally:

[TABLE]

As the $dB_{i}^{A}$ are independent and independent of $\nabla_{x}\phi(X^{A}_{i})$, it can be shown that in the limit of a large number of particles this term can be neglected [7]. We get:

[TABLE]

Now, exchanging the sums in $i$ and $k_{1}$ and $i$ and $k_{3}$ in the previous equation, one obtains:

[TABLE]

From the symmetry of $\Phi^{AA}$, we have:

[TABLE]

leading to:

[TABLE]

or again:

[TABLE]

where we have formally integrated by parts to obtain the second equality. We then exchange the order of integration and pass to the limit $N_{A},N_{B}\rightarrow\infty$. If $f^{A}_{N}\rightarrow f^{A}$, $g^{AA}_{N}\rightarrow g^{AA},g^{AB}_{N}\rightarrow g^{AB}$, then:

[TABLE]

where:

[TABLE]

Similarly, we can show that if $f^{B}_{N}\rightarrow f^{B}$, $g^{BB}_{N}\rightarrow g^{BB},g_{N}^{BA}\rightarrow g^{BA}$ as $N_{A},N_{B}\rightarrow\infty$ with $N_{A}/N_{B}\to r_{AB}$, we get:

[TABLE]

with:

[TABLE]

We now turn towards the computation of the inter- and intra- species link distributions.

A.2 Evolution equation for the particle links

Here we develop only the computations for one intraspecies link distribution, namely $g^{AA}_{N}$, the computation of $g^{BB}_{N}$ being similar. From its asymmetry, the computation of the interspecies link distribution needs special treatment and we will develop the computation of $g^{AB}_{N}$ further. We remark that the noise in (3) transforms directly into a linear diffusion term for $f^{A}$, all other contributions, analogously as in 53, vanish in the large $N_{A}$ limit. It is not difficult to see that the same simplification takes place for $g^{AA}_{N}$ in the $N_{A}\rightarrow\infty$ limit. Thus, to reduce the computations we will first use (3) without noise, and reintroduce the diffusion term in the end.

Intraspecies link distribution

Following the same principle as for $f^{A}_{N}$, one can write:

[TABLE]

where $E_{k}$ corresponds to the $k$-th line of (56). For the sake of simplicity, the computation of $E_{1}$ only is developed here. The computation of the other ones are similar and omitted.

From Eq. (3), one obtains:

[TABLE]

Now, exchanging the sums in $k$ and $k_{1}$ and $k$ and $k_{3}$ and using the symmetry of $\Phi^{AA}$, one obtains:

[TABLE]

In the first two lines of the above expression, $k=k_{1}$ in the internal sums play a special role: indeed we know that the $k=k_{1}$ term, at variance with the other terms, always contribute to the sum. On the contrary, $k=k_{3}$ in the last line above should not be distinguished. So, we get:

[TABLE]

Because there is no restriction on the number of links per particle, the sums over $k$ cannot be simplified in this case. In order to express the terms in the $k\neq k_{1}$ sum, we define the number of intra-species links connected to a particle:

[TABLE]

where Card denotes the cardinal of a set. Then, as $N_{A}\rightarrow\infty$, we assume that the following mean-field approximation holds for any chosen link $k_{1}$:

[TABLE]

where

[TABLE]

is the conditional probability of finding an intraspecies link conditioned on the fact that one of the particles of this link has the same location as $i(k_{1})$. Then, as $N_{A}\rightarrow\infty$ , $C_{i}^{k_{1}}$ is the mean number of intraspecies links per particle. The mean number of intraspecies links AA in the volume $dX^{A}_{i(k_{1})}$ is $N_{A}\int g^{AA}(X^{A}_{i(k_{1})},x_{2})dx_{2}$ and the mean number of particles of type A in $dX^{A}_{i(k_{1})}$ is $N_{A}f^{A}(X^{A}_{i(k_{1})})$. Thus:

[TABLE]

So, we get:

[TABLE]

Inserting these expressions in Eq. (57), one obtains, when $N_{A},N_{B}$ are large:

[TABLE]

where

[TABLE]

Finally, we find:

[TABLE]

After the same treatment for $E_{2}$ of Eq. (56) and in the limit $N_{A},N_{B}\rightarrow\infty,\frac{N_{A}}{N_{B}}\rightarrow r_{AB}>0$, one obtains the final equation for $g^{AA}$:

[TABLE]

where,

[TABLE]

Integrating by parts, changing the variables and order of integrals we easily obtain:

[TABLE]

Finally, restoring the noise, we obtain the final equation for $g^{AA}$:

[TABLE]

where

[TABLE]

Eq. (63) does not take into account the phenomena of creation and destruction of intraspecies links. According to the description at the beginning of this paper, our model describes a process of creation of links with rate $\nu^{AA}_{c,N,\varepsilon}$, provided the two type $A$ particles are sufficiently close to each others. Hence, the number of new intraspecies links will be proportional to the number of pairs of particles such that one of them is close to $x_{1}$ and the other one is close to $x_{2}$, whose distance is less than $R$:

[TABLE]

where $h^{AA}(x_{1},x_{2},t)=\lim_{N_{A}\rightarrow\infty}h^{AA}_{N}(x_{1},x_{2},t)$ and where $h^{AA}_{N}(x_{1},x_{2},t)$ is the two particle distribution function for particles of type A defined by Eq.(12). This number has to be decreased by the number of pairs of particles of the same type that are already connected by existing links

[TABLE]

Therefore, the average number of new intraspecies links created in the interval $[t,t+\Delta t[$ is equal to

[TABLE]

Dividing this expression by $N_{A}$ and using (18), the rate of creation of new intraspecies link at $(x_{1},x_{2})$ is in the limit $N_{A}\to\infty$:

[TABLE]

Notice the scaling of $\nu^{AA}_{c,N,\varepsilon}$ in (18): it ensures that among the $O(N_{A}^{2})$ possible links, only $O(N_{A})$ are effectively present, and the total number of AA links $K^{AA}=O(N_{A})$. The rate of destruction of existing intraspecies link at $(x_{1},x_{2})$ is:

[TABLE]

Including these source terms in (63), we obtain

[TABLE]

Quite straightforwardly, we can show that in the limit $N_{A},N_{B}\rightarrow\infty,\frac{N_{A}}{N_{B}}\rightarrow r_{AB}>0$, $g^{BB}(x_{1},x_{2},t)$ solves:

[TABLE]

where we have assumed the scaling (18) for $\nu^{BB}_{c,N,\varepsilon}N_{B}$, and

[TABLE]

We stress the fact that $\Phi^{AB}(x_{1},x_{2})$ (force of a particle of type $B$ close to $x_{2}$ exerted on a particle $A$ close to $x_{1}$) is not necessarily equal to $\Phi^{BA}(x_{2},x_{1})$ (force of a particle of type $A$ close to $x_{1}$ exerted on a particle $B$ close to $x_{2}$).

Computations of the interspecies link distribution

Here, we develop the computations for the interspecies link distribution $g^{AB}(x_{1},x_{2},t)$. Proceeding as before, we write:

[TABLE]

where, ignoring the noise and using Eq.(3):

[TABLE]

Exchanging the sums in the first two terms and using the symmetry of $\Phi^{AA}$, we have

[TABLE]

Now, in the same spirit as before we define the number of interspecies links linked to a particle of type $A$:

[TABLE]

and we make the following mean-field assumption

[TABLE]

where $P^{AB}(X^{A}_{i(k_{3})},x_{2})$ is the conditional probability of finding an interspecies link conditioned on the fact that the type $A$ particle of the link has the same location as $i(k_{3})$:

[TABLE]

Now as $N_{A},N_{B}\rightarrow\infty$, $C_{i,A}^{k_{3}}$ is the mean number of interspecies links per particle of type A. The mean number of interspecies links the type A particle of which belonging to the volume $dX^{A}_{i(k_{3})}$ is $N_{A}\int g^{AB}(X^{A}_{i(k_{3})},x_{2})dx_{2}$, and the mean number of particles of type A is $N_{A}f^{A}(X^{A}_{i(k_{3})})$. Therefore,

[TABLE]

leading, when $N_{A},N_{B}$ are large, to

[TABLE]

and

[TABLE]

and altogether:

[TABLE]

Now, exchanging the sums in the last term of $E_{1}$, we obtain

[TABLE]

where, as previously, we distinguish the link $k_{3}$, that always contributes to the inner sum. Using the theory previously developed, we can write, when $N_{A},N_{B}$ are large:

[TABLE]

while

[TABLE]

Altogether, performing the same computations for $E_{2}$ in (65), we obtain:

[TABLE]

By carefully performing integration by parts and change of order of integrals, we can obtain:

[TABLE]

Finally, note that in the case of interspecies links, the average number of new links created in the interval $[t,t+\Delta t[$ is equal to

[TABLE]

Assuming for $\nu_{c,N,\varepsilon}^{AB}$ the scaling (18) as the number of particles tends to infinity and adding the noise, we obtain the final equation for $g^{AB}$:

[TABLE]

where $h^{AB}(x_{1},x_{2},t)=\lim_{N_{A},N_{B}\rightarrow\infty}h^{N_{AB}}(x_{1},x_{2},t)$ with $h^{N_{AB}}(x_{1},x_{2},t)$ defined by Eq. (14).

Appendix B Scaling of the kinetic model

B.1 Dimensionless Equations

In order to express the problem in dimensionless variables, we denote by $t_{0}$ the unit of time and $x_{0}$, $f^{S}_{0}=\frac{1}{x_{0}^{2}}$, $g^{ST}_{0}=\frac{1}{x_{0}^{4}}$, $g^{ST}_{0}=\frac{1}{x_{0}^{4}}$ the units of space and distribution functions, where $S$ and $T$ can be either $A$ or $B$ and refer to the particle type. The scaling of $f^{S}(x,\theta)$, $g^{ST}(x_{1},x_{2})$ and $h^{ST}(x_{1},x_{2})$ comes from the fact that they are probability distribution functions on a 2D domain. The following dimensionless variables are defined:

[TABLE]

and the following dimensionless parameters are introduced:

[TABLE]

where we assumed that the potential scales as the potential energy $\frac{x_{0}^{2}}{t_{0}^{2}}$. We first have:

[TABLE]

and:

[TABLE]

In this new set of variables, choosing $S=A,T=A$ (the same scaling apply for the other equations) and omitting the primes and bar for clarity, Eqs. (51)-(52) become :

[TABLE]

and

[TABLE]

Finally, we choose the space and time scales $x_{0}$, $t_{0}$ such that $\mu=1$.

B.2 Scaled equations

In order to describe the system at a macroscopic scale, a small parameter $\varepsilon\ll 1$ is introduced and the space and time units are set to $\tilde{x_{0}}=\varepsilon^{-1/2}x_{0}$, $\tilde{t_{0}}=\varepsilon^{-1}t_{0}$. The variables $x$, $t$, $R$ and unknowns $f$ and $g$ are then correspondingly changed to $\tilde{x}=\sqrt{\varepsilon}x$, $\tilde{t}=\varepsilon t$, $\tilde{R}=\sqrt{\varepsilon}R$. Therefore, $\tilde{f}^{S}(\bar{x})=\varepsilon^{-1}f^{S}(x)$, $\tilde{g}^{ST}(\tilde{x}_{1},\tilde{x}_{2},\tilde{t})=\varepsilon^{-2}g(x_{1},x_{2},t)$ and $\tilde{h}^{ST}(\tilde{x}_{1},\tilde{x}_{2},\tilde{t})=\varepsilon^{-2}h(x_{1},x_{2},t)$. The diffusion constant is supposed to be of order 1, $D^{S}=\tilde{D}^{S}$, and we suppose that the interaction potentials scale as $\Phi^{ST}(x_{1},x_{2})=\tilde{\Phi}^{ST}(\tilde{x}_{1},\tilde{x}_{2})$. Then,

[TABLE]

and with $\mu=1$, we obtain the same equation for $\tilde{f}^{S}$ (for $S$=A for instance):

[TABLE]

In order to simplify the analysis of the system, the process of linking/unlinking is supposed to occur at a very fast time scale: this is the meaning of the $\varepsilon$-scaling of the rates in (18). As $\chi_{|x_{1}-x_{2}|\leq R}=\chi_{|\tilde{x}_{1}-\tilde{x}_{2}|\leq\tilde{R}}$, we have (for $S,T=$(A,A)):

[TABLE]

Now, we aim to pass to the limit $\varepsilon\rightarrow 0$. Denoting $f_{\varepsilon}^{S}=\tilde{f}^{S}$, $g_{\varepsilon}^{ST}=\tilde{g}^{ST}$ and $h_{\varepsilon}^{ST}=\tilde{h}^{ST}$, we want to derive the same system of macroscopic equations as in Proposition 1 of section 2.2.

*Proof. * From Eq.(69) generalized to $S,T$ and using the assumption $h^{ST}_{\varepsilon}(x_{1},x_{2},t)=f_{\varepsilon}^{S}(x_{1},t)f_{\varepsilon}^{T}(x_{2},t)$, and dropping the tildas, we have:

[TABLE]

Therefore in the limit $\varepsilon\rightarrow 0$, we have that

[TABLE]

for all $S,T=A,B$. Plugging the expression of $g^{ST}$ into the equation for $f^{S}$, we have

[TABLE]

where

[TABLE]

Therefore, if the potentials $\Phi^{ST}(x_{1},x_{2})=U^{ST}(|x_{1}-x_{2}|)$, we can write

[TABLE]

for some potentials $\tilde{\Phi}^{ST}$ such that:

[TABLE]

Recall that due to the scaling of linking frequencies with $N_{A},N_{B}$ we obtain in the limit that $\frac{\nu^{BA}_{c}}{\nu^{BA}_{d}}=r_{AB}\frac{\nu^{AB}_{c}}{\nu^{AB}_{d}}$. With this observation the proof of Proposition 1 is complete.

Appendix C Numerical data visualisation

C.1 Macroscopic model

For the numerical simulations of the macroscopic model we consider a periodic square $[-7.5,\;7.5]\times[-7.5,\;7.5]$ discretized with a space step $\Delta x=0.3$ and the time step $\Delta t=5\cdot 10^{-3}$. As in the microscopic case, we fix the radii of all types of interactions to $R=1$ and consider the same diffusion coefficient for the species $D=10^{-4}$. We use the numerical scheme introduced in [17] that was developed in our recent work [8] to study aggregation-diffusion equation for single species.

We consider all four cases of interaction intensities $\kappa^{ST}$ as in the microscopic case, as summarized in Table 1. The initial densities for species $A$ and species $B$ are random perturbations of constant functions, The constants are chosen such that the total mass is equal to 1, i.e.

[TABLE]

where $X_{S}(x,y)$ is uniformly distributed random variable between 0 and 1.

To visualise the macroscopic simulations we plot green or red balls in the regions where the densities of green or red cells, respectively, dominate. The balls are of radius $0.3$ with a center at the center of the corresponding pixel perturbed by uniformly distributed random number from interval $[-0.2,0.2]$ in $x$ and $y$ direction. We determine that the concentration of green cells dominates over the red ones if the difference between their densities is larger than parameter $tres_{1}$. It is equal to $5.0e-4$ times the maximum of the densities of green or red cells at the final time of simulation. If the difference between the densities of cells is less than $tres_{1}$ but more than $tres_{2}=tres_{1}-5.0e-6$ we plot randomly red or green ball. The ”empty” black regions on the figures (see for example Figure 6 and last row of Figure 7) correspond to the case when the difference between the densities of the cells is less than $tres_{2}$.

C.2 Image processing

In order to compare quantitatively the microscopic and macroscopic simulations, we use image processing tools to define quantifiers of the structures observed in the simulation images. As we will restrict ourselves to regimes where the $B$-family clusterize (represented in green), the developed tools aim to describe type B clusters but can easily be adapted to detect red clusters. Given a RGB image transformed into a binary image via thresholding of the intensity, we use morphological operators to erase the isolated pixels and dilate the image using a binary gradient mask. Interior gaps are then filled and the new image is smoothed via morphological operations. The boundaries of the isolated clusters are then detected using the Moore-Neighbor tracing algorithm modified by Jacob’s stopping criteria, implemented in the intrinsic matlab function bwconncomp that we couple with a function that appropriately converts the output of bwconncomp to take into account periodic boundary conditions. Once the clusters are appropriately separated and borders are detected, we finally use the matlab intrinsic function regionprops to measure each 8-connected object (region) of the image. This image processing enables us to compute the number $N_{g}$ of green clusters in the image. The elongation of each cluster is given by the eccentricity of the ellipse that has the same second-moments as the region (cluster). Finally, we define the overall overlapping amount $Q$ as:

[TABLE]

where $P^{i}_{g}$ (resp. $P^{i}_{r}$) is equal to 1 if pixel $i$ has a non-zero green component (resp. red) and $N_{P}$ is the total number of pixels in the image. Thus defined, $Q\approx 1$ when the two families are perfectly mixed (corresponding to all pixels having both red and green components), and $Q=0$ correspond to completely separated phases (where each pixel is either green or red). Therefore, parameter $Q$ enables to distinguish between homogeneous and segregated states. However, note that $Q$ does not give any information on the form of the clusters.