AlgorithmexpList of Algorithms

# Modelling adhesion in stochastic and mean-field models of cell migration

Shahzeb Raja Noureen<EMAIL_ADDRESS>Centre for Mathematical Biology,
University of Bath, Claverton Down, Bath, BA2 7AY, UK.  Richard L. Mort
Division of Biomedical and Life Sciences, Faculty of Health and Medicine,
Furness Building, Lancaster University, Bailrigg, Lancaster, LA1 4YG, UK.
Christian A. Yates  Centre for Mathematical Biology, University of Bath,
Claverton Down, Bath, BA2 7AY, UK.

(28th August 2024)

###### Abstract

Adhesion between cells plays an important role in many biological processes
such as tissue morphogenesis and homeostasis, wound healing and cancer cell
metastasis. From a mathematical perspective, adhesion between multiple cell
types has been previously analysed using discrete and continuum models
including the Cellular Potts models and partial differential equations (PDEs).
While these models can represent certain biological situations well, Cellular
Potts models can be computationally expensive and continuum models only
capture the macroscopic behaviour of a population of cells, ignoring
stochasticity and the discrete nature of cell dynamics. Cellular automaton
models allow us to address these problems and can be used for a wide variety
of biological systems. In this paper, we consider a cellular automaton
approach and develop an on-lattice agent-based model (ABM) for cell migration
and adhesion in a population composed of two cell types. By deriving and
comparing the corresponding PDEs to the ABM, we demonstrate that cell
aggregation and cell sorting are not possible in the PDE model. Therefore, we
propose a set of stochastic mean equations (SMEs) which better capture the
behaviour of the ABM in one and two dimensions.

Keywords: cell-cell adhesion, lattice-based model, agent-based model, cell
migration, pattern formation, continuum model.

##  I Introduction

Collective migration of cells is an important biological process responsible
for the correct formation of tissues during development and their maintenance
throughout life. Although each cell is an individual unit, its movement cannot
always be considered independent of other cells in its environment. Cells have
been observed to exhibit a wide range of interactions with other cells in
their neighbourhood. Some of the well-known interactions include volume
exclusion (cells prevent other cells from occupying the same space) (Loeb,
1921; Abercrombie and Heaysman, 1954; Abercrombie, 1979, 1997), contact
inhibition of locomotion (when two cells collide they change directions or
stop moving) (Lee  _et al._ , 1995; Carmona-Fontaine  _et al._ , 2008; Deutsch
_et al._ , 2021), pushing (a moving cell can push cells in front) (Yates  _et
al._ , 2015), pulling (a moving cell can pull neighbouring cells with it)
(Yamanaka and Kondo, 2014; Davis and Trinkaus, 1981; Chappelle and Yates,
2019), and cell-cell swapping (two neighbouring cells can swap positions with
each other) (Raja Noureen  _et al._ , 2023; Osborne  _et al._ , 2017). Another
interaction to attract significant biological and mathematical interest is
cell-cell adhesion (Foty and Steinberg, 2005; Schmidt and Friedl, 2010) in
which cells tend to ‘stick’ to one another through surface binding proteins
known as cell adhesion molecules (CAMs) (Armstrong  _et al._ , 2006; Foty and
Steinberg, 2005). Cell-cell adhesion is implicated in a wide range of
important biological processes including the formation of tissues during
development and their stability post-development; the invasion and metastasis
of tumour cells; and developmental processes such as gastrulation and
vasculogenesis (Armstrong  _et al._ , 2006; Foty and Steinberg, 2005).

Several mathematical models have been successfully employed to represent cell-
cell adhesion. Discrete lattice-based stochastic models and deterministic
continuum models based on partial differential equations (PDEs) are both
commonly used paradigms. Armstrong  _et al._ (2006) propose a continuum
description with non-local adhesive interactions to model differential
adhesion, representing the average behaviour of two well-mixed populations of
cells. Differential adhesion describes the adhesion of two or more populations
of cells each of which has different adhesion properties. The differential
adhesion hypothesis (DAH) views the two populations of cells as immiscible
fluids with different surface tensions. It states that when two populations of
cells (such as embryonic cells) with different adhesion properties are mixed,
the energetically favourable configuration of the system in the steady state
depends on the relative adhesion strengths of the cells (Foty and Steinberg,
2005, 2004; Steinberg, 1963; Armstrong  _et al._ , 2006). Armstrong  _et al._
(2006) use their model to replicate experimentally observed cell sorting
patterns in single-species and two-species scenarios as predicted by the DAH
(Foty and Steinberg, 2005, 2004; Steinberg, 1963; Schötz  _et al._ , 2008;
Armstrong  _et al._ , 2006). One limitation of this continuum model is that it
relies upon characterising the mathematical functions governing the adhesive
forces acting on cells, which can be difficult to discern.

Anguige and Schmeiser (2008) develop a continuum model with local adhesive
interactions. Their approach considers cells on a one-dimensional discrete
lattice in which each lattice site can be occupied by multiple cells at the
same time and cells adhere to one another with strength, $\alpha>0$. They find
that for $\alpha<0.75$ cells tend to diffuse out towards a uniform steady
state as time progresses and for $\alpha>0.75$ aggregations of cells in space
can be observed.

Graner and Glazier (1992) and Glazier and Graner (1993) use an extended
Cellular Potts model – a discrete approach in which each cell is represented
by multiple lattice sites, allowing the cell to change shape according to a
given energy minimisation function (Graner and Glazier, 1992; Glazier and
Graner, 1993; Osborne  _et al._ , 2017) – to replicate cell sporting patterns
in biology from a uniformly mixed initial configuration of two populations.
Cellular Potts models can incorporate mechanical features such as membrane
tension and cell-cell and cell-substrate adhesion which, in conjunction with
the shape-changing property of the cells, adds an extra layer of realism to
these models (Osborne  _et al._ , 2017). For large cell populations,
simulating these models can be computationally expensive.

Simpson  _et al._ (2010a) use a cellular automaton model to represent the
adhesion of breast cancer cells. We will review this model in detail here
since the two-species adhesion model we propose later in this paper is based
on the single-species model of Simpson  _et al._ (2010a).

Modelling a monolayer of breast cancer cells, the authors consider agents on a
regular two-dimensional lattice with $L_{x}$ compartments in the horizontal
direction and $L_{y}$ compartments in the vertical direction. Each lattice
site is square with length $\Delta$ such that a site can only contain one
agent at a time making the model a volume exclusion process. An agent
occupying site $(i,j)$ can attempt to move to a neighbouring site (up, down,
left or right) with movement rate $m$ such that $m\delta t$ is the probability
of the agent attempting a move during the infinitesimally small time window of
duration $\delta t$. If the neighbouring site chosen to move into is empty,
the agent successfully moves and if the neighbouring site is already occupied
by another agent the move is aborted due to volume exclusion. Movement is
further restricted by adhesion to neighbours in the Von-Neumann neighbourhood.
The strength of the adhesive bond that two neighbouring agents share is
denoted as $q\in[0,1]$ with $q=0$ meaning no adhesion and $q=1$ meaning
maximal adhesion (i.e. no movement is possible). Consequently, the probability
of movement depends on the agent breaking current bonds with its neighbours.
Once selected to move into an empty site, the probability of overcoming the
adhesive forces to complete the move is $(1-q)^{\Sigma_{z\in Z_{ij}}C_{z}}$
where $Z_{ij}=\\{(i-1,j),(i+1,j),(i,j-1),(i,j+1)\\}$ is the set of the indices
corresponding to the four neighbouring sites of the site $(i,j)$ and $C_{z}$
is a binary indicator whereby $C_{z}=1$ if the neighbouring site $z\in Z_{ij}$
is occupied or $C_{z}=0$ otherwise.

To derive the continuum model, let $\hat{C}_{ij}(t)$ be the continuous
approximation of the density at the site $(i,j)$ at time $t$. We can write
down the master equation for the $\hat{C}_{ij}(t+\delta t)$ in terms of the
density of the neighbouring sites at time $t$,

$\displaystyle\hat{C}_{ij}(t+\delta t)$
$\displaystyle=\hat{C}_{ij}(t)+\frac{m}{4}\delta
t\hat{C}_{i-1,j}(t)(1-\hat{C}_{ij}(t))(1-q)^{\Sigma_{z\in
Z_{i-1,j}}\hat{C}_{z}}$ $\displaystyle\quad+\frac{m}{4}\delta
t\hat{C}_{i+1,j}(t)(1-\hat{C}_{ij}(t))(1-q)^{\Sigma_{z\in
Z_{i+1,j}}\hat{C}_{z}}$ $\displaystyle\quad+\frac{m}{4}\delta
t\hat{C}_{i,j-1}(t)(1-\hat{C}_{ij}(t))(1-q)^{\Sigma_{z\in
Z_{i,j-1}}\hat{C}_{z}}$ $\displaystyle\quad+\frac{m}{4}\delta
t\hat{C}_{i,j+1}(t)(1-\hat{C}_{ij}(t))(1-q)^{\Sigma_{z\in
Z_{i,j+1}}\hat{C}_{z}}$ $\displaystyle\quad-\frac{m}{4}\delta
t\hat{C}_{ij}(t)(1-\hat{C}_{i-1,j}(t))(1-q)^{\Sigma_{z\in
Z_{i,j}}\hat{C}_{z}}$ $\displaystyle\quad-\frac{m}{4}\delta
t\hat{C}_{ij}(t)(1-\hat{C}_{i+1,j}(t))(1-q)^{\Sigma_{z\in
Z_{i,j}}\hat{C}_{z}}$ $\displaystyle\quad-\frac{m}{4}\delta
t\hat{C}_{ij}(t)(1-\hat{C}_{i,j-1}(t))(1-q)^{\Sigma_{z\in
Z_{i,j}}\hat{C}_{z}}$ $\displaystyle\quad-\frac{m}{4}\delta
t\hat{C}_{ij}(t)(1-\hat{C}_{i,j+1}(t))(1-q)^{\Sigma_{z\in
Z_{i,j}}\hat{C}_{z}}$ (1)

The RHS of Equation (I) corresponds to the possible ways in which the
occupancy of the site $(i,j)$ can change. Either:

  1. 1.

nothing happens and the occupancy at time $t+\delta t$ is equal to the
occupancy at time $t$ (first term on the RHS), or;

  2. 2.

the possibility that the site $(i,j)$ is empty and an agent from a
neighbouring site moves into the site $(i,j)$ breaking its existing bonds with
neighbouring cells (if any) in doing so (second term on the RHS of line 1 and
the terms in lines 2 to 4), or;

  3. 3.

the possibility that the site $(i,j)$ is occupied at time $t$ and the
occupying agent moves out of site $(i,j)$ to one of the empty neighbouring
sites, breaking its existing bonds with neighbouring agents in doing so (lines
5 to 8 on the RHS).

Taylor expanding the terms in Equation (I) around site $(i,j)$ keeping terms
up to second order and taking the limit as $\Delta,\delta t\to 0$ gives us the
PDE for the mean continuum density $C(t,x,y)$,

${}\frac{\partial C}{\partial t}(t,x,y)=D_{0}\nabla\cdot[D(C(t,x,y))\nabla
C(t,x,y)],$ (2)

where,

$D_{0}=\lim_{\Delta,\delta t\to 0}\frac{m\Delta^{2}}{4\delta t},$ (3)

is the diffusion constant given that $\Delta^{2}/\delta t$ is held constant in
the diffusive limit and,

$D(C)=(1-q)^{4C}[1-4C(C-1)\ln(1-q)],$ (4)

is the density-dependent diffusion coefficient.

Figure 1: Comparison between the column-averaged ABM (black) and the
vertically-averaged solution of PDE (2) (red) for $q=0,0.25,0.5,0.9$ at times
$t=$ 0, 100, 1,000 with zero-flux boundary conditions. The movement rate $m=1$
in all cases. The domain is initialised with $L_{x}=400$ and $L_{y}=20$.
Agents are initialised on the domain such that all the sites between
$181\leqslant x\leqslant 220$ are occupied by an agent. The agents evolve
according to the ABM with zero-flux boundary conditions and column-averaged
density of agents at $t=$ 0, 100 and 1,000 are shown in black. The black
arrows indicate the direction of increasing time. Note: in the $q=0.9$ case
the PDE became unstable hence the numerical solution is absent from the panel
1.

We compare the column-averaged density of cells in the ABM,

$\bar{C}_{i}(t)=\frac{1}{L_{y}}\sum_{j=1}^{L_{y}}\langle C_{ij}(t)\rangle,$
(5)

where $\langle C_{ij}(t)\rangle$ is the average over $K$ repeats of the ABM,
to the vertically-averaged solution of Equation (2) given by,

$\bar{C}(t,x)=\frac{1}{L_{y}}\int_{y=0}^{y=L_{y}}C(t,x,y).$ (6)

We initialise the domain with $L_{x}=400$ and $L_{y}=20$ with zero-flux
boundaries such that the middle 10% of the lattice, corresponding to
$181\leqslant i\leqslant 220$ is occupied for all $j$. We let the positions of
the agents evolve according to the ABM with $m=1$ for $q\in[0,0.25,0.5,0.9]$.
Equation (2) is solved numerically with the same parameters under the same
initial and boundary conditions as the ABM and the solution is averaged over
the vertical extent of the domain (as in Equation (6)). The results are shown
in Figure 1.

We note that the agents move less as we increase the adhesion strength. This
agrees with our understanding of the ABM: the higher the adhesion the less
likely agents are to break free of their neighbours and hence the less likely
agents are to move. For $q=0$ the agreement between the ABM and PDE is
excellent, as expected, since in the absence of adhesion, the movement process
becomes a simple volume exclusion process (Liggett, 1999; Simpson  _et al._ ,
2007, 2009a, 2009b, 2010b, 2010a; Cai  _et al._ , 2007). As we increase the
adhesion strength to $q=0.25$ and $q=0.5$ slight deviations can be seen
between the ABM and the PDE for large times. Nevertheless in general, the
agreement remains promising. This disparity in agreement between the ABM and
the PDE for higher adhesion strengths is due to the spatial correlations
between lattice site occupancies in the ABM which are unaccounted for by the
PDE approximation of the ABM.

Figure 2: Colourmap showing the dependence of the PDE diffusion coefficient
$D(C)$ (see Equation (4)), on the adhesion strength $q$ and the species
density $C$. The dashed curve is $q^{\star}$ (see Equation (7)) which bounds
the region in $(C,q)$ space for which $D(C)<0$. The smallest value of $q$ for
which $D(C)$ can be 0 is $q=1-\mathrm{e}^{-1}$ which permits $D(C)=0$ for
$C=0.5$.

Increasing the adhesion strength further to $q=0.9$ the PDE becomes unstable
as the macroscopic diffusion coefficient $D(C)$ defined in Equation (4)
becomes negative for some regions of our space. The ABM shows minimal
dispersion of agents which makes the density profiles at $t=0$, 100 and 1000
appear nearly indistinguishable from each other. The region in the $(C,q)$
space for which $D(C)<0$ is illustrated as a heatmap in Figure 2. The black
dashed curve is defined by,

$q^{\star}=1-\exp\left\\{\frac{1}{4C(C-1)}\right\\},$ (7)

and bounds the region in $(C,q)$ space for which $D(C)<0$. A positive
diffusion coefficient implies classical diffusion (in which mass descends
gradients of mass) and a negative value implies aggregation (as mass climbs
gradients of mass). The lowest value of adhesion strength which admits a non-
positive diffusion coefficient is $q=1-\mathrm{e}^{-1}$. As a result, the PDE
becomes unstable which suggests that adhesion-aided cell sorting into
aggregation patterns such as engulfment cannot be observed using the PDEs.

In the next section, we develop a one-dimensional adhesion model for two
species with cell-cell swapping. We derive the continuum model and show that
it does a poor job of approximating the mean behaviour of the ABM when the
adhesion between agents is strong. We then propose an alternative way to model
the mean behaviour of the two species ABM which not only better represents the
mean behaviour of the ABM than the continuum model but also allows us to
replicate patterns observed in developmental biology.

##  II Two-species model

In this section, we extend the single-species adhesion model to a population
of two species, A and B. We model cell movement on a one-dimensional lattice,
derive the corresponding one-dimensional PDEs and compare the average
behaviour of the ABM to the numerical solution of the PDEs. We also explore
the possibility of self-organisation into patterns as a result of cell
sorting.

We initialise the lattice with $L$ sites of width $\Delta$. Let $m_{A}$ and
$m_{B}$ be the rates of movement of species A and B, respectively, such that
$m_{A}\delta t$ and $m_{B}\delta t$ are the probabilities that a type-A or
type-B agent moves to either of its two neighbouring sites during a small time
window of duration $\delta t$. If the target site is vacant the agent
successfully moves. If, however, the target site is already occupied by
another agent we attempt to swap the positions of the two agents with
probability $\rho$ (Raja Noureen  _et al._ , 2023). A successful swap results
in the two neighbouring agents exchanging positions with each other and a
failed swap results in the event being aborted in which case the two agents do
not move.

To model cell-cell adhesion we introduce the adhesion strengths $p$, $q$ and
$r$ ($0\leqslant p,q,r\leqslant 1$) where $p$ models adhesion between two
type-A agents, $q$ models adhesion between two type-B agents and $r$ models
adhesion between a type-A agent and a type-B agent. Defining
$Z_{i}=\\{i-1,i+1\\}$ as the set of positions of sites in the neighbourhood of
site $i$, the probability that a type-A agent at site $i$ breaks its existing
connections given that it is attempting to move is given by,

$p_{\text{break}}^{\text{agent}}=(1-p)^{\mathrm{\sum}_{z\in
Z_{i}}A_{z}}(1-r)^{\sum_{z\in Z_{i}}B_{z}},$ (8)

where $A_{z}$ is a binary indicator such that $A_{z}=1$ if the neighbour at
position $z\in Z_{i}$ is a type-A agent and 0 otherwise. Similarly, $B_{z}$ is
a binary indicator such that $B_{z}=1$ if the neighbouring agent at position
$z\in Z_{i}$ is a type-B agent and 0 otherwise. Similarly, if the agent at
site $i$ is a type-B agent then,

$p_{\text{break}}^{\text{agent}}=(1-q)^{\mathrm{\sum}_{z\in
Z_{i}}B_{z}}(1-r)^{\sum_{z\in Z_{i}}A_{z}}.$ (9)

If the target site is already occupied by another agent and the agents are
attempting to swap positions then the probability that the target agent breaks
its existing adhesive bonds given the agent at the target site is of type A is
given by,

$p_{\text{break}}^{\text{target}}=(1-p)^{\mathrm{\sum}_{y\in
Y_{z}}A_{y}}(1-r)^{\sum_{y\in Y_{z}}B_{z}},$ (10)

where $Y_{z}$ is the set of sites in the neighbourhood of the target site.
Likewise, if the agent at the target site is a type-B agent,

$p_{\text{break}}^{\text{target}}=(1-q)^{\mathrm{\sum}_{y\in
Y_{z}}B_{y}}(1-r)^{\sum_{y\in Y_{z}}A_{z}}.$ (11)

Altogether the probability of movement, given an agent and a target site have
been chosen, is $p_{move}=p^{\text{agent}}_{\text{break}}(1-c_{y})+\rho
p^{\text{agent}}_{\text{break}}p^{\text{target}}_{\text{break}}c_{y}$ where
$c_{y}\in\\{0,1\\}$ is the occupancy of the target site with 0 meaning the
target site is vacant and 1 meaning that the site is occupied. The
probabilities $p^{\text{agent}}_{\text{break}}$ and
$p^{\text{target}}_{\text{break}}$ depend on whether the focal agent and, in
the case of cell-cell swapping the agent being swapped with, are type-A or
type-B agents, as indicated above in Equations (8)-(11). Note that in the case
of two swapping agents, we are modelling the breaking and reforming of the
adhesive bonds between the two swapping agents explicitly although it would
also be a valid choice to ignore the bonds between the two swapping agents
altogether as they would be reformed following a successful swap.

In Figure 33 we show an example of an agent at site $i$ attempting to move
into an already occupied site $i+1$. The focal agent is of type A with a
type-A agent on its left (site $i-1$) and a type-B agent on its right (site
$i+1$), thus $p^{\text{break}}_{\text{agent}}=(1-p)(1-r)$. The agent at the
target site is a type-B agent with a type-A agent on its left (site $i$) and a
type-B agent on its right (site $i+2$) and therefore
$p^{\text{break}}_{\text{target}}=(1-q)(1-r)$. Hence, the probability of a
successful swap between the pair of agents on this occasion is
$p_{\text{move}}=\rho(1-p)(1-q)(1-r)^{2}$. If the swap is successful then the
two agents make new connections with their neighbours following the move (see
3).

Figure 3: Schematic illustrating an example of the two species ABM with
adhesion. A red site represents a type-A agent, a green site represents a
type-B agent and a black site represents an empty site. Adhesion strengths
between neighbouring pairs of agents are shown as letters $p$, $q$ and $r$
alongside a dash representing a bond between a couple of agents. The initial
configuration of the lattice is shown in 3. The arrow symbolises the red agent
at site $i$ being chosen to attempt to move into the already occupied site
$i+1$. In 3 we show the result of a successful swap. Note that the labels
denoting the adhesive bonds between pairs of agents have been updated in 3.

The agent-based model helps us understand and model interactions between
individual agents. The average behaviour of a large number of cells can be
found by averaging the ABM over a large number of repeats which can be
computationally expensive. PDE-based continuum models provide us with a cheap
and efficient way to approximate the mean behaviour of the discrete ABM. In
the next section, we derive the continuum model corresponding to the one-
dimensional ABM and we compare the two modelling paradigms to assess how well
the continuum model approximates the mean behaviour of the ABM.

###  II.1 One-dimensional continuum model

Let $\hat{A}_{i}(t)$ be the continuous approximation of the mean density of
species A at the site with index $i$ at time $t$ and let $\hat{B}_{i}(t)$ be
the same for species B. Then the continuous approximation of the mean
densities of species A and species B, respectively, at the site $i$ at time
$t+\delta t$, where $\delta t$ is infinitesimally small change in time, are
given by,

$\displaystyle\hat{A}_{i}(t+\delta t)=\hat{A}_{i}+\frac{m_{A}}{2}\delta
t\hat{A}_{i-1}(1-\hat{A}_{i}-\hat{B}_{i})(1-p)^{\Sigma_{z\in
Z_{i-1}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i-1}}\hat{B}_{z}}$
$\displaystyle\quad+\frac{m_{A}}{2}\delta
t\hat{A}_{i+1}(1-\hat{A}_{i}-\hat{B}_{i})(1-p)^{\Sigma_{z\in
Z_{i+1}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i+1}}\hat{B}_{z}}$
$\displaystyle\quad-\Big{(}\frac{m_{A}}{2}\delta
t(1-\hat{A}_{i-1}-\hat{B}_{i-1})(1-p)^{\Sigma_{z\in
Z_{i}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i}}\hat{B}_{z}}$
$\displaystyle\quad+\frac{m_{A}}{2}\delta
t(1-\hat{A}_{i+1}-\hat{B}_{i+1})(1-p)^{\Sigma_{z\in
Z_{i}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i}}\hat{B}_{z}}\Big{)}$
$\displaystyle\quad+\frac{(m_{A}+m_{B})}{2}\delta
t\rho\hat{A}_{i-1}\hat{B}_{i}(1-p)^{\Sigma_{z\in
Z_{i-1}}\hat{A}_{z}}(1-q)^{\Sigma_{y\in Y_{i}}\hat{B}_{y}}(1-r)^{\Sigma_{z\in
Z_{i-1}}\hat{B}_{z}}(1-r)^{\Sigma_{y\in Y_{i}}\hat{A}_{y}}$
$\displaystyle\quad+\frac{(m_{A}+m_{B})}{2}\delta
t\rho\hat{A}_{i+1}\hat{B}_{i}(1-p)^{\Sigma_{z\in
Z_{i+1}}\hat{A}_{z}}(1-q)^{\Sigma_{y\in Y_{i}}\hat{B}_{y}}(1-r)^{\Sigma_{z\in
Z_{i+1}}\hat{B}_{z}}(1-r)^{\Sigma_{y\in Y_{i}}\hat{A}_{y}}$
$\displaystyle\quad-\frac{(m_{A}+m_{B})}{2}\delta
t\rho\hat{A}_{i}\hat{B}_{i-1}(1-p)^{\Sigma_{z\in
Z_{i}}\hat{A}_{z}}(1-q)^{\Sigma_{y\in Y_{i-1}}\hat{B}_{y}}(1-r)^{\Sigma_{z\in
Z_{i}}\hat{B}_{z}}(1-r)^{\Sigma_{y\in Y_{i-1}}\hat{A}_{y}}$
$\displaystyle\quad-\frac{(m_{A}+m_{B})}{2}\delta
t\rho\hat{A}_{i}\hat{B}_{i+1}(1-p)^{\Sigma_{z\in
Z_{i}}\hat{A}_{z}}(1-q)^{\Sigma_{y\in Y_{i+1}}\hat{B}_{y}}(1-r)^{\Sigma_{z\in
Z_{i}}\hat{B}_{z}}(1-r)^{\Sigma_{y\in Y_{i+1}}\hat{A}_{y}},$ (12)

$\displaystyle\hat{B}_{i}(t+\delta t)=\hat{B}_{i}+\frac{m_{B}}{2}\delta
t\hat{B}_{i-1}(1-\hat{B}_{i}-\hat{A}_{i})(1-q)^{\Sigma_{z\in
Z_{i-1}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i-1}}\hat{A}_{z}}$
$\displaystyle\quad+\frac{m_{B}}{2}\delta
t\hat{B}_{i+1}(1-\hat{B}_{i}-\hat{A}_{i})(1-q)^{\Sigma_{z\in
Z_{i+1}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i+1}}\hat{A}_{z}}$
$\displaystyle\quad-\Big{(}\frac{m_{B}}{2}\delta
t(1-\hat{B}_{i-1}-\hat{A}_{i-1})(1-q)^{\Sigma_{z\in
Z_{i}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i}}\hat{A}_{z}}$
$\displaystyle\quad+\frac{m_{B}}{2}\delta
t(1-\hat{B}_{i+1}-\hat{A}_{i+1})(1-q)^{\Sigma_{z\in
Z_{i}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i}}\hat{A}_{z}}\Big{)}$
$\displaystyle\quad+\frac{(m_{B}+m_{A})}{2}\delta
t\rho\hat{B}_{i-1}\hat{A}_{i}(1-q)^{\Sigma_{z\in
Z_{i-1}}\hat{B}_{z}}(1-p)^{\Sigma_{y\in Y_{i}}\hat{A}_{y}}(1-r)^{\Sigma_{z\in
Z_{i-1}}\hat{A}_{z}}(1-r)^{\Sigma_{y\in Y_{i}}\hat{B}_{y}}$
$\displaystyle\quad+\frac{(m_{B}+m_{A})}{2}\delta
t\rho\hat{B}_{i+1}\hat{A}_{i}(1-q)^{\Sigma_{z\in
Z_{i+1}}\hat{B}_{z}}(1-p)^{\Sigma_{y\in Y_{i}}\hat{A}_{y}}(1-r)^{\Sigma_{z\in
Z_{i+1}}\hat{A}_{z}}(1-r)^{\Sigma_{y\in Y_{i}}\hat{B}_{y}}$
$\displaystyle\quad-\frac{(m_{B}+m_{A})}{2}\delta
t\rho\hat{B}_{i}\hat{A}_{i-1}(1-q)^{\Sigma_{z\in
Z_{i}}\hat{B}_{z}}(1-p)^{\Sigma_{y\in Y_{i-1}}\hat{A}_{y}}(1-r)^{\Sigma_{z\in
Z_{i}}\hat{A}_{z}}(1-r)^{\Sigma_{y\in Y_{i-1}}\hat{B}_{y}}$
$\displaystyle\quad-\frac{(m_{B}+m_{A})}{2}\delta
t\rho\hat{B}_{i}\hat{A}_{i+1}(1-q)^{\Sigma_{z\in
Z_{i}}\hat{B}_{z}}(1-p)^{\Sigma_{y\in Y_{i+1}}\hat{A}_{y}}(1-r)^{\Sigma_{z\in
Z_{i}}\hat{A}_{z}}(1-r)^{\Sigma_{y\in Y_{i+1}}\hat{B}_{y}}.$ (13)

Taylor expanding terms in Equations (II.1) and (II.1) around the site $i$ and
taking $\Delta,\;\delta t\to 0$, such that $\Delta^{2}/\delta t$ remains
constant, leads to the coupled one-dimensional PDEs,

$\displaystyle\frac{\partial A}{\partial t}$
$\displaystyle=\frac{\partial}{\partial
x}\Bigg{[}D_{A}\Big{(}D_{1}(A,B)\frac{\partial B}{\partial
x}+D_{2}(A,B)\frac{\partial A}{\partial
x}\Big{)}+(D_{A}+D_{B})\rho\Big{(}D_{3}(A,B)\frac{\partial A}{\partial
x}-D_{4}(A,B)\frac{\partial B}{\partial x}\Big{)}\Bigg{]},$ (14)

$\displaystyle\frac{\partial B}{\partial t}$
$\displaystyle=\frac{\partial}{\partial
x}\Bigg{[}D_{B}\Big{(}D_{5}(A,B)\frac{\partial A}{\partial
x}+D_{6}(A,B)\frac{\partial B}{\partial
x}\Big{)}+(D_{A}+D_{B})\rho\Big{(}D_{7}(A,B)\frac{\partial B}{\partial
x}-D_{8}(A,B)\frac{\partial A}{\partial x}\Big{)}\Bigg{]},$ (15)

where,

$D_{A}=\lim_{\begin{subarray}{c}\Delta\to 0\\\ \delta t\to
0\end{subarray}}\frac{m_{A}\Delta^{2}}{2\delta t},\qquad
D_{B}=\lim_{\begin{subarray}{c}\Delta\to 0\\\ \delta t\to
0\end{subarray}}\frac{m_{B}\Delta^{2}}{2\delta t},$ (16)

and,

$\displaystyle D_{1}(A,B)$
$\displaystyle=(1-p)^{2A}(1-r)^{2B}(2\ln(1-r)(1-A-B)+1)A,$ (17) $\displaystyle
D_{2}(A,B)$ $\displaystyle=(1-p)^{2A}(1-r)^{2B}(2\ln(1-p)A(1-A-B)+1-B),$ (18)
$\displaystyle D_{3}(A,B)$
$\displaystyle=(1-p)^{2A}(1-q)^{2B}(1-r)^{2A+2B}B(1-2\ln(1-r)A+2\ln(1-p)A),$
(19) $\displaystyle D_{4}(A,B)$
$\displaystyle=(1-p)^{2A}(1-q)^{2B}(1-r)^{2A+2B}A(1-2\ln(1-q)B+2\ln(1-r)B),$
(20) $\displaystyle D_{5}(A,B)$
$\displaystyle=(1-q)^{2B}(1-r)^{2A}(2\ln(1-r)(1-A-B)+1)B,$ (21) $\displaystyle
D_{6}(A,B)$ $\displaystyle=(1-q)^{2B}(1-r)^{2A}(2\ln(1-q)B(1-A-B)+1-A),$ (22)
$\displaystyle D_{7}(A,B)$
$\displaystyle=(1-q)^{2B}(1-p)^{2A}(1-r)^{2A+2B}A(1-2\ln(1-r)B+2\ln(1-q)B),$
(23) $\displaystyle D_{8}(A,B)$
$\displaystyle=(1-q)^{2B}(1-p)^{2A}(1-r)^{2A+2B}B(1-2\ln(1-p)A+2\ln(1-r)A).$
(24)

Figure 4: Comparison between the average two-species ABM and the PDEs (14) and
(15) for set of parameters 4 $p=r=0,\;q=0.25,\;\rho=0.25$ 4,
$p=q=0.25,\;r=0,\;\rho=0.75$ and 4 $p=q=r=0.25,\;\rho=1$. For all the three
figures presented here the movement rates are $m_{A}=m_{B}=1$ and we impose
zero-flux boundary conditions in both the ABM and the PDE. The ABM is shown as
black noisy curves and the PDE solution is presented in red for species A and
green for species B. The domain of length $L=200$ is initialised such that the
central 20% of it has density 1 of type A agents and the rest of the domain is
initialised with type-B agents with a density of 0.5. The ABM is averaged over
500 repeats and we show the density of species A and species B at $t=$ 0, 100
and 1,000. The red and green arrows show the direction of increasing time.

In order to compare our ABM and PDE solutions, we consider a domain of length
$L=200$ and initialise it such that the central 20% of the domain has density
1 of type A agents and the rest of the has density 0.5 of type-B agents, i.e.,

$c_{A}=\begin{cases}1&\text{if }81\leqslant x\leqslant 120,\\\
0&\text{otherwise},\end{cases}$ (25)

and,

$c_{B}=\begin{cases}0.5&\text{if }1\leqslant x\leqslant 80\text{ and
}121\leqslant x\leqslant 200,\\\ 0&\text{otherwise}.\end{cases}$ (26)

Figure 5: Poor agreement between the two-species ABM and the PDEs (14) and
(15). Figure descriptions are as in Figure 4. Parameter values used here are:
5 $p=0.50,\;q=0,\;r=0.25,\;\rho=0.50$ , 5 $p=r=0.25,\;q=0,\;\rho=0.75$, 5
$p=0,\;q=0.5,\;r=0,\;\rho=1$ and 5 $p=0.90,\;q=0.80,\;r=0,\;\rho=1$. The red
and green arrows show the direction of increasing time in panels 5-5.

We let the positions of the agents evolve according to the two-species ABM and
solve the corresponding PDEs numerically with zero flux boundary conditions.
In Figure 4 we present the average density of agents in ABM and the continuum
model at $t=$ 0, 100 and 1,000 for three different sets of low adhesion
strengths and a range of swapping probabilities. Good agreement between the
average ABM and the continuum model can be observed in each case. In Figure
55-5 we see that the agreement is not as good for some other adhesion
strengths due to correlations between sites’ occupancies which the continuum
model ignores. Furthermore, for very high self-adhesion (panel 5) in the ABM,
the overall movement of agents is significantly reduced as a result of fewer
bonds breaking between pairs of neighbouring agents. The PDE, on the other
hand, is unstable for such strong adhesion and consequently we have not been
able to plot its solution. In this particular case the diffusion/cross-
diffusion coefficients $D_{1}(A,B),\;D_{3}(A,B),\;D_{4}(A,B),\;D_{6}(A,B)$ and
$D_{7}(A,b)$ all become negative which implies the potential for aggregation.
While the PDE mimics the average behaviour of the ABM well for some parameter
regimes, it fails to represent the true average dynamics of the ABM in other
cases. In addition, the continuum model is not an appropriate approach for
modelling patterns that emerge as a result of adhesion-mediated cell sorting.

Instead, we can use stochastic mean equations (SMEs) to represent the mean
behaviour of the agents (Thompson  _et al._ , 2012; Anguige and Schmeiser,
2008). We can derive the SMEs directly from Equation (II.1) and (II.1) with
some algebraic manipulation and taking $\delta t\to 0$ which gives us an ODE
for each lattice site $i$ making a system of $L$ equations (one for each site)
for each species:

$\displaystyle\frac{d\hat{A}_{i}(t)}{dt}$ $\displaystyle=\frac{m_{A}}{2}\delta
t\hat{A}_{i-1}(1-\hat{A}_{i}-\hat{B}_{i})(1-p)^{\Sigma_{z\in
Z_{i-1}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i-1}}\hat{B}_{z}}$
$\displaystyle\quad+\frac{m_{A}}{2}\delta
t\hat{A}_{i+1}(1-\hat{A}_{i}-\hat{B}_{i})(1-p)^{\Sigma_{z\in
Z_{i+1}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i+1}}\hat{B}_{z}}$
$\displaystyle\quad-\Big{(}\frac{m_{A}}{2}\delta
t(1-\hat{A}_{i-1}-\hat{B}_{i-1})(1-p)^{\Sigma_{z\in
Z_{i}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i}}\hat{B}_{z}}$
$\displaystyle\quad+\frac{m_{A}}{2}\delta
t(1-\hat{A}_{i+1}-\hat{B}_{i+1})(1-p)^{\Sigma_{z\in
Z_{i}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i}}\hat{B}_{z}}\Big{)}$
$\displaystyle\quad+\frac{(m_{A}+m_{B})}{2}\delta
t\rho\hat{A}_{i-1}\hat{B}_{i}(1-p)^{\Sigma_{z\in
Z_{i-1}}\hat{A}_{z}}(1-q)^{\Sigma_{y\in Y_{i}}\hat{B}_{y}}(1-r)^{\Sigma_{z\in
Z_{i-1}}\hat{B}_{z}}(1-r)^{\Sigma_{y\in Y_{i}}\hat{A}_{y}}$
$\displaystyle\quad+\frac{(m_{A}+m_{B})}{2}\delta
t\rho\hat{A}_{i+1}\hat{B}_{i}(1-p)^{\Sigma_{z\in
Z_{i+1}}\hat{A}_{z}}(1-q)^{\Sigma_{y\in Y_{i}}\hat{B}_{y}}(1-r)^{\Sigma_{z\in
Z_{i+1}}\hat{B}_{z}}(1-r)^{\Sigma_{y\in Y_{i}}\hat{A}_{y}}$
$\displaystyle\quad-\frac{(m_{A}+m_{B})}{2}\delta
t\rho\hat{A}_{i}\hat{B}_{i-1}(1-p)^{\Sigma_{z\in
Z_{i}}\hat{A}_{z}}(1-q)^{\Sigma_{y\in Y_{i-1}}\hat{B}_{y}}(1-r)^{\Sigma_{z\in
Z_{i}}\hat{B}_{z}}(1-r)^{\Sigma_{y\in Y_{i-1}}\hat{A}_{y}}$
$\displaystyle\quad-\frac{(m_{A}+m_{B})}{2}\delta
t\rho\hat{A}_{i}\hat{B}_{i+1}(1-p)^{\Sigma_{z\in
Z_{i}}\hat{A}_{z}}(1-q)^{\Sigma_{y\in Y_{i+1}}\hat{B}_{y}}(1-r)^{\Sigma_{z\in
Z_{i}}\hat{B}_{z}}(1-r)^{\Sigma_{y\in Y_{i+1}}\hat{A}_{y}},$ (27)

$\displaystyle\frac{d\hat{B}_{i}(t)}{dt}$ $\displaystyle=\frac{m_{B}}{2}\delta
t\hat{B}_{i-1}(1-\hat{B}_{i}-\hat{A}_{i})(1-q)^{\Sigma_{z\in
Z_{i-1}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i-1}}\hat{A}_{z}}$
$\displaystyle\quad+\frac{m_{B}}{2}\delta
t\hat{B}_{i+1}(1-\hat{B}_{i}-\hat{A}_{i})(1-q)^{\Sigma_{z\in
Z_{i+1}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i+1}}\hat{A}_{z}}$
$\displaystyle\quad-\Big{(}\frac{m_{B}}{2}\delta
t(1-\hat{B}_{i-1}-\hat{A}_{i-1})(1-q)^{\Sigma_{z\in
Z_{i}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i}}\hat{A}_{z}}$
$\displaystyle\quad+\frac{m_{B}}{2}\delta
t(1-\hat{B}_{i+1}-\hat{A}_{i+1})(1-q)^{\Sigma_{z\in
Z_{i}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i}}\hat{A}_{z}}\Big{)}$
$\displaystyle\quad+\frac{(m_{B}+m_{A})}{2}\delta
t\rho\hat{B}_{i-1}\hat{A}_{i}(1-q)^{\Sigma_{z\in
Z_{i-1}}\hat{B}_{z}}(1-p)^{\Sigma_{y\in Y_{i}}\hat{A}_{y}}(1-r)^{\Sigma_{z\in
Z_{i-1}}\hat{A}_{z}}(1-r)^{\Sigma_{y\in Y_{i}}\hat{B}_{y}}$
$\displaystyle\quad+\frac{(m_{B}+m_{A})}{2}\delta
t\rho\hat{B}_{i+1}\hat{A}_{i}(1-q)^{\Sigma_{z\in
Z_{i+1}}\hat{B}_{z}}(1-p)^{\Sigma_{y\in Y_{i}}\hat{A}_{y}}(1-r)^{\Sigma_{z\in
Z_{i+1}}\hat{A}_{z}}(1-r)^{\Sigma_{y\in Y_{i}}\hat{B}_{y}}$
$\displaystyle\quad-\frac{(m_{B}+m_{A})}{2}\delta
t\rho\hat{B}_{i}\hat{A}_{i-1}(1-q)^{\Sigma_{z\in
Z_{i}}\hat{B}_{z}}(1-p)^{\Sigma_{y\in Y_{i-1}}\hat{A}_{y}}(1-r)^{\Sigma_{z\in
Z_{i}}\hat{A}_{z}}(1-r)^{\Sigma_{y\in Y_{i-1}}\hat{B}_{y}}$
$\displaystyle\quad-\frac{(m_{B}+m_{A})}{2}\delta
t\rho\hat{B}_{i}\hat{A}_{i+1}(1-q)^{\Sigma_{z\in
Z_{i}}\hat{B}_{z}}(1-p)^{\Sigma_{y\in Y_{i+1}}\hat{A}_{y}}(1-r)^{\Sigma_{z\in
Z_{i}}\hat{A}_{z}}(1-r)^{\Sigma_{y\in Y_{i+1}}\hat{B}_{y}}.$ (28)

Figure 6: Comparison between the average two-species ABM and the numerical
solution of the SMEs (II.1) and (II.1) for the same parameters as Figure 55-5
shown here in panels 6-6, respectively. The initial condition is the same as
Figure 5.

Figure 7: Pattern formation in a single realisation of the ABM. The one-
dimensional lattice with $L=100$ sites is initialised with agents of type A
and type B. In 7 species A and species B are initially segregated. Each site
on the left half of the lattice (corresponding to $1\leqslant i\leqslant 50$)
is seeded with a type-A agent with probability $\pi_{1}=0.9$ or the site is
left empty otherwise. Each site on the right half of the lattice
(corresponding to $51\leqslant i\leqslant 100$) is seeded with a type-B agent
with probability $\pi_{2}=0.9$ or left empty otherwise. In 7 and 7 each site
is occupied with a type-A agent with probability $\pi_{1}=0.45$ or a type-B
agent with probability $\pi_{2}=0.45$. A site is left empty with probability
$1-\pi_{1}-\pi_{2}$. We let the positions of the agents evolve according to
the two-species ABM and present the densities of species A and B at $t=$[0,
1,000, 100,000]. The three simulations evolve to three distinct types of
patterns: mixing for $p=0.6,\,q=0.55,\,r=0.5$ in panels 7-7, aggregation for
$p=0.9,\,q=0.3,\,r=0$ in panels 8-8 and one-dimensional checkerboard (i.e.
alternating) for $p=0,\,q=0,\,r=0.9$ in panels 8-8. In all cases, we assume
the rate of movements $m_{A}=m_{B}=1$ and $\rho=1$ with reflecting boundaries.

In Figure 6 we compare the average behaviour of the ABM with the numerical
solution of the SMEs for the same parameters as Figure 55-5 (for which we saw
poor agreement between the ABM and the PDE previously). We can immediately see
that the agreement between the ABM and SMEs is better than the PDEs as the
SMEs more closely resemble the average behaviour of the ABM.

Figure 8: Pattern formation in the SMEs. The one-dimensional domain with
$L=100$ sites is initialised with mean densities of species A and B. In 7
species A and species B are initially segregated where each site on the left
half of the lattice (corresponding to $1\leqslant i\leqslant 50$) is given a
mean density $c_{A}=0.9$ of species A and each site on the right half of the
lattice (corresponding to $51\leqslant i\leqslant 100$) is assigned a mean
density $c_{B}=0.9$ of species B. In 7 and 7 each lattice site is assigned a
mean density of $c_{A}=0.45$ and $c_{B}=0.45$. In all three scenarios, the
initial density was perturbed by adding small independent normally distributed
random variables (with mean 0 and variance $0.01$) to the mean density of each
species at each lattice site (for symmetry breaking). The SMEs were solved
numerically with zero-flux boundary conditions and densities of species A and
B at $t=$[0, 1,000, 100,000] showing evolution to three types of patterns are
presented here: mixing for $p=0.6,\,q=0.55,\,r=0.5$ in panels 7-7, aggregation
for $p=0.9,\,q=0.3,\,r=0$ in panels 8-8 and checkerboard for
$p=0,\,q=0,\,r=0.9$ in panels 8-81. In all cases, we assume the rate of
movements $m_{A}=m_{B}=1$ and $\rho=1$.

To explore whether our ABM and SMEs are capable of replicating the sorting
behaviour of biological cells, we investigate the ABM and the SMEs with
different sets of adhesion parameters. We let the positions of the agents
evolve according to the ABM and solve the corresponding SMEs numerically with
reflecting boundaries. In Figures 7 and 8 we exemplify three types of
behaviours shown by our model in the ABM and the SME (initial data and
parameters in figure caption). In panels 8-8 subject to a mixture of high to
moderate adhesion strengths ($p=0.6,\;q=0.55,\;r=0.5$) we see mixing of the
two populations to a uniform steady state. More complex patterns are observed
in panels 8-8, in which high self-adhesion between species A and low self-
adhesion between species B gives rise to out-of-phase peaks and troughs of
density characterising aggregation. In panels 8-8 high cross adhesion
($r=0.90$) shows cell sorting behaviour in which the density of lattice sites
alternates between species A and B resembling a one-dimensional analogue of a
checkerboard pattern. In all three cases we have used a swapping probability
of 1 since we found that varying the swapping probability did not
fundamentally change the pattern in the SMEs and a high swapping probability
led to faster convergence to a pattern.

11footnotetext: In all the panels red denotes the mean density of species A
and light green denotes that of species B. In panels 7-7 and 7-7 dark green
represents the overlap between the mean densities of the two species.

###  II.2 Two-dimensional model

In the two-dimensional ABM the agents move on a two-dimensional grid that has
$L_{x}$ sites in the horizontal direction and $L_{y}$ sites in the vertical
direction. Naturally, agents on the two-dimensional grid have more neighbours
(up, down, left and right) than on a one-dimensional lattice (left and right)
which in turn reduces the probability of an agent breaking free from their
neighbours due to more adhesive forces from the additional neighbours.

Figure 9: A schematic illustrating an example of the two-dimensional ABM with
adhesion. A red site is occupied by a type-A agent, a green site represents a
type-B agent and a black site is empty. Adhesion strengths between
neighbouring pairs of agents are shown as letters $p$, $q$ and $r$ alongside a
dash representing a bond between a pair of agents. For clarity, we have only
shown labels for three pairs but adhesion is assumed between all neighbouring
pairs of agents. The initial configuration of the lattice is shown in panel 9.
The arrow symbolises the red agent being chosen to move to a site to the
right. In panel 9 we show the result of a successful swap. Note the labels
denoting the adhesive bond between pairs of neighbouring agents have also been
updated.

In Figure 9 we show an example of an agent swapping into an already occupied
site to the right of it on a two-dimensional lattice. The moving agent is of
type A (red) with one type-A agent and two type-B agents (green) in its
neighbourhood, thus the probability of it breaking adhesive bonds with its
neighbours is $(1-p)(1-r)^{2}$. The agent at the target site is a type-B agent
with 1 type-B agent and two type-A agents in its neighbourhood and therefore
the probability of the agent at the target site breaking adhesive bonds with
its neighbours is $(1-q)(1-r)^{2}$. Thus the probability that the swap is
successful is $\rho(1-p)(1-q)(1-r)^{4}$. In the case of a successful swap the
two agents make new connections with their new neighbours (see panel 9).

The PDEs for the densities of species A and species B in the two-dimensional
setting can be derived from the corresponding two-dimensional master equations
in a similar way to the one-dimensional PDEs. The two-dimensional master
equations for species A and B can be found in Appendix A.

The continuum PDEs describing the evolution of the densities of species A and
species B are given by,

$\displaystyle\frac{\partial A}{\partial
t}=\nabla\cdot\Bigg{[}D_{A}\Big{(}D_{1}^{\text{2D}}(A,B)\nabla
B+D_{2}^{\text{2D}}(A,B)\nabla A\Big{)}$
$\displaystyle\quad+(D_{A}+D_{B})\rho\Big{(}D_{3}^{\text{2D}}(A,B)\nabla
A-D_{4}^{\text{2D}}(A,B)\nabla B\Big{)}\Bigg{]},$ (29)

$\displaystyle\frac{\partial B}{\partial
t}=\nabla\cdot\Bigg{[}D_{B}\Big{(}D_{5}^{\text{2D}}(A,B)\nabla
A+D_{6}^{\text{2D}}(A,B)\nabla B\Big{)}$
$\displaystyle\quad+(D_{A}+D_{B})\rho\Big{(}D_{7}^{\text{2D}}(A,B)\nabla
B-D_{8}^{\text{2D}}(A,B)\nabla A\Big{)}\Bigg{]},$ (30)

where,

$D_{A}=\lim_{\begin{subarray}{c}\Delta\to 0\\\ \delta t\to
0\end{subarray}}\frac{m_{A}}{4}\frac{\Delta^{2}}{\delta t},\quad
D_{B}=\lim_{\begin{subarray}{c}\Delta\to 0\\\ \delta t\to
0\end{subarray}}\frac{m_{B}}{4}\frac{\Delta^{2}}{\delta t},$ (31)

and,

$\displaystyle D_{1}^{\text{2D}}(A,B)$
$\displaystyle=(1-p)^{4A}(1-r)^{4B}(4\ln(1-r)(1-A-B)+1)A,$ $\displaystyle
D_{2}^{\text{2D}}(A,B)$
$\displaystyle=(1-p)^{4A}(1-r)^{4B}(4\ln(1-p)A(1-A-B)+1-B),$ $\displaystyle
D_{3}^{\text{2D}}(A,B)$
$\displaystyle=(1-p)^{4A}(1-q)^{4B}(1-r)^{4A+4B}B(1-4\ln(1-r)A+4\ln(1-p)A),$
$\displaystyle D_{4}^{\text{2D}}(A,B)$
$\displaystyle=(1-p)^{4A}(1-q)^{4B}(1-r)^{4A+4B}A(1-4\ln(1-q)B+4\ln(1-r)B),$
$\displaystyle D_{5}^{\text{2D}}(A,B)$
$\displaystyle=(1-q)^{4B}(1-r)^{4A}(4\ln(1-r)(1-A-B)+1)B,$ $\displaystyle
D_{6}^{\text{2D}}(A,B)$
$\displaystyle=(1-q)^{4B}(1-r)^{4A}(4\ln(1-q)B(1-A-B)+1-A),$ $\displaystyle
D_{7}^{\text{2D}}(A,B)$
$\displaystyle=(1-q)^{4B}(1-p)^{4A}(1-r)^{4A+4B}(A(1-4\ln(1-r)B+4\ln(1-q)B),$
$\displaystyle D_{8}^{\text{2D}}(A,B)$
$\displaystyle=(1-q)^{4B}(1-p)^{4A}(1-r)^{4A+4B}B(1-4\ln(1-p)A+4\ln(1-r)A).$
(32)

Figure 10: Comparison between the two species column-averaged ABM and the
vertically-averaged solution of PDEs (II.2) and (II.2) at $t=0,\,100,\,1000$
for $m_{A}=m_{B}=1$ and 10 $(p,q,r,\rho)=(0,0.25,0,0.25)$, 10
$(p,q,r,\rho)=(0.25,0.5,0.25,1)$ and 10 $(p,q,r,\rho)=(0.25,0.5,0.5,0.75)$.
The lattice has dimensions $L_{x}=200$ columns and $L_{y}=20$ rows with zero-
flux boundary conditions and is initialised such that all the sites for
$81\leqslant x\leqslant 120$ have density 1 for species A and the remainder of
the domain is initialised with species B at density 0.5.

Figure 11: Engulfment 11-11 and checkerboard 11-11 patterns in the two-
dimensional ABM model on a 20 by 20 lattice from a uniformly random initial
condition with reflecting boundaries and $m_{A}=m_{B}=\rho=1$. Each lattice
site was seeded with a type-A agent with probability $\pi_{1}$ or a type-B
agent with probability $\pi_{2}=1-\pi_{1}$. For engulfment in 11 $\pi_{1}=0.7$
and hence $\pi_{2}=0.3$ and for the checkerboard simulation in 11
$\pi_{1}=\pi_{2}=0.5$. We present lattice snapshots at $t=$[0, 10,000,
100,000] for engulfment 11-11 for $(p,q,r)=(0.2,0.9,0)$ and checkerboard 11-11
for $(p,q,r)=(0,0,0.8)$.

In Figure 10 we compare the average column density of the ABM with the
vertically-averaged solution of PDEs (II.2) and (II.2) over a lattice with
$L_{x}=200$ and $L_{y}=20$ with zero-flux boundaries and movement rates
$m_{A}=m_{B}=1$. Lattice is initialised such that all the sites in the middle
20% of the domain are fully occupied by type-A agents and the rest of the
sites are seeded with type-B agents at density 0.5. For ease of comparison, at
$t=0,\,100$ and 1,000 the densities of the two-dimensional ABM were averaged
over the vertical extent of the domain to obtain a one-dimensional
representation of the densities of species A and B.

Likewise, the two-dimensional PDEs (II.2) and (II.2) were solved numerically
for the densities of both species which were then vertically averaged to
obtain one-dimensional analogues. We notice good agreement between the ABM and
the PDE for a range of swapping and adhesion parameters. However, we also
found that for very strong adhesion, as in the one-dimensional scenario, the
PDEs become unstable (figure not shown). Consequently, it is not possible to
observe patterns that require strong adhesion in the PDE model. As with the
one-dimensional model, we can use the SMEs corresponding to the two-
dimensional ABM to analyse our model for pattern formation.

Similar to the one-dimensional model, the SMEs can be derived directly from
the two-dimensional master equations for species A and B (see Appendix B). In
Figures 11 and 12 we present examples of engulfment and checkerboard patterns
obtained with our adhesion model under the ABM and the SMEs, respectively, on
a 20 by 20 domain with zero-flux boundaries.

Figure 12: Engulfment and checkerboard patterns in the two-dimensional SME
model on a 20 by 20 lattice from a uniformly random initial condition with
reflecting boundaries with $m_{A}=m_{B}=\rho=1$. In both cases, each lattice
site was initialised with a constant mean density $c_{A}$ of species A, which
was perturbed by adding independent normally distributed random variables
(with zero mean and variance 0.01) to give the perturbed density of species A,
$\tilde{c}_{A}$. The density of species B at each site was then assigned as
$c_{B}=1-\tilde{c}_{A}$. For engulfment, $c_{A}=0.7$ and for the checkerboard
$c_{A}=0.5$. We present the solution of the SMEs at $t=$[0, 10,000, 100,000]
for engulfment 12-12 for $(p,q,r)=(0.2,0.9,0)$ and the checkerboard 12-12 for
$(p,q,r)=(0,0,0.8)$. Due to space constraints, the colour bar is exclusively
displayed in the rightmost panel for each scenario. Nevertheless, it remains
consistent across all panels for that specific scenario.

In panels 11-11 of Figures 11 and 12 strong self-adhesion between agents of
type A, low self-adhesion between agents of type B and no cross-adhesion leads
to engulfment of species A by species B over time. Starting from uniformly
random densities, species A slowly engulfs species B over time. In panels
11-11 high cross-adhesion and no self-adhesion leads to the self-organisation
of cells into a checkerboard pattern. The agreement between the ABM and the
SME patterns is commendable for the intermediate as well as (and especially)
the final pattern. The engulfment and checkerboard patterns seen here are
comparable to the type of engulfment observed in the cultures of dissociated
chick embryo pigmented retinal epithelial and neural retinal cells in which
the neural retinal cells engulf the pigmented retinal epithelial cells
(Armstrong, 1971) and the checkerboard tiling on the luminal surface of the
mature oviduct epithelium in adult Japanese quail (Honda  _et al._ , 1986).

##  III Conclusion

In this paper, we have developed a model for cell migration incorporating
cell-cell adhesion. By first describing the ABM for the dynamics of individual
cells, we derived a set of partial differential equations (PDEs) describing
the macroscopic behaviour of the ABM. The work presented here is primarily a
response to the lack of multispecies discrete and continuum models with local
adhesive interactions mediating contact inhibition and serves to elucidate the
proof-of-concept model we have developed. In both, one-dimensional and two-
dimensional models, presented in Sections II.1 and II.2, respectively, we
showed that the microscopic and macroscopic modelling paradigms agree well
with each other for some parameter values but do not for others. Particularly,
we notice that the PDE descriptions are incapable of replicating adhesion-
mediated cell sorting patterns. To overcome this limitation, we derived a set
of stochastic mean equations (SMEs) to represent the mean density of cells as
an alternative to the PDEs. We found that not only were the SMEs able to
reproduce biologically realistic patterns, but they also showed agreement
between the ABM and the SME for certain parameter regimes in which the ABM and
the PDEs did not agree well.

We found that our discrete model is capable of replicating biologically
realistic patterns such as engulfment and checker-board which are comparable
to the engulfment of chick embryo pigmented retinal epithelial cells by neural
retinal cells in culture (Armstrong, 1971) and the checkerboard pattern in
Japanese quail oviduct epithelium eggs (Honda  _et al._ , 1986), respectively.
Our results accord with the differential adhesion hypothesis (Steinberg,
1963); when two populations with different of adhesive strengths interact,
cells with stronger adhesion will prefer to adhere to each other and exclude
the species with weaker adhesion as this results in more energetically
favourable configurations (Dickinson, 2020). In the case of the checkerboard
pattern (see Figure 11 and 12), the higher cross-adhesion between the two
populations leads to the preferential attachment of one cell type to the other
indirectly repelling cells of the same type leading to the energetically
favourable state resembling a checkerboard.

Our model has the potential to elucidate cell adhesion dynamics in a two-
species environment, applicable across a spectrum of biological contexts where
cell-cell adhesion plays a role in driving cell cohesiveness such as wound
healing, angiogenesis and the separation of cells in germ layers endoderm,
mesoderm and ectoderm in the embryo (Dickinson, 2020; Schötz  _et al._ ,
2008).

Although we have included realism in our models by incorporating cell-cell
interactions such as adhesion, volume filling and cell-cell swapping, we are
neglecting the fact that cell populations are not constant over time. A
natural extension of this study could be to allow cell division and analyse
its effects. This could be achieved using a multistage cell-cycle model for
improved realism (Yates  _et al._ , 2017; Smith and Martin, 1973; Weber  _et
al._ , 2014; Mort  _et al._ , 2014; Howard and Pelc, 1951). Another limitation
is the fixed size of the domain over which cell movement is modelled. In
reality, the domain (representing a developing embryo or growing tissue, for
example) might grow over time. Incorporating domain growth and cell division
into our model would improve the applicability of our work and help us
understand the role of adhesion in a wider variety of biological processes.

##  Appendix A Two-dimensional partial differential equations (PDEs)

Let $\hat{A}_{ij}(t)$ be the continuous approximation of the mean density of
species A at the site $(i,j)$ at time $t$ and let $\hat{B}_{ij}(t)$ be the be
the equivalent for species B. Then the master equation for the continuous
approximation of the mean occupancy at the site $(i,j)$ and time $t+\delta t$
where $\delta t$ is an infinitesimally small change in time is given by,

$\displaystyle\hat{A}_{ij}(t+\delta t)=A_{ij}+\frac{m_{A}}{4}\delta
t(1-\hat{A}_{ij}-\hat{B}_{ij})\Big{[}\hat{A}_{i-1,j}(1-p)^{\Sigma_{z\in
Z_{i-1,j}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i-1,j}}\hat{B}_{z}}$
$\displaystyle\quad+\hat{A}_{i+1,j}(1-p)^{\Sigma_{z\in
Z_{i+1,j}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in
Z_{i+1,j}}\hat{B}_{z}}+\hat{A}_{i,j-1}(1-p)^{\Sigma_{z\in
Z_{i,j-1}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j-1}}\hat{B}_{z}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\hat{A}_{i,j+1}(1-p)^{\Sigma_{z\in
Z_{i,j+1}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j+1}}\hat{B}_{z}}\Big{]}$
$\displaystyle\quad-\frac{m_{A}}{4}\delta t\hat{A}_{ij}(1-p)^{\Sigma_{z\in
Z_{ij}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in
Z_{ij}}\hat{B}_{z}}\Big{[}4-\hat{A}_{i-1,j}-\hat{B}_{i-1,j}-\hat{A}_{i+1,j}-\hat{B}_{i+1,j}$
$\displaystyle\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\hat{A}_{i,j-1}-\hat{B}_{i,j-1}-\hat{A}_{i,j+1}-\hat{B}_{i,j+1}\Big{]}$
$\displaystyle\quad+\frac{(m_{A}+m_{B})}{4}\delta t\rho
B_{ij}(1-q)^{\Sigma_{y\in Y_{ij}}\hat{B}_{z}}(1-r)^{\Sigma_{y\in
Y_{ij}}\hat{A}_{z}}\Big{[}\hat{A}_{i-1,j}(1-p)^{\Sigma_{z\in
Z_{i-1,j}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i-1,j}}\hat{B}_{z}}$
$\displaystyle\quad+\hat{A}_{i+1,j}(1-p)^{\Sigma_{z\in
Z_{i+1,j}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in
Z_{i+1,j}}\hat{B}_{z}}+\hat{A}_{i,j-1}(1-p)^{\Sigma_{z\in
Z_{i,j-1}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j-1}}\hat{B}_{z}}$
$\displaystyle\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\hat{A}_{i,j+1}(1-p)^{\Sigma_{z\in
Z_{i,j+1}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j+1}}\hat{B}_{z}}\Big{]}$
$\displaystyle\quad-\frac{(m_{A}+m_{B})}{4}\delta
t\rho\hat{A}_{i,j}(1-p)^{\Sigma_{z\in Z_{ij}}\hat{A}_{z}}(1-r)^{\Sigma_{y\in
Y_{ij}}\hat{B}_{z}}\Big{[}\hat{B}_{i-1,j}(1-q)^{\Sigma_{y\in
Y_{i-1,j}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i-1,j}}\hat{A}_{z}}$
$\displaystyle\quad+\hat{B}_{i+1,j}(1-q)^{\Sigma_{y\in
Y_{ij}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in
Z_{i+1,j}}\hat{A}_{z}}+\hat{B}_{i,j-1}(1-q)^{\Sigma_{y\in
Y_{ij}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j-1}}\hat{A}_{z}}$
$\displaystyle\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\hat{B}_{i,j+1}(1-q)^{\Sigma_{y\in
Y_{ij}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j+1}}\hat{A}_{z}}\Big{]},$ (33)

$\displaystyle\hat{B}_{ij}(t+\delta t)=B_{ij}+\frac{m_{B}}{4}\delta
t(1-\hat{B}_{ij}-\hat{A}_{ij})\Big{[}\hat{B}_{i-1,j}(1-q)^{\Sigma_{z\in
Z_{i-1,j}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i-1,j}}\hat{A}_{z}}$
$\displaystyle\quad+\hat{B}_{i+1,j}(1-q)^{\Sigma_{z\in
Z_{i+1,j}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in
Z_{i+1,j}}\hat{A}_{z}}+\hat{B}_{i,j-1}(1-q)^{\Sigma_{z\in
Z_{i,j-1}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j-1}}\hat{A}_{z}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\hat{B}_{i,j+1}(1-q)^{\Sigma_{z\in
Z_{i,j+1}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j+1}}\hat{A}_{z}}\Big{]}$
$\displaystyle\quad-\frac{m_{B}}{4}\delta t\hat{B}_{ij}(1-q)^{\Sigma_{z\in
Z_{ij}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in
Z_{ij}}\hat{A}_{z}}\Big{[}4-\hat{B}_{i-1,j}-\hat{A}_{i-1,j}-\hat{B}_{i+1,j}-\hat{A}_{i+1,j}$
$\displaystyle\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\hat{B}_{i,j-1}-\hat{A}_{i,j-1}-\hat{B}_{i,j+1}-\hat{A}_{i,j+1}\Big{]}$
$\displaystyle\quad+\frac{(m_{B}+m_{A})}{4}\delta t\rho
A_{ij}(1-p)^{\Sigma_{y\in Y_{ij}}\hat{A}_{z}}(1-r)^{\Sigma_{y\in
Y_{ij}}\hat{B}_{z}}\Big{[}\hat{B}_{i-1,j}(1-q)^{\Sigma_{z\in
Z_{i-1,j}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i-1,j}}\hat{A}_{z}}$
$\displaystyle\quad+\hat{B}_{i+1,j}(1-q)^{\Sigma_{z\in
Z_{i+1,j}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in
Z_{i+1,j}}\hat{A}_{z}}+\hat{B}_{i,j-1}(1-q)^{\Sigma_{z\in
Z_{i,j-1}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j-1}}\hat{A}_{z}}$
$\displaystyle\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\hat{B}_{i,j+1}(1-q)^{\Sigma_{z\in
Z_{i,j+1}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j+1}}\hat{A}_{z}}\Big{]}$
$\displaystyle\quad-\frac{(m_{B}+m_{A})}{4}\delta
t\rho\hat{B}_{i,j}(1-q)^{\Sigma_{z\in Z_{ij}}\hat{B}_{z}}(1-r)^{\Sigma_{y\in
Y_{ij}}\hat{A}_{z}}\Big{[}\hat{A}_{i-1,j}(1-p)^{\Sigma_{y\in
Y_{i-1,j}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i-1,j}\hat{B}_{z}}}$
$\displaystyle\quad+\hat{A}_{i+1,j}(1-p)^{\Sigma_{y\in
Y_{ij}\hat{A}_{z}}}(1-r)^{\Sigma_{z\in
Z_{i+1,j}\hat{B}_{z}}}+\hat{A}_{i,j-1}(1-p)^{\Sigma_{y\in
Y_{ij}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j-1}}\hat{B}_{z}}$
$\displaystyle\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\hat{A}_{i,j+1}(1-p)^{\Sigma_{y\in
Y_{ij}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j+1}}\hat{B}_{z}}\Big{]}.$ (34)

Taylor expanding the terms in the above equations around $i$ and taking the
diffusive limit leads to the coupled PDEs,

$\displaystyle\frac{\partial A}{\partial
t}=\nabla\cdot\Bigg{[}D_{A}\Big{(}D_{1}^{\text{2D}}(A,B)\nabla
B+D_{2}^{\text{2D}}(A,B)\nabla A\Big{)}$
$\displaystyle\quad+(D_{A}+D_{B})\rho\Big{(}D_{3}^{\text{2D}}(A,B)\nabla
A-D_{4}^{\text{2D}}(A,B)\nabla B\Big{)}\Bigg{]},$

$\displaystyle\frac{\partial B}{\partial
t}=\nabla\cdot\Bigg{[}D_{B}\Big{(}D_{5}^{\text{2D}}(A,B)\nabla
A+D_{6}^{\text{2D}}(A,B)\nabla B\Big{)}$
$\displaystyle\quad+(D_{A}+D_{B})\rho\Big{(}D_{7}^{\text{2D}}(A,B)\nabla
B-D_{8}^{\text{2D}}(A,B)\nabla A\Big{)}\Bigg{]},$

where $D_{A}$, $D_{B}$ are given by Equation (31) and
$D_{1}^{\text{2D}},...,D_{8}^{\text{2D}}$ by Equation set (32) in the main
text.

##  Appendix B Two-dimensional stochastic mean equations

The two-dimensional stochastic mean equations (SMEs) can be derived directly
from the master equations (A) and (A) by taking away $A_{ij}$ from both sides,
dividing both sides by $\delta t$ and taking the limit as $\delta t\to 0$.
This leads to $L_{x}\times L_{y}$ ordinary differential equations (ODEs) (one
for each lattice site $(i,j)$),

$\displaystyle\frac{\mathrm{d}\hat{A}_{ij}(t)}{\mathrm{d}t}=\frac{m_{A}}{4}(1-\hat{A}_{ij}-\hat{B}_{ij})\Big{[}\hat{A}_{i-1,j}(1-p)^{\Sigma_{z\in
Z_{i-1,j}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i-1,j}}\hat{B}_{z}}$
$\displaystyle\quad+\hat{A}_{i+1,j}(1-p)^{\Sigma_{z\in
Z_{i+1,j}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in
Z_{i+1,j}}\hat{B}_{z}}+\hat{A}_{i,j-1}(1-p)^{\Sigma_{z\in
Z_{i,j-1}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j-1}}\hat{B}_{z}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\hat{A}_{i,j+1}(1-p)^{\Sigma_{z\in
Z_{i,j+1}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j+1}}\hat{B}_{z}}\Big{]}$
$\displaystyle\quad-\frac{m_{A}}{4}\hat{A}_{ij}(1-p)^{\Sigma_{z\in
Z_{ij}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in
Z_{ij}}\hat{B}_{z}}\Big{[}4-\hat{A}_{i-1,j}-\hat{B}_{i-1,j}-\hat{A}_{i+1,j}-\hat{B}_{i+1,j}$
$\displaystyle\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\hat{A}_{i,j-1}-\hat{B}_{i,j-1}-\hat{A}_{i,j+1}-\hat{B}_{i,j+1}\Big{]}$
$\displaystyle\quad+\frac{(m_{A}+m_{B})}{4}\rho B_{ij}(1-q)^{\Sigma_{y\in
Y_{ij}}\hat{B}_{z}}(1-r)^{\Sigma_{y\in
Y_{ij}}\hat{A}_{z}}\Big{[}\hat{A}_{i-1,j}(1-p)^{\Sigma_{z\in
Z_{i-1,j}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i-1,j}}\hat{B}_{z}}$
$\displaystyle\quad+\hat{A}_{i+1,j}(1-p)^{\Sigma_{z\in
Z_{i+1,j}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in
Z_{i+1,j}}\hat{B}_{z}}+\hat{A}_{i,j-1}(1-p)^{\Sigma_{z\in
Z_{i,j-1}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j-1}}\hat{B}_{z}}$
$\displaystyle\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\hat{A}_{i,j+1}(1-p)^{\Sigma_{z\in
Z_{i,j+1}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j+1}}\hat{B}_{z}}\Big{]}$
$\displaystyle\quad-\frac{(m_{A}+m_{B})}{4}\rho\hat{A}_{i,j}(1-p)^{\Sigma_{z\in
Z_{ij}}\hat{A}_{z}}(1-r)^{\Sigma_{y\in
Y_{ij}}\hat{B}_{z}}\Big{[}\hat{B}_{i-1,j}(1-q)^{\Sigma_{y\in
Y_{i-1,j}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i-1,j}}\hat{A}_{z}}$
$\displaystyle\quad+\hat{B}_{i+1,j}(1-q)^{\Sigma_{y\in
Y_{ij}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in
Z_{i+1,j}}\hat{A}_{z}}+\hat{B}_{i,j-1}(1-q)^{\Sigma_{y\in
Y_{ij}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j-1}}\hat{A}_{z}}$
$\displaystyle\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\hat{B}_{i,j+1}(1-q)^{\Sigma_{y\in
Y_{ij}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j+1}}\hat{A}_{z}}\Big{]},$ (35)

$\displaystyle\frac{\mathrm{d}\hat{B}_{ij}(t)}{\mathrm{d}t}=\frac{m_{B}}{4}(1-\hat{B}_{ij}-\hat{A}_{ij})\Big{[}\hat{B}_{i-1,j}(1-q)^{\Sigma_{z\in
Z_{i-1,j}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i-1,j}}\hat{A}_{z}}$
$\displaystyle\quad+\hat{B}_{i+1,j}(1-q)^{\Sigma_{z\in
Z_{i+1,j}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in
Z_{i+1,j}}\hat{A}_{z}}+\hat{B}_{i,j-1}(1-q)^{\Sigma_{z\in
Z_{i,j-1}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j-1}}\hat{A}_{z}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\hat{B}_{i,j+1}(1-q)^{\Sigma_{z\in
Z_{i,j+1}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j+1}}\hat{A}_{z}}\Big{]}$
$\displaystyle\quad-\frac{m_{B}}{4}\hat{B}_{ij}(1-q)^{\Sigma_{z\in
Z_{ij}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in
Z_{ij}}\hat{A}_{z}}\Big{[}4-\hat{B}_{i-1,j}-\hat{A}_{i-1,j}-\hat{B}_{i+1,j}-\hat{A}_{i+1,j}$
$\displaystyle\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\hat{B}_{i,j-1}-\hat{A}_{i,j-1}-\hat{B}_{i,j+1}-\hat{A}_{i,j+1}\Big{]}$
$\displaystyle\quad+\frac{(m_{B}+m_{A})}{4}\rho A_{ij}(1-p)^{\Sigma_{y\in
Y_{ij}}\hat{A}_{z}}(1-r)^{\Sigma_{y\in
Y_{ij}}\hat{B}_{z}}\Big{[}\hat{B}_{i-1,j}(1-q)^{\Sigma_{z\in
Z_{i-1,j}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i-1,j}}\hat{A}_{z}}$
$\displaystyle\quad+\hat{B}_{i+1,j}(1-q)^{\Sigma_{z\in
Z_{i+1,j}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in
Z_{i+1,j}}\hat{A}_{z}}+\hat{B}_{i,j-1}(1-q)^{\Sigma_{z\in
Z_{i,j-1}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j-1}}\hat{A}_{z}}$
$\displaystyle\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\hat{B}_{i,j+1}(1-q)^{\Sigma_{z\in
Z_{i,j+1}}\hat{B}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j+1}}\hat{A}_{z}}\Big{]}$
$\displaystyle\quad-\frac{(m_{B}+m_{A})}{4}\rho\hat{B}_{i,j}(1-q)^{\Sigma_{z\in
Z_{ij}}\hat{B}_{z}}(1-r)^{\Sigma_{y\in
Y_{ij}}\hat{A}_{z}}\Big{[}\hat{A}_{i-1,j}(1-p)^{\Sigma_{y\in
Y_{i-1,j}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i-1,j}\hat{B}_{z}}}$
$\displaystyle\quad+\hat{A}_{i+1,j}(1-p)^{\Sigma_{y\in
Y_{ij}\hat{A}_{z}}}(1-r)^{\Sigma_{z\in
Z_{i+1,j}\hat{B}_{z}}}+\hat{A}_{i,j-1}(1-p)^{\Sigma_{y\in
Y_{ij}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j-1}}\hat{B}_{z}}$
$\displaystyle\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\hat{A}_{i,j+1}(1-p)^{\Sigma_{y\in
Y_{ij}}\hat{A}_{z}}(1-r)^{\Sigma_{z\in Z_{i,j+1}}\hat{B}_{z}}\Big{]}.$ (36)

## References

  * Loeb (1921) L. Loeb, Science 53, 261 (1921), 1645041 .
  * Abercrombie and Heaysman (1954) M. Abercrombie and J. E. M. Heaysman, Exp. Cell Res. 6, 293 (1954).
  * Abercrombie (1979) M. Abercrombie, Nature 281, 259 (1979).
  * Abercrombie (1997) M. Abercrombie, Proc. R. Soc. Lond. B Biol. Sci. 207, 129 (1997).
  * Lee  _et al._ (1995) Y. Lee, S. Kouvroukoglou, L. V. McIntire, and K. Zygourakis, Biophys J 69, 1284 (1995).
  * Carmona-Fontaine  _et al._ (2008) C. Carmona-Fontaine, H. K. Matthews, S. Kuriyama, M. Moreno, G. A. Dunn, M. Parsons, C. D. Stern, and R. Mayor, Nature 456, 957 (2008).
  * Deutsch  _et al._ (2021) A. Deutsch, J. M. Nava-Sedeño, S. Syga, and H. Hatzikirou, PLoS Comput. Biol. 17, e1009066 (2021).
  * Yates  _et al._ (2015) C. A. Yates, A. Parker, and R. E. Baker, Phys. Rev. E 91, 052711 (2015).
  * Yamanaka and Kondo (2014) H. Yamanaka and S. Kondo, Proc. Natl. Acad. Sci. U. S. A. 111, 1867 (2014).
  * Davis and Trinkaus (1981) E. M. Davis and J. P. Trinkaus, Development 63, 29 (1981).
  * Chappelle and Yates (2019) G. Chappelle and C. A. Yates, Phys. Rev. E 99, 062413 (2019).
  * Raja Noureen  _et al._ (2023) S. Raja Noureen, J. P. Owen, R. L. Mort, and C. A. Yates, Phys. Rev. E 107, 044402 (2023).
  * Osborne  _et al._ (2017) J. M. Osborne, A. G. Fletcher, J. M. Pitt-Francis, P. K. Maini, and D. J. Gavaghan, PLoS Comput. Biol. 13, e1005387 (2017).
  * Foty and Steinberg (2005) R. A. Foty and M. S. Steinberg, Dev. Biol. 278, 255 (2005).
  * Schmidt and Friedl (2010) S. Schmidt and P. Friedl, Cell Tissue Res 339, 83 (2010).
  * Armstrong  _et al._ (2006) N. J. Armstrong, K. J. Painter, and J. A. Sherratt, J. Theor. Biol. 243, 98 (2006).
  * Foty and Steinberg (2004) R. A. Foty and M. S. Steinberg, Int. J. Dev. Biol. 48, 397 (2004).
  * Steinberg (1963) M. S. Steinberg, Science 141, 401 (1963), 1711447 .
  * Schötz  _et al._ (2008) E.-M. Schötz, R. D. Burdine, F. Jülicher, M. S. Steinberg, C.-P. Heisenberg, and R. A. Foty, HFSP J. 2, 42 (2008).
  * Anguige and Schmeiser (2008) K. Anguige and C. Schmeiser, J. Math. Biol. 58, 395 (2008).
  * Graner and Glazier (1992) F. Graner and J. A. Glazier, Phys. Rev. Lett. 69, 2013 (1992).
  * Glazier and Graner (1993) J. Glazier and F. Graner, Phys. Rev. E 47, 2128 (1993).
  * Simpson  _et al._ (2010a) M. J. Simpson, C. Towne, D. S. McElwain, and Z. Upton, Phys. Rev. E 82, 041901 (2010a).
  * Liggett (1999) T. M. Liggett,  _Stochastic Interacting Systems: Contact, Voter and Exclusion Processes_ , Grundlehren Der Mathematischen Wissenschaften, Vol. 324 (Springer, Berlin, Heidelberg, 1999).
  * Simpson  _et al._ (2007) M. J. Simpson, A. Merrifield, K. A. Landman, and B. D. Hughes, Phys. Rev. E 76, 021918 (2007).
  * Simpson  _et al._ (2009a) M. J. Simpson, B. D. Hughes, and K. A. Landman, Australas. J. Eng. Educ. 15, 59 (2009a).
  * Simpson  _et al._ (2009b) M. J. Simpson, K. A. Landman, and B. D. Hughes, Phys. A 388, 399 (2009b).
  * Simpson  _et al._ (2010b) M. J. Simpson, K. A. Landman, and B. D. Hughes, Phys. A 389, 3779 (2010b).
  * Cai  _et al._ (2007) A. Q. Cai, K. A. Landman, and B. D. Hughes, J. Theor. Biol. 245, 576 (2007).
  * Thompson  _et al._ (2012) R. N. Thompson, C. A. Yates, and R. E. Baker, Bull. Math. Biol. 74, 2793 (2012).
  * Armstrong (1971) P. B. Armstrong, W. Roux’ Archiv f. Entwicklungsmechanik 168, 125 (1971).
  * Honda  _et al._ (1986) H. Honda, H. Yamanaka, and G. Eguchi, J Embryol Exp Morphol 98, 1 (1986).
  * Dickinson (2020) D. J. Dickinson, Curr. Biol. 30, R1071 (2020).
  * Yates  _et al._ (2017) C. A. Yates, M. J. Ford, and R. L. Mort, Bull. Math. Biol. 79, 2905 (2017).
  * Smith and Martin (1973) J. A. Smith and L. Martin, Proc. Natl. Acad. Sci. U.S.A. 70, 1263 (1973).
  * Weber  _et al._ (2014) T. S. Weber, I. Jaehnert, C. Schichor, M. Or-Guil, and J. Carneiro, PLoS Comput. Biol. 10, e1003616 (2014).
  * Mort  _et al._ (2014) R. L. Mort, M. J. Ford, A. Sakaue-Sawano, N. O. Lindstrom, A. Casadio, A. T. Douglas, M. A. Keighren, P. Hohenstein, A. Miyawaki, and I. J. Jackson, Cell Cycle 13, 2681 (2014).
  * Howard and Pelc (1951) A. Howard and S. R. Pelc, Exp. Cell Res. 2, 178 (1951).