C symbol=C, format=0 reset=section G symbol=G, format=0 reset=section M
symbol=M, format=0 reset=section B symbol=B, format=0 reset=section xi
symbol=ξ, format=0 reset=section

# Modelling non-local cell-cell adhesion: a multiscale approach

Anna Zhigun and Mabel Lizzy Rajendran School of Mathematics and Physics,
Queen’s University Belfast, University Road, Belfast BT7 1NN, Northern
Ireland, UK<EMAIL_ADDRESS>of Mathematics, Watson Building,
University of Birmingham, Edgbaston, Birmingham B15 2TT, UK,
<EMAIL_ADDRESS>
###### Abstract

Cell-cell adhesion plays a vital role in the development and maintenance of
multicellular organisms. One of its functions is regulation of cell migration,
such as occurs, e.g. during embryogenesis or in cancer. In this work, we
develop a versatile multiscale approach to modelling a moving self-adhesive
cell population that combines a careful microscopic description of a
deterministic adhesion-driven motion component with an efficient mesoscopic
representation of a stochastic velocity-jump process. This approach gives rise
to mesoscopic models in the form of kinetic transport equations featuring
multiple non-localities. Subsequent parabolic and hyperbolic scalings produce
general classes of equations with non-local adhesion and myopic diffusion, a
special case being the classical macroscopic model proposed in [4]. Our
simulations show how the combination of the two motion effects can unfold.

Cell-cell adhesion relies on the subcellular cell adhesion molecule binding.
Our approach lends itself conveniently to capturing this microscopic effect.
On the macroscale, this results in an additional non-linear integral equation
of a novel type that is coupled to the cell density equation.

Keywords: cadherin binding, cell adhesion molecule binding, cell-cell
adhesion, cell movement, diffusion-adhesion equations, hyperbolic scaling,
kinetic transport equations, multiscale modelling, myopic diffusion, non-local
models, parabolic scaling
MSC 2020: 35B27 35Q49 35Q92 45K05 92C17

##  1 Introduction

###  1.1 Biological background

Development and functioning of multicellular organisms crucially depend on
cell-cell adhesion (CCA). This is the process of cells binding to their
neighbours to form multicellular complexes by building cell-cell junctions.
Formation of new tissues and organs during embryogenesis as well as their
maintenance, be it as part of homeostasis or during wound healing, all rely on
CCA. Alteration of CCA is linked to cancer invasion and metastasis [21, 18].

Various types of cell-cell junctions exist, each enabling a specific adhesion
functionality. Mainly responsible for keeping cells together are adherens
junctions [1, Chapter 19]. They are facilitated by a particular type of cell
adhesion molecule (CAM), an adhesion-mediating transmembrane protein, called
cadherin. Cadherins require extracellular calcium ($Ca^{2+}$) in order to form
junctions. Through catenins, a family of intracellular proteins, cadherins are
indirectly connected to actin filaments that are part of the cell
cytoskeleton. Since individual cadherin bonds are rather weak, many such bonds
need to be established in parallel to secure a strong anchoring junction
between two cells. A comprehensive description of these junctions can be found
in [1, Chapter 19].

Depending on other chemical factors, cadherins either suppress migration out
of resting epithelium or support collective invasion. While E-cadherin is
responsible for strong bonds in the former case, various members of the
cadherin family that have weaker adhesion strengths than E-cadherin, such as
N-cadherin, are mainly observed when migration occurs. Loss of E-cadherin, the
main cell-cell binding CAM in epithelial cells, is believed to be a
fundamental event in the epithelial-mesenchymal transition (EMT), a process by
which cells switch from epithelial to mesenchymal stem type. In cancer, EMT
enables invasion, the precursor of metastasis. Another key adhesion type is
cell-tissue adhesion, yet in this work we concentrate solely on CCA. We refer
to [18, 21], as well as references in these sources, for further details on
the role of cell adhesion in cancer.

Motivated by the biological observations outlined above, our main aim in this
paper is to derive two prototypical classes of continuum mathematical models
for a diffusion-advection-driven motion of a self-adhesive cell population in
a heterogeneous environment. The model in [4] is regained as a special case of
one of these models, 3.14. This equation does not account for calcium-mediated
cadherin binding dynamics, whereas our novel model 4.34 includes such
dynamics.

###  1.2 Modelling background

Reaction-diffusion-advection equations with a density-dependent non-local
advection velocity in the form of a spatial integral are a popular choice when
it comes to modelling adhesion on the level of population densities. Such non-
local terms indirectly account for cell-cell interactions through the effect
that they have on the bulk motion. Starting from the integro-partial
differential equation (IPDE) that was proposed in [4] many extensions of that
adhesion model were developed and treated mathematically rigorously and
numerically, see [11] and references therein. Numerical simulations confirm
that models of this sort reproduce aggregate formation caused by CCA.

However, it is hardly possible to capture important information that needs to
be passed from lower scales to obtain a realistic model if modelling is done
directly on the level of densities. More accurate models are obtained by
zooming to the cellular or even subcellular levels and/or the level of cell
density distributions and subsequently performing an upscaling. Several
approaches to such derivations were adopted in the context of IPDE
(diffusion)-adhesion models. We briefly review them.

  1. M1

One possible starting point is a system of a large number of first-order
stochastic differential equations (SDEs). Each of the SDEs describes the
temporal evolution of the spatial position of a single population member that
interacts with other individuals and is influenced by stochastic fluctuations,
typically in the form of a Gaussian white noise. Interaction, which is often a
combination of multiple effects, is characterised by an appropriate
interaction potential. In the cases where one preliminarily considers a system
of second order SDEs for spatial positions, simplifying assumptions are made
in order to reduce to the first-order as above.

Using empirical measures and Itô’s formula, one can deduce the corresponding
mean field stochastic IPDE for population density that contains a non-locality
in the drift term. Choosing suitable scalings in each part of the interaction
potential allows to remove randomness: as the number of particles tends to
infinity, a fully deterministic macroscopic mean field IPDE is obtained. The
interaction potential contributes with non-local advection and/or the non-
linear part of diffusion. Stochastic fluctuations produce a linear diffusion
component. For a detailed discussion of this approach to derivation of non-
local diffusion-advection models spanning from modelling aspects to rigorous
mathematical treatment we refer to [32, 31] as well as references therein.

In [30], two classes of mean field IPDE models of adhesion were derived while
keeping the number of individuals finite. This was achieved by imposing
closing relations. The mean-field approximation yielded an IPDE of the same
form as in [4], whereas the Kirkwood superposition approximation produced a
non-standard system of two strongly coupled IPDEs.

No application of the described approach has dealt with CAM binding dynamics
or comparable effects.

  2. M2

Assuming that cell motion follows a space-jump random process, the evolution
of the population density can be modelled by a Master equation. The time
derivative of the density is then given by a spatial integral operator
governed by a redistribution kernel that describes the probability of jumping
from one position to another. Under suitable assumptions on the redistribution
kernel one can rescale the equation and then pass to the limit letting the
jump length tend to zero, with the result being a diffusion-advection PDE.
This approach allowed to formally derive IPDE diffusion-adhesion models in
[7], see also references therein on further details regarding derivation and
scaling of Master equations.

In [7], the redistribution kernel was split into the even and odd components,
leading to a myopic diffusion and advection, respectively. The odd component
was generated by the so-called cell’s polarization vector. This vector was
assumed to be a superposition of local adhesion strengths generated in small
test volumes in the cell’s environment. The local adhesion strength was
assumed to be proportional to: the distance to the cell, the available space,
and the amount of bound CAMs. Depending on the way the latter evolution was
modelled, one obtained adhesion velocities proportional to either single or
double spatial integrals. A special case of the first option led to the model
in [4]. The subcellular binding-unbinding dynamics of CAMs was described by
quickly equilibrating ordinary differential equations (ODEs). Their
coefficients were obtained from local densities or their integrals.

  3. M3

If it can be assumed that cell motion follows a velocity-jump process, then a
kinetic transport equation (KTE) lends itself to the description of the
evolution of the mesoscopic cell density with respect to time, position,
velocity, and, possibly, other variables (see M4 below). In the absence of
source terms, it takes the form of an IPDE where the differential transport
term captures the deterministic directed movement and a velocity integral term
cumulates the effect of stochastic fluctuations due to switching from one
velocity to another. The later term is the turning operator based on a turning
kernel that gives the probability of velocity switches. Removing the turning
operator would leave us with a mesoscopic mean field equation of the form of a
conservative transport equation (CTE). This equation can be obtained from a
microscopic ODE system describing deterministic cell movement using empirical
measures, i.e. in the same fashion as described in M1. It is commonly assumed
that velocity changes are purely stochastic in nature, implying the absence of
velocity-induced transport. However, a few works have introduced such a
transport term into the model [41, 13, 12, 10]. We mention in passing that in
the physical context similar equations exist. For instance, the linear
Boltzmann-Maxwell equation can describe a gas of charged particles moving
under the influence of an external field through an unchanging background of
another type of particles. Yet cell interactions with other cells and lifeless
matter in their surroundings are unlike collisions of physical particles with
each other or their background. Hence, different kind of forces and
interaction kernels need to be considered.

Upscaling, i.e. a suitable rescaling of a KTE using a small scaling parameter
and a limit procedure while this parameter is being sent to zero, yields a
mean field equation for the macroscopic population density. Standard scalings
are the hyperbolic and the parabolic ones. As a rule, the hyperbolic scaling
is available without preconditions, whereas the parabolic scaling is only
possible if the first velocity moment of the turning operator tends to zero.
In contrast to the parabolic upscaling where diffusion arises directly,
higher-order correction terms need to be included in order to produce
diffusion in the hyperbolic scaling case. We refer to [24] for a detailed
discussion of the two scaling types in the context of movement of living
organisms, cells in particular.

This approach was applied in the contexts of self-organised animal
aggregations [9, 16] and cell dispersal mediated by non-local sensing [28,
29]. In these works, formal derivations of advection(-diffusion) equations
were based on KTEs containing transport with respect to the spatial variable
alone. To produce non-localities on the macroscale, turning operators were
chosen that contained both velocity and spatial averaging. In the first
setting, the focus was on the interaction inside a population in one [9, 16]
and two [9] dimensions (see also references in these works with regard to
previous mesoscopic modelling in this context). This multilayered effect was
modelled by a turning kernel that is split into constants and samplings, over
the whole domain, of several interaction sub-kernels. Along with density- and
distance-dependent weightings, the latter involves functions that measure
differences between the previous velocity direction and the directions of: the
future velocity and the neighbours’ velocity and relative position. The
modelling in [28, 29] had no dimension restrictions and aimed at describing
cell polarization. There, the turning kernel samples density- and distance-
weighted measurements of a macroscopic quantity at positions along the future
velocity direction. Choosing this quantity to be the cell population density,
non-local CCA can be captured.

A suitable rescaling in [9] kept the lowest-order part of the kernel velocity-
independent, allowing for a parabolic scaling that yielded a PDE with non-
localities in both the diffusion coefficient and the adhesion velocity. Due to
the structure of the kernel, the parabolic scaling in [28, 29] was only
possible under an assumption that removed the non-locality.

CAM binding dynamics or comparable effects were not considered in either of
these works.

In the present work we develop an alternative approach to modelling CCA. Our
main goal is to construct a flexible multiscale modelling framework that
captures better the biological observations described in Subsection 1.1 above.
To be precise, we want to ensure that on the macroscale:

  1. (i)

the diffusion term originates from stochastic fluctuations;

  2. (ii)

the adhesion term is non-local,

  3. (iii)

microscopic CAM binding dynamics being its source.

Method Origin of | adhesion term | diffusion term | CAM dynamics
---|---|---|---
M1 | deterministic | deterministic & stochastic | none
M2 | stochastic & CAMs | stochastic | deterministic
M3 | stochastic | stochastic | none
Section 4 | deterministic & CAMs | stochastic | deterministic
Table 1.1: Comparison of M1-M3 and our approach.

Table 1.1 compares the above outlined approaches and our approach in Section 4
(discussed below). The comparison is based on criteria that are related to the
declared objectives (i)-(iii). Each of M1-M3 meets the requirements partially,
yet fails to meet them all. Note that even though M2 incorporates CAM binding
dynamics, that dynamics is purely deterministic (as it is generated by ODEs),
whereas the adhesion term still has a stochastic origin, as it is derived from
a redistribution kernel.

Since we aim at carefully modelling the CAM binding dynamics, we need to
recall a relevant extension of the KTE framework.

  1. M4

The kinetic theory of active particles (KTAP) [6] extends M3 to the settings
where there are non-physical ’activity’ variables that characterise the state
of a cell along with its position and velocity. In the context of cell
adhesion, this framework allowed to incorporate integrin binding dynamics.
This class of CAMs mediates cell-tissue rather than cell-cell interactions.
The corresponding models were first developed in [27, 26] to describe cancer
invasion, for which cell-tissue interactions are a prerequisite. Viewing the
proportion of bound integrins of a cell to molecules of a signal as an
activity variable, the authors derived multiscale systems that couple a KTE
for mesoscopic cell density and macroscopic reaction(-diffusion) equations for
chemical signals. This approach was taken further in [17], where a formal
upscaling was performed. Binding and unbinding of integrins was assumed to
equilibrate very quickly. For other, much slower, processes, a standard
parabolic scaling was adopted. The result was an equation containing myopic
diffusion and local advection.

For free-swimming cells moving in response to a chemical signal, similar KTEs
were constructed with active variables being certain characteristics of cell
internal state [34, 35]. Experimenting with different types of terms and
scalings led to non-standard terms on the macroscale, such as, e.g. fractional
diffusion [34] or flux-limited chemotaxis [35]. In these works, the upscaling
was done in a rigorous way.

So far, no non-local interactions of activity variables on the microscale have
been considered.

Inspired by M1-M4 as well as approaches to microscale modelling of the
deterministic portion of cell motion preceding a KTE in [13, 41] and to KTE
upscaling in [41], we use a multiscale approach to formally derive two classes
of non-local CCA models: firstly, a basic model without CAM binding dynamics
in Section 3 and, secondly, a considerably more involved model which includes
such dynamics in Section 4. Our derivations go through the following series of
steps.

  1. (1)

Develop a detailed microscopic description of the deterministic part of the
evolution of cells and, in the case of the second model, also of their CAMs;

  2. (2)

lift the modelling to the mesoscopic level of a CTE using empirical measures;

  3. (3)

introduce a turning kernel to account for stochastic fluctuations, yielding a
KTE;

  4. (4)

perform the parabolic and hyperbolic upscalings to obtain macroscopic IPDEs.

In both cases, the resulting macroscopic cell density satisfies an IPDE with
myopic diffusion and non-local adhesion. In our second model 4.34 the adhesion
strength is proportional to the fraction of bound CAMs. This quantity
satisfies, together with the cell density, a novel non-linear integral
equation.

Our strategy benefits from the accuracy and flexibility allowed by both the
microscale modelling of the deterministic motion component and the mesoscale
modelling of stochastic velocity changes. In particular, it allows to avoid
direct modelling and handling of stochasticity. They are part of method M1 and
are often challenging.

We would like to stress that this work aims at developing a modelling
framework and an understanding of what type of CCA models we should expect on
the macroscale. We do not address a specific situation that would correspond
to a concrete experiment. Another word of caution concerns the upscaling
procedures, which are only done formally. A rigorous verification, such as was
carried out in [41], is beyond the scope of the present work.

The remainder of the paper is organised as follows. After introducing some
notation in Section 2, we first derive a basic model without CAM binding
dynamics in Section 3 and then a more involved model that includes such
dynamics in Section 4. In Section 5, we present and discuss the results of
one-dimensional simulations for the model in Section 3. Finally, we discuss
and summarise our findings in Section 6.

##  2 Notation

In this Section we introduce some notations that are used throughout this
paper.

  * •

We denote by $|\cdot|$ the length of a vector but also the volume of a set in
$\mathbb{R}^{d}$, $d\in\mathbb{N}$.

  * •

For $\rho>0$ we set

$\displaystyle B_{\rho}:=$
$\displaystyle\\{\theta\in\mathbb{R}^{d}:|\theta|<\rho\\},$
$\displaystyle\overline{B}_{\rho}:=$
$\displaystyle\\{\theta\in\mathbb{R}^{d}:|\theta|\leq\rho\\},$ $\displaystyle
S_{\rho}:=$ $\displaystyle\\{\theta\in\mathbb{R}^{d}:|\theta|=\rho\\}.$

  * •

Several physical variables have their traditional meaning, i.e.
$t\in[0,\infty)$, $x\in\mathbb{R}^{d}$, and $v\in\mathbb{R}^{d}$ stand for
time, position in space, and velocity, respectively, the space dimension being
$d\in\mathbb{N}$. In the context where these and, in Section 4, yet another
variable, $y$, serve as independent variables, we refer to $t$ and $x$ as
macroscopic variables, $v$ and $y$ being referred to as non-macroscopic or
mesoscopic.

  * •

Convolution with respect to variable $x$ is denoted by $\star$.

  * •

When integrating with respect to a variable $w\in W\subset\mathbb{R}^{k}$,
$k\in\mathbb{N}$, we use the notation

$\overline{c}^{w}:=\int_{W}c\,dw$

if $c$ is a function defined on $W$ and it is evident from the context what
$W$ is. Similarly, if $c$ is a measure on $W$, we set

$\overline{c}^{w}:=\int_{W}c(dw)$

to be the integral with respect to that measure. If $w$ is a vector consisting
of all non-macroscopic variables, we omit the upper index and write
$\overline{c}$ instead. We refer to a density moment that is obtained through
such integration as macroscopic moment.

  * •

If $z_{0}$ is a point in $\mathbb{R}^{k}$, $k\in\mathbb{N}$, then
$\delta_{z_{0}}$ denotes the Dirac delta distribution centred at $z_{0}$.

  * •

We omit arguments of functions in many instances in order to simplify the
notation.

##  3 Modelling without CAM binding dynamics

In this Section we formally derive a basic non-local diffusion-adhesion model
under the assumption that adhesion originates directly from cell-cell
interactions, thus ignoring the subcellular CAM binding dynamics at this
stage. We follow steps (1)-(4) outlined in Subsection 1.2.

###  3.1 Microscale model

Similar to [41, 13, 12], we begin with a detailed description of the
deterministic part of the cell movement on the microscale which includes
acceleration due to external forces. Let a population of a large number

$1\ll N\in\mathbb{N}$

of cells be modelled as points with position and velocity coordinates

$\displaystyle(x_{i},v_{i})\in\mathbb{R}^{d}\times\mathbb{R}^{d},\qquad
i\in\\{1,\dots,N\\}.$

Following Newton’s second law, we set up an initial value problem (IVP) for an
ODE system that describes their motion:

$\displaystyle\frac{dx_{i}}{dt}=v_{i},$ (3.1a)
$\displaystyle\frac{dv_{i}}{dt}=-av_{i}+\chi(\cdot,x_{i})\frac{1}{N}\sum_{\underset{j\neq
i}{j=1}}^{N}\nabla_{x}H_{r}(x_{i}-x_{j}),$ (3.1b)
$\displaystyle(x_{i},v_{i})(0)=(x_{i0},v_{i0}),\ $ (3.1c) where
$(x_{i0},v_{i0})\in\mathbb{R}^{d}\times\mathbb{R}^{d}$, $i\in\\{1,\dots,N\\}$,
$a,r>0$,  $\displaystyle
H_{r}(x):=\frac{1}{r|B_{r}|}\int_{\min\\{|x|,r\\}}^{r}F(s)\,ds,$ (3.2)

and

$\displaystyle
F:[0,r]\rightarrow[0,\infty),\qquad\chi:[0,\infty)\times\mathbb{R}^{d}\rightarrow[0,\infty),\
\chi=\chi(t,x),$

are some continuous functions.

As in [13, 41], the term $(-av)$ on the right-hand side of 3.1b is included to
describe the acceleration (rather, deceleration) due to the viscous force.
Following Stokes’ law, we take it to be proportional to the velocity of the
cell.

The reminder of the right-hand side of 3.1b describes the acceleration due to
CCA forces. This diverts from the choices made in [41, 13, 12] where the
external forces that acted on cells were local and solely due to macroscopic
signals. It is also different from M1 because no simplifying assumptions are
made that would allow to reduce the ODE system 3.1 to a single first-order ODE
for $x_{i}$.

The adhesion force, a special case of the interaction force, is the sum of
forces due to interaction with individual cells within reach. The scaling by
$1/N$ before the sum in 3.1b corresponds to the mean field assumption that we
adopt here. Similar to [4, 30, 7] and many other works, we assume that
interaction occurs only within a sensing region that has the form of a ball of
a fixed sensing radius $r$ centred at the cell’s position. The case of a more
realistic sensing radius that could be a function of the physical variables,
as proposed in [29], is not considered here. To account for a possible spatial
heterogeneity of sensitivity to adhesion, we multiply by a parameter function
$\chi$ instead. Note that in general the resulting interaction kernel

$\displaystyle K(t,x,x^{\prime}):=\chi(t,x)\nabla_{x}H_{r}(x-x^{\prime})$
(3.3)

is not skew-symmetric with respect to the spatial variables $x$ and
$x^{\prime}$ and depends on $t$. Unlike the collision of physical particles
that follows Newton’s third law, here we allow for non-mechanical influences
on the strength of the force one cell exerts on the other.

As is standard practice, we assume the interaction force between two
individual to be proportional to the gradient of a potential, $H_{r}$, which
we refer to as adhesion potential. Function $F$ describes the dependence of
the adhesion strength upon the distance relative to $r$. The chosen domain of
integration in 3.2 ensures that no interaction occurs outside the sensing
region. The gradient of $H_{r}$ computes to

$\displaystyle\nabla_{x}H_{r}(x)=\begin{cases}-\frac{1}{r|B_{r}|}\frac{x}{|x|}F(|x|)&\text{in
}B_{r}\backslash\\{0\\},\\\ 0&\text{in
}\mathbb{R}^{d}\backslash\overline{B}_{r},\end{cases}$ (3.4)

and, unless $F(0)=F(r)=0$, it fails to exist at $0$ and/or on $S_{r}$. One
could avoid this problem by replacing $\nabla H_{r}$ in 3.1b by a function
that extends it to the whole of $\mathbb{R}^{d}$. We ignore the issue,
assuming that cells do not accumulate on lower-dimensional sets such as points
and spheres of radius $r$. Even if a small proportion of cells happens to be
at a distance exactly zero or $r$ from a certain cell at some point in time,
the corresponding contribution to the right-hand side of 3.1b would then be
small due to the factor $1/N$.

As seen in Subsection 3.3 below, the $r$-dependent coefficient $1/(r|B_{r}|)$
in 3.4 appears before the non-local adhesion term on the macroscale. We give
our motivation for its inclusion as well as argue that $F(0)$ needs to be non-
zero later in Subsection 3.3.

Finally, since cell speeds cannot become arbitrary large, we want to ensure
that they are contained in the ball ${{B}_{s}}$ for some unattainable upper
bound $s>0$. A suitable rescaling turns $s$ into $1$. Thus, from now on we
require that $v_{i}$’s do not leave the velocity space

$\displaystyle V:={B_{1}}.$

Basic ODE theory guaranties this under the condition

$\displaystyle\frac{1}{r|B_{r}|}\sup_{[0,\infty)\times\mathbb{R}^{d}}\chi\sup_{[0,r]}F\leq
a$ (3.5)

due to 3.4.

######  Remark 3.1 (Well-posedness for 3.1).

System 3.1 can be rewritten in the form

$\displaystyle\frac{dz}{dt}=g(\cdot,z),$ $\displaystyle z(0)=z_{0},$

where

$\displaystyle
z:=(x_{1}^{T},v_{1}^{T},\dots,x_{N}^{T},v_{N}^{T})^{T}:[0,\infty)\to(\mathbb{R}^{d}\times
V)^{N},$

and

$g:[0,\infty)\times(\mathbb{R}^{d}\times
V)^{N}\to[0,\infty)\times(\mathbb{R}^{d}\times V)^{N}$

is obtained by copying into a vector the right-hand sides of the equations in
3.1 in the correct order. Let us now assume that both $F$ and $\chi$ are
Lipschitz. As discussed above, $\nabla_{x}H_{r}$ may fail to exist on a lower-
dimensional set in $\mathbb{R}^{d}$, implying that the classical well-
posedness theory of first-order ODE systems cannot be used. Still, we verify
in a coming paper that $\nabla H_{r}$ belongs to the class of the vector-
valued functions of bounded variation (BV), i.e. it possesses derivatives that
are signed Radon measures. Using the chain rule for BV functions [3, Chapter 3
§3.10 Theorem 3.96]), one can deduce that $g$ inherits this property on
$(\mathbb{R}^{d}\times V)^{N}$. Furthermore, it is evident that $g$ is
essentially bounded and satisfies

$\displaystyle\nabla_{z}\cdot g(z)\equiv-dNa,$

so that the divergence is bounded. Therefore, the theory developed in [2]
provides well-posedness of 3.1 in a certain generalised sense. We do not
pursue this further here.

###  3.2 Mesoscale model

Our next step is to lift the microscopic model 3.1 to the mesoscale and extend
it to a full KTE that also accounts for stochastic velocity changes.

For each $N\in\mathbb{N}$, we introduce the time-dependent empirical measure

$\displaystyle
c_{N}(t,\cdot,\cdot):=\frac{1}{N}\sum_{i=1}^{N}\delta_{(x_{i},v_{i})(t)},\qquad
t\in[0,\infty),$ (3.6)

where $(x_{i},v_{i})$ is the trajectory that the $i$th cell follows in the
space-velocity space. This measure-valued function is an appropriate
description of the mesoscopic population density rescaled so that the total
mass is normalized to one. Each distribution $\delta_{(x_{i},v_{i})}$ models a
point mass concentrated at $(x_{i},v_{i})$, i.e. the density of a cell at
$x_{i}$ with velocity $v_{i}$.

Let us assume for a moment that $H_{r}$ is sufficiently regular. In this case,
the classical ODE theory provides the well-posedness of 3.1. Moreover, the
empirical measure $c_{N}$ corresponding to the solution of 3.1 satisfies in
the distributional sense the mean field PDE

$\displaystyle\nabla_{(t,x,v)}\cdot\left(\left(1,v,-av+\chi\nabla_{x}H_{r}\star\overline{c_{N}}\right)c_{N}\right)=0$
(3.7)

and, obviously, also the initial condition

$\displaystyle
c_{N}(0,\cdot,\cdot)=\frac{1}{N}\sum_{j=1}^{N}\delta_{(x_{i0},v_{i0})}.$

For constant $\chi$, this is well-known, see, e.g. [20]. The general case
follows with Lemma A.1 in the Appendix. However, if $H_{r}$ is not regular
enough, 3.7 may fail to make sense. In particular,
$(\nabla_{x}H_{r}\star\bar{c})c$ is in general not well-defined if
$\nabla_{x}H_{r}$ is not continuous and $c$ is a discrete, hence singular,
measure. For the reason detailed in Subsection 3.1 we ignore this issue here.

Since the population number is assumed to be large, we are interested in the
mean field limit as $N\to\infty$. This allows to deal with less concentrated,
hence less singular, solutions to 3.7 that are functions and not discrete
measures. Since 3.7 does not depend on $N$, it is reasonable to expect that
this is the equation that is obtained in the limit, i.e. that $c_{N}$
converges to some $c$ that satisfies

$\displaystyle\nabla_{(t,x,v)}\cdot\left(\left(1,v,-av+\chi\nabla_{x}H_{r}\star\overline{c}\right)c\right)=0.$
(3.8)

The CTE 3.8 provides the description of the deterministic part of cell
movement driven by 3.1 on the mesoscopic level. To complete the modelling, we
still need to include a term that accounts for stochastic perturbations. Since
adhesion is particularly relevant in cancer (see Subsection 1.1), we include a
turning operator that accounts for chaotic realignment with tissue fibers.
Following [41, 24], we choose a basic turning operator

$\displaystyle c\mapsto dq\overline{c}-c$

to illustrate our approach. Here $q$ models the orientational distribution of
tissue fibers and satisfies the following assumptions:

######  Assumptions 3.2.

  1. 1.

$q:\mathbb{R}^{d}\times V\to[0,\infty)$ and only depends on $x$ and
$\frac{v}{|v|}$;

  2. 2.

$\bar{q}=\frac{1}{d}$.

This kind of turning operator has been used in many models for cancer
migration, see, e.g. references given in [41]. The resulting mesoscopic
equation is

$\displaystyle\nabla_{(t,x,v)}\cdot\left(\left(1,v,-av+\chi\nabla_{x}H_{r}\star\overline{c}\right)c\right)=dq\overline{c}-c.$
(3.9)

It is a blend of a KTE and a mean field equation. This doubly non-local IPDE
accounts for both the deterministic cell-cell and stochastic cell-tissue
interactions.

Due to 3.5 and 3.4 we have that

$\displaystyle(-av+\chi\nabla_{x}H_{r}\star u)\cdot v\leq$
$\displaystyle-a|v|^{2}+|v|\sup_{[0,\infty)\times\mathbb{R}^{d}}\chi\,\underset{B_{r}}{\sup}|\nabla
H_{r}|\|u\|_{L^{1}(\mathbb{R}^{d})}$ $\displaystyle=$
$\displaystyle-a+\sup_{[0,\infty)\times\mathbb{R}^{d}}\chi\,\underset{B_{r}}{\sup}|\nabla
H_{r}|$ $\displaystyle\leq$ $\displaystyle 0\qquad\text{in
}[0,\infty)\times\mathbb{R}^{d}\times S_{1}\qquad\text{for
}\|u\|_{L^{1}(\mathbb{R}^{d})}=1.$ (3.10)

Consequently, the characteristics of the transport part of equation 3.9 that
start in $\mathbb{R}^{d}\times V$ do not leave this set. Hence,

$\displaystyle c=0\qquad\text{in }[0,\infty)\times\mathbb{R}^{d}\times S_{1}$
(3.11)

are admissible boundary conditions for equation 3.9

######  Remark 3.3 (Mean field limit for 3.7).

For constant $\chi$, the fact that 3.8 is obtained from 3.7 in the mean field
limit is a direct consequence of the results in [14], provided that $H_{r}$ is
smooth, and [25] if it is not, at least for $a=0$.

######  Remark 3.4 (Solvability of 3.9).

We are not aware of results on solvability for such doubly non-local non-
linear equations as 3.9.

###  3.3 Macroscale model

In this Subsection we upscale 3.9, 3.11 to obtain equations for the
macroscopic cell density. To begin with, we introduce a macroscopic rescaling
of time and space and of functions depending on them: for
$\varepsilon\in(0,1]$ let

$\displaystyle\hat{t}=\varepsilon^{\kappa}t,\qquad\kappa\in\\{1,2\\},$
$\displaystyle\hat{x}=\varepsilon x,\qquad\hat{r}=\varepsilon r,$
$\displaystyle\hat{\chi}(\hat{t},\hat{x}):=\chi(t,x),\qquad\hat{F}(\hat{s}):=F(s),\qquad\hat{q}(\hat{x},v):=q(x,v),$
$\displaystyle c^{\varepsilon}(\hat{t},\hat{x},v):=c(t,x,v).$

The values $\kappa=1$ and $\kappa=2$ correspond to the usual hyperbolic and
parabolic scalings, respectively. Under the proposed scaling we have

$\displaystyle\nabla_{x}H_{r}\star\overline{c}(t,x)=$
$\displaystyle-\frac{1}{r|B_{r}|}\int_{B_{r}}\frac{y}{|y|}F(|y|)\overline{c}(t,x-y)\,dy$
$\displaystyle=$
$\displaystyle-\frac{\varepsilon^{d+1}}{\hat{r}|B_{\hat{r}}|}\int_{B_{\frac{\hat{r}}{\varepsilon}}}\frac{y}{|y|}\hat{F}(\varepsilon|y|)\overline{c^{\varepsilon}}\left(\hat{t},\hat{x}-\varepsilon
y\right)\,dy$ $\displaystyle=$
$\displaystyle-\frac{\varepsilon}{\hat{r}|B_{\hat{r}}|}\int_{B_{\hat{r}}}\frac{y}{|y|}F(\hat{y})\overline{c^{\varepsilon}}\left(\hat{t},\hat{x}-\hat{y}\right)\,d\hat{y}$
$\displaystyle=$
$\displaystyle\varepsilon\nabla_{\hat{x}}H_{\hat{r}}\star\overline{c^{\varepsilon}}(\hat{t},\hat{x}).$
(3.12)

Rescaling 3.9 and 3.11, using 3.12, and dropping the hats leads to

$\displaystyle\nabla_{(t,x,v)}\cdot\left(\left(\varepsilon^{\kappa},\varepsilon
v,-a(v-v^{\varepsilon}_{*})\right)c^{\varepsilon}\right)=dq\overline{c^{\varepsilon}}-c^{\varepsilon},$
(3.13a) $\displaystyle c^{\varepsilon}=0\qquad\text{in }\mathbb{R}^{d}\times
S_{1},$ (3.13b) $\displaystyle
v^{\varepsilon}_{*}:=\varepsilon\frac{1}{a}\chi\nabla_{x}H_{r}\star\overline{c^{\varepsilon}}.$
(3.13c)

Set

$\displaystyle c^{0}:=\underset{\varepsilon\rightarrow
0}{\lim}\,c^{\varepsilon},$ $\displaystyle
c^{0}_{1}:={\underset{\varepsilon\rightarrow
0}{\lim}\,\partial_{\varepsilon}c^{\varepsilon}},$ $\displaystyle
c_{01}^{\varepsilon}:=c^{0}+\varepsilon c^{0}_{1}.$

We are interested in obtaining equations for the macroscopic zero- and first-
order approximations of $c^{\varepsilon}$, i.e. $\overline{c^{0}}$ and
$\overline{c_{01}^{\varepsilon}}$. Equation 3.13a has the same form as
equation (3.3) in [41]. However, the term $v^{\varepsilon}_{*}$ is not exactly
of the form we considered in [41]. Indeed, it depends on variable $t$ and
lacks saturation. Still, since it vanishes at $\varepsilon=0$, the very same
formal derivation as was done in that work can be carried out in the present
case. In particular, one obtains equations

$\displaystyle(a+1)\partial_{t}\overline{c^{0}}=\frac{1}{2a+1}\frac{d}{d+2}\nabla_{x}\nabla_{x}^{T}:\left(\mathbb{D}[q]\overline{c^{0}}\right)-\nabla_{x}\cdot(\overline{c^{0}}\chi\nabla_{x}H_{r}\star\overline{c^{0}})\qquad\text{if
}\kappa=2\text{ and }\mathbb{E}[q]\equiv 0$ (3.14)

and

$\displaystyle(a+1)\partial_{t}\overline{c_{01}^{\varepsilon}}+\frac{d}{d+1}\nabla_{x}\cdot\left(\mathbb{E}[q]\overline{c_{01}^{\varepsilon}}\right)$
$\displaystyle=$
$\displaystyle{\varepsilon}\left(\frac{1}{2a+1}\frac{d}{d+2}\nabla_{x}\nabla_{x}^{T}:\left(\mathbb{D}[q]\overline{c_{01}^{\varepsilon}}\right)-\frac{1}{(a+1)^{2}}\frac{d^{2}}{(d+1)^{2}}\nabla_{x}\cdot\left(\mathbb{E}[q]\nabla_{x}\cdot\left(\overline{c_{01}^{\varepsilon}}\mathbb{E}[q]\right)\right)\right)$
$\displaystyle-\varepsilon\nabla_{x}\cdot\left(\overline{c_{01}^{\varepsilon}}\chi\nabla_{x}H_{r}\star\overline{c_{01}^{\varepsilon}}\right)$
$\displaystyle+O\left(\varepsilon^{2}\right)\qquad\text{if }\kappa=1,$ (3.15)

where

$\displaystyle\mathbb{E}[q]:=\int_{S_{1}}\theta q(\theta)\,d\theta,$
$\displaystyle\mathbb{D}[q]:=\int_{S_{1}}\theta\theta^{T}q(\theta)\,d\theta.$

Recalling the adhesion operator

$\displaystyle{\cal A}_{r}u(x)=$
$\displaystyle\frac{1}{r|B_{r}|}\int_{B_{r}}u(x+\xi)\frac{\xi}{|\xi|}F(|\xi|)\,d\xi$

from [4] and noticing that

$\displaystyle{\cal A}_{r}u\equiv$ $\displaystyle\nabla_{x}H_{r}\star u,$
(3.17)

we can alternatively rewrite 3.14 and 3.15 as follows:

$\displaystyle(a+1)\partial_{t}\overline{c^{0}}=\frac{1}{2a+1}\frac{d}{d+2}\nabla_{x}\nabla_{x}^{T}:\left(\mathbb{D}[q]\overline{c^{0}}\right)-\nabla_{x}\cdot(\overline{c^{0}}\chi{\cal
A}_{r}\overline{c^{0}})\qquad\text{if }\kappa=2\text{ and }\mathbb{E}[q]\equiv
0$ (3.18)

and

$\displaystyle(a+1)\partial_{t}\overline{c_{01}^{\varepsilon}}+\frac{d}{d+1}\nabla_{x}\cdot\left(\mathbb{E}[q]\overline{c_{01}^{\varepsilon}}\right)$
$\displaystyle=$
$\displaystyle{\varepsilon}\left(\frac{1}{2a+1}\frac{d}{d+2}\nabla_{x}\nabla_{x}^{T}:\left(\mathbb{D}[q]\overline{c_{01}^{\varepsilon}}\right)-\frac{1}{(a+1)^{2}}\frac{d^{2}}{(d+1)^{2}}\nabla_{x}\cdot\left(\mathbb{E}[q]\nabla_{x}\cdot\left(\overline{c_{01}^{\varepsilon}}\mathbb{E}[q]\right)\right)\right)$
$\displaystyle-\varepsilon\nabla_{x}\cdot\left(\overline{c_{01}^{\varepsilon}}\chi{\cal
A}_{r}\overline{c_{01}^{\varepsilon}}\right)$
$\displaystyle+O\left(\varepsilon^{2}\right)\qquad\text{if }\kappa=1.$ (3.19)

Both 3.18 and 3.19 contain the same $q$-dependent terms, such as, e.g. the
myopic diffusion

$\displaystyle\nabla_{x}\nabla_{x}^{T}:\left(\mathbb{D}[q]u\right),$

as the corresponding equations (3.18) and (3.51) from [41]. We refer to that
work for a discussion of these terms, as well as for the formulas for the
mesoscopic approximations $c^{0}$ and $c_{1}^{0}$. In contrast to [41], our
new equations 3.18 and 3.19 also contain the non-local advection term

$\displaystyle-\nabla_{x}\cdot(u\chi{\cal A}_{r}u)$

of the form originally proposed in [4] to model CCA. That model is a special
case of 3.18 and corresponds to $\chi$ and $q$ being constant.

Another way to rewrite, e.g. 3.18 is by decomposing the spatial motion into an
anisotropic diffusion in the divergence form and advection:

$\displaystyle(a+1)\partial_{t}\overline{c^{0}}=\frac{1}{2a+1}\frac{d}{d+2}\nabla_{x}\cdot\left(\mathbb{D}[q]\nabla_{x}\overline{c^{0}}\right)+\nabla_{x}\cdot\left(\overline{c^{0}}\left(\frac{1}{2a+1}\frac{d}{d+2}\nabla_{x}\cdot\mathbb{D}[q]-\chi{\cal
A}_{r}\overline{c^{0}}\right)\right).$ (3.20)

In Section 5 below we show results of some numerical simulations for this
model.

The heuristic analysis in [19] as well as the rigorous study in [15] showed
that in the limit as $r\to 0$ the non-local adhesion operator ${\cal A}_{r}$
approaches the (local) spatial gradient, provided that $F(0)=d+1$. This is the
expected limit behaviour and the reason for the factor $1/(r|B_{r}|)$ in
Subsection 3.3 and the observation that $F(0)\neq 0$, see Subsection 3.1.

######  Remark 3.5 (Non-local operator $\mathring{\nabla}_{r}$).

Since our approach to the derivation of the mesoscopic equation 3.9 from the
microscopic ODE system 3.1 is based on empirical measures, it limits the
admissible choices of the interaction potential. If we were to start directly
on the mesoscopic level and would only accept integrable densities $c$ rather
than discrete measures, then we could choose a discontinuous potential such as

$\displaystyle H_{r}(x):=\begin{cases}1&\text{in }B_{r},\\\ 0&\text{in
}\mathbb{R}^{d}\backslash\overline{B}_{r}.\end{cases}$ (3.21)

In this case, the gradient of $H_{r}$ is a measure. In higher dimensions
$d\geq 2$ it is given by

$\displaystyle\nabla_{x}H_{r}=-id_{x}\,dS_{r},$

which leads to

$\displaystyle\nabla_{x}H_{r}\star u=$
$\displaystyle\int_{S_{r}}u(x+\xi)\xi\,d\xi$ $\displaystyle=$
$\displaystyle\mathring{\nabla}_{r}u.$ (3.22)

The latter is the non-local operator that was previously introduced in [33] to
describe non-local chemotaxis.

######  Remark 3.6 (Solvability of 3.18 and 3.19).

Several works established solvability of non-local diffusion-adhesion
equations as well as systems containing them, see [11] and references therein.
Yet none of them included the case of a myopic diffusion. Solvability in the
presence of a scalar [40] or a tensor [22] myopic diffusion has so far been
accomplished for models with advection due to haptotaxis, i.e. a directed
movement along the spatial gradient of the macroscopic tissue density, rather
than adhesion.

######  Remark 3.7 (Rigorous upscaling of 3.9).

The presented meso-to-macro upscaling is formal. Unlike the case that was
handled in [41], 3.9 is a non-linear equation and, as previously observed in
Subsection 3.2, the term $u{\cal A}_{r}u$ is not defined for singular measures
$u$. Thus, the approach that we developed in [41] is not directly applicable.
Indeed, there we relied on the linearity of the KTE and on the possibility of
considering measure-valued solutions. A rigorous upscaling for 3.9 remains an
open question.

##  4 Modelling with CAM binding dynamics

In this Section we derive a new non-local diffusion-adhesion model that takes
into account subcellular CAM binding dynamics. We adhere to steps (1)-(4)
outlined in Subsection 1.2.

###  4.1 Microscale model

The basic microscopic model 3.1 neglects the CAMs binding dynamics, which, as
pointed out in the Introduction, is the underlying mechanism of cell-cell
binding. In [7], this mechanism was taken into account. There it was assumed
that at each time and position in space a single cell is moving while the rest
of the population in its background is effectively standing still. The
interactions between the cell and the background population were described by
reversible reactions that either discriminate between bound/free CAMs or not.

In this Subsection, we exploit the microscopic approach that allows to
describe mutual interactions between the CAMs of a pair of cells. For this
purpose, we construct a system of ODEs that includes equations not just for
$x_{i}$ and $v_{i}$ but also for $y_{i}$, the proportion of bound CAMs of
$i$th cell. This bears resemblance to modelling in [27, 26], although, in our
case, the interactions occur inside the population rather than with an
external signal.

For the reader’s convenience, we first collect all involved model variables
and parameters, including those previously introduced in Section 3:

  * •

$1\ll N\in\mathbb{N}$: population number;

  * •

$t\in[0,\infty)$: time;

  * •

$d\in\mathbb{N}$: space dimension;

  * •

$x_{i}:[0,\infty)\rightarrow\mathbb{R}^{d}$, $x_{i}=x_{i}(t)$: position of
$i$th cell in space;

  * •

$V=B_{1}$: velocity space;

  * •

$v_{i}:[0,\infty)\rightarrow V$, $v_{i}=v_{i}(t)$: velocity of $i$th cell;

  * •

$a\in(0,\infty)$: cell deceleration rate;

  * •

$r\in(0,\infty)$: cell sensing radius;

  * •

$R\in(0,\infty)$: total number of CAMs of a cell;

  * •

$y_{i}:[0,\infty)\rightarrow(0,1)$, $y_{i}=y_{i}(t)$: proportion of bound CAMs
of $i$th cell;

  * •

$S:[0,\infty)\times\mathbb{R}^{d}\rightarrow[0,\infty)$, $S=S(t,x)$:
concentration of a chemical on which the likeliness to bind/unbind depends;

  * •

$k^{+}/k^{-}:[0,\infty)\times[0,\infty)\rightarrow[0,\infty)$ : CAM
binding/unbinding rate constants, depend on $S$ and the distance between
interacting cells;

  * •

$F:[0,r]\rightarrow[0,\infty)$: distance-dependent component of adhesion
force;

  * •

$H_{r}$: adhesion potential as defined by 3.2;

  * •

$\chi:[0,\infty)\times\mathbb{R}^{d}\rightarrow[0,\infty)$, $\chi=\chi(t,x)$:
adhesion sensitivity.

We begin by describing the cell motion: for $i\in\\{1,\dots,N\\}$

$\displaystyle\frac{dx_{i}}{dt}=$ $\displaystyle v_{i},$ (4.1a)
$\displaystyle\frac{dv_{i}}{dt}=$ $\displaystyle-
av_{i}+\chi(\cdot,x_{i})y_{i}\frac{1}{N}\sum_{\underset{j\neq
i}{j=1}}^{N}\nabla_{x}H_{r}(x_{i}-x_{j}).$ (4.1b)

Unlike 3.1, in 4.1 the total adhesion force acting on the $i$th cell is taken
to be proportional to $y_{i}$. This implies that the more bonds a cell has,
the stronger this force is.

We make the following assumptions on the CAM binding:

  1. (Ai)

a single type of CAMs influences cell motion;

  2. (Aii)

all cells have exactly the same number of these CAMs;

  3. (Aiii)

each cell has both bound and free CAMs at every point in time;

  4. (Aiv)

a CAM of a cell can only bind to a distinct CAM of another cell;

  5. (Av)

CAMs of a pair of cells can bind only if the distance between the cells is
smaller than the sensing radius $r$;

  6. (Avi)

if the distance approaches $r$, then their CAMs unbind;

  7. (Avii)

binding/unbinding kinetics obeys an analog of the mass action law;

  8. (Aviii)

binding is reversible; the corresponding rates for a pair of cells depends on
the distance between them and the concentration of a chemical at the middle
distance.

An example of an application that we have in mind here is a simplified
description of the formation of adherens junctions through the calcium-
mediated cadherin binding, see Subsection 1.2 above.

Let ${\cal F}_{i}$ and ${\cal B}_{i}$ denote respectively the free and bound
CAMs of the $i$th cell. Then, the above assumptions can be described by the
following ‘reactions’: for all $i,j\in\\{1,\dots,N\\}$, $i\neq j$,

$\displaystyle\ce{{\cal F}_{i}+{\cal
F}_{j}<=>[$k^{+}\left(S\left(t,\frac{1}{2}(x_{i}+x_{j})\right),|x_{i}-x_{j}|\right)$][$k^{-}\left(S\left(t,\frac{1}{2}(x_{i}+x_{j})\right),|x_{i}-x_{j}|\right)$]{\cal
B}_{i}+{\cal B}_{j}},$ (4.2)

where

$\displaystyle k^{\pm}\equiv 0\qquad\text{in }[0,\infty)\times[r,\infty).$
(4.3)

It is reasonable to suppose in this context that the following quantities
replace the standard chemical concentrations:

$\displaystyle[{\cal B}_{i}]:=$ $\displaystyle\frac{[\text{number of bound
CAMs of }i\text{th cell}]}{[\text{total number of CAMs in
population}]\cdot[\text{volume of the sensing region}]}$ $\displaystyle=$
$\displaystyle\frac{y_{i}R}{NR|B_{r}|}$ $\displaystyle=$
$\displaystyle\frac{y_{i}}{N|B_{r}|}$

and, similarly,

$\displaystyle[{\cal F}_{i}]:=\frac{1-y_{i}}{N|B_{r}|}.$

We provide an example of binding/unbinding rates $k^{+}$/$k^{-}$.

######  Example 4.1.

For some non-decreasing $K^{+}:[0,\infty)\to[0,\infty)$ and non-increasing
$K^{-}:[0,\infty)\to[0,\infty)$ and constants $a^{\pm},b^{\pm}\in(0,\infty)$
set

$\displaystyle k^{\pm}(S,\rho):=K^{\pm}(S)\varphi^{\pm}(\rho)\qquad\text{for
}S,\rho\in[0,\infty),$
$\displaystyle\varphi^{\pm}(\rho):=\begin{cases}\left(r^{b^{\pm}}-\rho^{b^{\pm}}\right)^{\pm
a^{\pm}}&\text{for }\rho\in[0,r),\\\ 0&\text{for
}\rho\in[r,\infty).\end{cases}$

This choice produces a speedy detachment when the cell distance is close to
$r$.

In order to estimate the binding/unbinding rates in a concrete type of cell-
cell binding, the general framework developed in [5] can be adopted. There,
the binding process is decomposed into two phases: the formation of the
encounter complex and the actual bond formation, both being reversible
reactions. The reaction rates for the former are described by functions of the
cell separating distance and the translational diffusion coefficient for CAM
motion in the cell membrane. Specifically in the case of cadherin binding,
higher calcium concentrations correlate with faster diffusion of cadherins in
the cell membrane [8]. As to the bond formation, it is well-understood [1,
Chapter 19] that calcium is indispensable for cadherins to achieve the rigid
structure that is necessary for them to bind. These are the observations that
have led us to assume that $k^{\pm}$ are functions of distance and
concentration of a chemical that mediates the binding of CAMs.

Applying the law of mass action to 4.2, we arrive at the following ODE system
for the dynamics of $y_{i}$’s:

$\displaystyle\frac{dy_{i}}{dt}=$
$\displaystyle\frac{1}{N}\sum_{\underset{j\neq
i}{j=1}}^{N}G_{r}[S](t,(x_{i},y_{i}),(x_{j},y_{j})),\qquad
i\in\\{1,\dots,N\\},$ (4.4a) where $\displaystyle
G_{r}[S](t,(x,y),(x^{\prime},y^{\prime})):=G_{r}^{+}[S](t,x,x^{\prime})(1-y)(1-y^{\prime})-G_{r}^{-}[S](t,x,x^{\prime})yy^{\prime},$
(4.4b) $\displaystyle
G_{r}^{{\pm}}[S](t,x,x^{\prime}):=\frac{1}{|B_{r}|}k^{{\pm}}\left(S\left(t,\frac{1}{2}(x+x^{\prime})\right),|x-x^{\prime}|\right).$
(4.4c)

System 4.1, 4.4 is our new microscopic model.

Since $k^{\pm}\geq 0$, we have that

$\displaystyle G_{r}[S](t,(x,0),(x^{\prime},y^{\prime}))\geq 0,\
G_{r}[S](t,(x,1),(x^{\prime},y^{\prime}))\leq 0\qquad\text{for all
}t\in[0,\infty),\ x,x^{\prime}\in\mathbb{R}^{d},\ y^{\prime}\in[0,1].$ (4.5)

Standard ODE theory implies that $y_{i}$’s do not leave $[0,1]$.

######  Remark 4.2 (Well-posedness of 4.1, 4.4).

The well-posedness of 4.1, 4.4 can be addressed in the same way as for 3.1,
see Remark 3.1. Depending on the choice of functions $k^{\pm}$ one could avail
of the classical ODE existence and uniqueness results or the more general
theory from [2]. In particular, since the divergence of the right-hand side of
the ODE system for
$(x_{1}^{T},v_{1}^{T},y_{1},\dots,x_{N}^{T},v_{N}^{T},y_{N})^{T}$ computes to

$\displaystyle-dNa-\frac{1}{N}\sum_{i=1}^{N}\sum_{\underset{j\neq
i}{j=1}}^{N}(G_{r}^{+}[S](t,x_{i},x_{j})(1-y_{j})+G_{r}^{-}[S](t,x_{i},x_{j})y_{j})$
$\displaystyle=$ $\displaystyle-
dNa-\frac{1}{N}\sum_{j=1}^{N}\left((1-y_{j})\sum_{\underset{j\neq
i}{i=1}}^{N}G_{r}^{+}[S](t,x_{i},x_{j})+y_{j}\sum_{\underset{j\neq
i}{i=1}}^{N}G_{r}^{-}[S](t,x_{i},x_{j})\right),$ (4.6)

the essential boundedness of $k^{\pm}$ is necessary for the results from [2]
to apply.

###  4.2 Mesoscale model

In this Subsection, we lift the microscopic model 4.1, 4.4 to the mesoscale
and extend it to a full KTE that includes stochastic velocity changes. Similar
to Subsection 3.2, we begin by introducing the empirical measures

$\displaystyle
c_{N}(t,\cdot,\cdot,\cdot):=\frac{1}{N}\sum_{j=1}^{N}\delta_{(x_{i},v_{i},y_{i})(t)}.$

A formal application of Lemma A.1 to system 4.1, 4.4 and
$Z:=(x^{T},v^{T},y)^{T}$ yields that $c_{N}$ solves in the distributional
sense the mean field IVP

$\displaystyle\nabla_{(t,x,v,y)}\cdot\left(\left(1,v,-av+\chi
y\nabla_{x}H_{r}\star\overline{c_{N}},{\cal
G}_{r}[S]\overline{c_{N}}^{v}-\frac{1}{N}\operatorname{tr}G_{r}[S]\right)c_{N}\right)=0,$
(4.7a) $\displaystyle
c_{N}(0,\cdot,\cdot,\cdot)=\frac{1}{N}\sum_{j=1}^{N}\delta_{(x_{i},v_{i},y_{i})(0)}.$

where

$\displaystyle{\cal
G}_{r}[S]u\,(t,x,y):=\int_{0}^{1}\int_{\mathbb{R}^{d}}G_{r}[S](t,(x,y),(x^{\prime},y^{\prime}))u(t,x^{\prime},y^{\prime})\,dx^{\prime}dy^{\prime},$
$\displaystyle\operatorname{tr}G_{r}[S](t,x,y):=G_{r}[S](t,(x,y),(x,y)).$

Here, as in Subsection 3.2, we ignore the potential discontinuities in the
kernels. Passing formally to the limit as $N\rightarrow\infty$ in 4.7a, we
arrive at the mean field limit equation

$\displaystyle\nabla_{(t,x,v,y)}\cdot\left(\left(1,v,-av+\chi
y\nabla_{x}H_{r}\star\overline{c},{\cal
G}_{r}[S]\overline{c}^{v}\right)c\right)=0.$ (4.8)

To account for chaotic interactions with tissue we use the same turning
operator as in 3.9. The resulting mesoscopic equation is thus:

$\displaystyle\nabla_{(t,x,v,y)}\cdot\left(\left(1,v,-av+\chi
y\nabla_{x}H_{r}\star\overline{c},{\cal
G}_{r}[S]\overline{c}^{v}\right)c\right)=dq\overline{c}^{v}-c.$ (4.9)

Due to 3.5 and 3.4 we have that

$\displaystyle(-av+\chi y\nabla_{x}H_{r}\star u)\cdot v\leq$
$\displaystyle-a|v|^{2}+|v|\underset{B_{r}}{\sup}|\nabla
H_{r}|\|u\|_{L^{1}(\mathbb{R}^{d})}$ $\displaystyle=$
$\displaystyle-a+\underset{B_{r}}{\sup}|\nabla H_{r}|$ $\displaystyle\leq$
$\displaystyle 0\qquad\text{in }[0,\infty)\times\mathbb{R}^{d}\times
S_{1}\times[0,1]\qquad\text{for }\|u\|_{L^{1}(\mathbb{R}^{d})}=1.$ (4.10)

Further, 4.5 implies that

$\displaystyle{\cal G}_{r}[S]u\geq 0\qquad\text{in
}[0,\infty)\times\mathbb{R}^{d}\times\\{0\\},\qquad\text{ for }u\geq 0,$
(4.11a) $\displaystyle{\cal G}_{r}[S]u\leq 0\qquad\text{in
}[0,\infty)\times\mathbb{R}^{d}\times\\{1\\},\qquad\text{ for }u\geq 0.$
(4.11b)

Combining 4.10 and 4.11, we conclude that the characteristics of the transport
part of equation 4.8 that start in $\mathbb{R}^{d}\times V\times[0,1]$ do not
leave this set. Hence,

$\displaystyle c=0\qquad\text{in
}[0,\infty)\times\mathbb{R}^{d}\times((\overline{B_{1}}\times\\{0,1\\})\cup(S_{1}\times(0,1)))$
(4.12)

are admissible boundary conditions for 4.9.

######  Remark 4.3 (Rigorous treatment of 4.7a).

We are not aware of rigorous results on well-posedness or mean field limit for
such CTEs as 4.7a. Note that unlike standard applications arising in physics,
the kernels of the integral operators we are dealing with here are not skew
symmetric.

######  Remark 4.4 (Solvability of 4.9).

Equation 4.9 is a generalisation of 3.9. As mentioned in Remark 3.4, the
solvability of the latter equation has not been addressed so far.

###  4.3 Macroscale model

In this Subsection we upscale 4.9, 4.12 to obtain equations for the
macroscopic cell density. As in Subsection 3.2, we begin by introducing a
macroscopic rescaling of time and space and of functions depending on them:
for $\varepsilon\in(0,1]$ let

$\displaystyle\hat{t}=\varepsilon^{\kappa}t,\qquad\text{for
}\kappa\in\\{1,2\\},$ $\displaystyle\hat{x}=\varepsilon
x,\qquad\hat{r}=\varepsilon r,\qquad\hat{d}=\varepsilon d,$
$\displaystyle\hat{S}(\hat{t},\hat{x}):=S(t,x),\qquad\hat{\chi}(\hat{t},\hat{x}):=\chi(t,x),\qquad\hat{F}(\hat{s}):=F(s),\qquad\hat{q}(\hat{x},v):=q(x,v),$
$\displaystyle\hat{k}^{\pm}(S,\hat{d}):=\varepsilon^{-\mu}k^{\pm}(S,d)\qquad\text{for
some }\mu>0,$ $\displaystyle
c^{\varepsilon}(\hat{t},\hat{x},v,y):=c(t,x,v,y).$

As before, we consider here hyperbolic ($\kappa=1$) and parabolic ($\kappa=2$)
space-time scalings. The rescaling of reaction rate constants means rescaling
of $dy/dt$. A negative epsilon power is chosen to reflect the fact that the
CAM binding and unbinding are the fastest among all included processes.
Rescaling 4.9 and 4.12 and dropping the hats leads to

$\displaystyle\nabla_{(t,x,v,y)}\cdot\left(\left(\varepsilon^{\kappa},\varepsilon
v,-av+\varepsilon\chi
y\nabla_{x}H_{r}\star\overline{c^{\varepsilon}},\varepsilon^{-\mu}{\cal
G}_{r}[S]\overline{c^{\varepsilon}}^{v}\right)c^{\varepsilon}\right)=dq\overline{c^{\varepsilon}}^{v}-c^{\varepsilon},$
(4.13a) $\displaystyle c^{\varepsilon}=0\qquad\text{in
}{[0,\infty)\times\mathbb{R}^{d}\times((\overline{B_{1}}\times\\{0,1\\})\cup(S_{1}\times(0,1)))}.$
(4.13b)

Following the approach from [41], we derive equations connecting some zero,
first, and second moments of $c^{\varepsilon}$. To begin with, we integrate
4.13a by parts with respect to $(v,y)$ over $V\times[0,1]$ using 4.13b and
then divide by $\varepsilon^{\kappa}$ in order to obtain an equation which
connects the macroscopic zero and first order $v$ moments:

$\displaystyle\partial_{t}\overline{c^{\varepsilon}}+\varepsilon^{1-\kappa}\nabla_{x}\cdot\overline{vc^{\varepsilon}}=0.$
(4.14)

Next, we multiply 4.13a by $v$ and once again integrate by parts over
$V\times[0,1]$ using 4.13b:

$\displaystyle\varepsilon^{\kappa}\partial_{t}\overline{vc^{\varepsilon}}+\varepsilon\nabla_{x}\cdot\overline{vv^{T}c^{\varepsilon}}+a\overline{vc^{\varepsilon}}-\varepsilon\chi\overline{yc^{\varepsilon}}\nabla_{x}H_{r}\star\overline{c^{\varepsilon}}=\frac{n}{n+1}\mathbb{E}[q]\overline{c^{\varepsilon}}-\overline{vc^{\varepsilon}}.$
(4.15)

Rearranging and dividing 4.15 by $\varepsilon^{\kappa-1}$ leads to

$\displaystyle-(a+1)\varepsilon^{1-\kappa}\overline{vc^{\varepsilon}}=$
$\displaystyle\varepsilon^{2-\kappa}\nabla_{x}\cdot\overline{vv^{T}c^{\varepsilon}}-\varepsilon^{2-\kappa}\left(\chi\overline{yc^{\varepsilon}}\nabla_{x}H_{r}\star\overline{c^{\varepsilon}}+\varepsilon^{-1}\frac{n}{n+1}\mathbb{E}[q]\overline{c^{\varepsilon}}\right)+\varepsilon\partial_{t}\overline{vc^{\varepsilon}}.$
(4.16)

Next, we apply $(\nabla_{x}\cdot)$ to both sides of 4.16 and plug the
expression on the right-hand side into 4.14. In order to eliminate the
resulting term with the mixed derivative $(\nabla_{x}\cdot)\partial_{t}$ we
apply $\varepsilon^{\kappa}\partial_{t}$ to both sides of 4.14. Thus we arrive
at the following differential equation for the macroscopic moments of zero and
second order:

$\displaystyle\varepsilon^{\kappa}\partial_{t^{2}}\overline{c^{\varepsilon}}+(a+1)\partial_{t}\overline{c^{\varepsilon}}=$
$\displaystyle\varepsilon^{2-\kappa}\nabla_{x}\nabla_{x}^{T}:\overline{vv^{T}c^{\varepsilon}}-\varepsilon^{2-\kappa}\nabla_{x}\cdot\left(\chi\overline{yc^{\varepsilon}}\nabla_{x}H_{r}\star\overline{c^{\varepsilon}}+\varepsilon^{-1}\frac{n}{n+1}\mathbb{E}[q]\overline{c^{\varepsilon}}\right).$
(4.17)

Now we are ready to start the limit procedure. Let

$c^{0}:=\underset{\varepsilon\rightarrow 0}{\lim}\,c^{\varepsilon}.$

Sending $\varepsilon$ to zero in 4.13b, we obtain

$\displaystyle c^{0}=0\qquad\text{in
}[0,\infty)\times\mathbb{R}^{d}\times((\overline{B_{1}}\times\\{0,1\\})\cup(S_{1}\times(0,1))).$
(4.18)

Multiplying 4.13a by $\varepsilon^{\mu}$ and passing to the limit as
$\varepsilon\to 0$ yields an equation for $c^{0}$:

$\displaystyle\partial_{y}\left({\cal
G}_{r}[S]\overline{c^{0}}^{v}c^{0}\right)=0.$ (4.19)

Integrating 4.19 using 4.18 yields

$\displaystyle{\cal G}_{r}[S]\overline{c^{0}}^{v}c^{0}\equiv 0.$ (4.20)

To resolve 4.20 with respect to $c^{0}$, we first need to study equation

$\displaystyle{\cal G}_{r}[S]u=0$ (4.21)

for a given function $u$. Combining 4.3 and 4.4b-4.4c, ${\cal G}_{r}[S]u$ can
be expressed in terms of moments:

$\displaystyle{\cal G}_{r}[S]u(\cdot,\cdot,y^{\prime})=$ $\displaystyle{\cal
G}_{r}^{+}[S](\overline{u}^{y}-\overline{yu}^{y})-y^{\prime}\left({\cal
G}_{r}^{+}[S](\overline{u}^{y}-\overline{yu}^{y})+{\cal
G}_{r}^{-}[S]\overline{yu}^{y}\right),$ (4.22a) $\displaystyle{\cal
G}_{r}^{{\pm}}[S]u(t,x):=$
$\displaystyle\int_{\mathbb{R}^{d}}G_{{r}}^{{\pm}}[S](t,x,x^{\prime})u(t,x^{\prime})\,dx^{\prime}.$
(4.22b)

Using 4.22a, we can resolve 4.21 with respect to variable $y$ and obtain

$\displaystyle{\cal G}_{r}[S]u(\cdot,\cdot,y_{*})=0\qquad\Leftrightarrow\qquad
y_{*}={\cal Y}_{r}[S]\left(\overline{u}^{y},\overline{yu}^{y}\right),$ (4.23a)
$\displaystyle{\cal Y}_{r}[S](\mu_{0},\mu_{1}):=\frac{{\cal
G}_{r}^{+}[S](\mu_{0}-\mu_{1})}{{\cal G}_{r}^{+}[S](\mu_{0}-\mu_{1})+{\cal
G}_{r}^{-}[S]\mu_{1}},$ (4.23b)

provided that the denominator of the fraction on the right-hand side of 4.23b
is non-zero. Given that $k^{\pm}\geq 0$, we have that

$\displaystyle{\cal G}_{r}^{+}[S](\mu_{0}-\mu_{1}),{\cal
G}_{r}^{-}[S]\mu_{1}\geq 0\qquad\text{for }\mu_{0}\geq\mu_{1}\geq 0,$ (4.24)

and, moreover, if

$\displaystyle{\cal G}_{r}^{+}[S](\mu_{0}-\mu_{1})>0\qquad\text{or}\qquad{\cal
G}_{r}^{-}[S]\mu_{1}>0$ (4.25)

as well, then ${\cal Y}_{r}[S](\mu_{0},\mu_{1})$ is well-defined and

$\displaystyle{\cal Y}_{r}[S](\mu_{0},\mu_{1})\in[0,1].$

Since $y\in[0,1]$, inequalities in 4.24 are satisfied for

$\displaystyle\mu_{0}:=\overline{u}^{y},\qquad\mu_{1}:=\overline{yu}^{y}.$

Assuming that for $u:=\overline{c^{0}}^{v}$ condition 4.25 is satisfied, i.e.
that

$\displaystyle\mu_{0}:=\overline{c^{0}},\qquad\mu_{1}:=\overline{yc^{0}}$

satisfy 4.25, we can resolve 4.20 with respect to the $y$-variable and obtain

$\displaystyle c^{0}(\cdot,\cdot,\cdot,dy)=\overline{c^{0}}^{y}\delta_{{\cal
Y}_{r}[S]\left(\overline{c^{0}},\overline{yc^{0}}\right)}(y).$ (4.26)

Multiplying 4.26 by $y$ and integrating over $V\times[0,1]$, we arrive at an
equation that connects the macroscopic zero and first order $y$ moments of
$c^{0}$:

$\displaystyle\overline{yc^{0}}=\overline{c^{0}}{\cal
Y}_{r}[S]\left(\overline{c^{0}},\overline{yc^{0}}\right).$ (4.27)

Next, we integrate 4.13a with respect to $y$ over $[0,1]$ using 4.13b and pass
to the limit as $\varepsilon\rightarrow 0$ to obtain

$\displaystyle-a\nabla_{v}\cdot\left(v\overline{c^{0}}^{y}\right)=dq\overline{c^{0}}-\overline{c^{0}}^{y}.$
(4.28)

This equation can be resolved with respect to $\overline{c^{0}}^{y}$ using the
method of characteristics (a similar case was treated in [41]). The solution
reads where

$\displaystyle\Cl[xi]{xi1}(v)=\begin{cases}\frac{d}{da-1}\left(|v|^{-d+\frac{1}{a}}-1\right)&\text{for
}a\neq\frac{1}{d},\\\ -\frac{d}{a}\ln|v|&\text{for }a=\frac{1}{d}.\end{cases}$
(4.29)

Altogether, combining 4.26 and LABEL:solc0, we obtain a formula for $c^{0}$ in
terms of some of its macroscopic moments:

$\displaystyle
c^{0}(\cdot,\cdot,\cdot,y)=\overline{c^{0}}q\Cr{xi1}\delta_{{\cal
Y}_{r}[S]\left(\overline{c^{0}},\overline{yc^{0}}\right)}(y).$ (4.30)

Next, we multiply 4.28 by $v$ and $vv^{T}$, respectively, and integrate by
parts over $V$ using 4.18 to obtain formulas for macroscopic moments of order
one,

$\displaystyle
a\overline{vc^{0}}=\frac{d}{d+1}\mathbb{E}[q]\overline{c^{0}}-\overline{vc^{0}}$
$\displaystyle\Leftrightarrow\qquad$
$\displaystyle\overline{vc^{0}}=\frac{1}{a+1}\frac{d}{d+1}\mathbb{E}[q]\overline{c^{0}},$
(4.31)

and two,

$\displaystyle
2a\overline{vv^{T}c^{0}}=d\overline{vv^{T}q}\,\overline{c^{0}}-\overline{vv^{T}c^{0}}$
$\displaystyle\Leftrightarrow\qquad$
$\displaystyle\overline{vv^{T}c^{0}}=\frac{1}{2a+1}\frac{d}{d+2}\mathbb{D}[q]\overline{c^{0}}.$
(4.32)

Passing formally to the limit in 4.14 and using 4.31 we arrive at the CTE

$\displaystyle(a+1)\partial_{t}\overline{c^{0}}+\frac{d}{d+1}\nabla_{x}\cdot\left(\overline{c^{0}}\mathbb{E}[q]\right)=0\qquad\text{if
}\kappa=1.$ (4.33)

This is typical for hyperbolic scaling. Now we address the case of parabolic
scaling. Passing formally to the limit in 4.17, using 4.32, and recalling
4.27, we arrive at

$\displaystyle(a+1)\partial_{t}\overline{c^{0}}=\frac{1}{2a+1}\frac{d}{d+2}\nabla_{x}\nabla_{x}^{T}:\left(\mathbb{D}[q]\overline{c^{0}}\right)-\nabla_{x}\cdot\left(\chi\overline{yc^{0}}\nabla_{x}H_{r}\star\overline{c^{0}}\right),$
(4.34a) $\displaystyle\overline{yc^{0}}=\overline{c^{0}}{\cal
Y}_{r}[S]\left(\overline{c^{0}},\overline{yc^{0}}\right)$ (4.34b)
$\displaystyle\text{if }\kappa=2\text{ and }\mathbb{E}[q]\equiv 0.$

Similar to 3.18, equation 4.34 for the macroscopic cell density,
$\overline{c^{0}}$, includes myopic diffusion and non-local adhesion. This
time, however, the sensitivity to the adhesion force acting on the cells is
proportional to the amount of bounded CAMs, $\overline{yc^{0}}$, and, thus, to
that of the adhesion bonds formed by the cells.

Equation 4.34b appears to be of a new type. In general, it is non-local and
non-linear and cannot be explicitly resolved for $\overline{yc^{0}}$. Our next
Example deals with a special case where the equation is local and easy to
solve.

######  Example 4.5.

Let us assume that $k^{\pm}$’s are singular measures rather than functions and
of the form

$\displaystyle k^{\pm}(S,\rho):=K^{\pm}(S)\delta_{0}(\rho).$ (4.35)

This choice would correspond to the impossible situation where binding would
only occur locally. For $k^{\pm}$ from 4.35, operators ${\cal
G}_{r}^{\pm}[S]$, as defined by 4.22b and 4.4c, turn into local multiplication
operators:

$\displaystyle{\cal G}_{r}^{\pm}[S]u=K^{\pm}(S)u,$

and equation 4.34b can be easily resolved:

$\displaystyle\overline{yc^{0}}=\frac{(K^{+})^{\frac{1}{2}}}{(K^{-})^{\frac{1}{2}}+(K^{+})^{\frac{1}{2}}}(S)\overline{c^{0}}.$
(4.36)

######  Remark 4.6 (Solvability of 4.34).

System 4.34 is strongly coupled. The presence of the non-linear non-local
equation 4.34b of a new type as well as of a potentially degenerate myopic
diffusion makes its analysis challenging. We establish a result on local well-
posedness for a variant of this system which includes a quasilinear degenerate
diffusion in divergence form rather than a myopic one in a coming paper. We
expect that our analysis there will indicate a reliable approach to a
numerical treatment of 4.34.

######  Remark 4.7 (Rigorous upscaling of 4.9).

As in the case of 3.9, the upscaling provided above for 4.9 is just formal,
see Remark 3.7 for the comparison with the case we studied earlier in [41].

##  5 Simulations

In this Section, we present some simulation results for equation 3.20 and its
non-myopic modification on a one-dimensional interval. Several options have
been considered in the literature regarding appropriate boundary conditions
and treatment of the integrand of ${\cal A}_{r}$ in the area where the
argument falls outside a bounded spatial domain [23]. Here we assume no-flux
boundary conditions and extend the integrand in ${\cal A}_{r}$ by zero.
Denoting $u:=\overline{c^{0}}$, we solve numerically the IBVP

$\displaystyle 2\partial_{t}u(t,x)=\frac{1}{9}\partial_{x}\left((D_{0}+\delta
x)\partial_{x}u(t,x)\right)$
$\displaystyle\phantom{2\partial_{t}u(t,x)}+\partial_{x}\left(u(t,x)\left(\theta\frac{1}{9}\delta-\frac{1}{2}\chi_{0}\int_{(-1,1)\cap(-x,-x+6)}u(t,x+\xi)\operatorname{sign}(\xi)\,d\xi\right)\right)$
$\displaystyle\text{ in }(0,25]\times(0,6),$ (5.1a) no-flux boundary
conditions $\displaystyle\text{ on }(0,25]\\!\times\\!\\{0,6\\},$ (5.1b)
$\displaystyle u_{0}=5$ $\displaystyle\text{ in }\\{0\\}\times(0,6).$ (5.1c)

Our setting here could mimic a hypothetical in vitro experiment that starts
with loss of strong cohesion in a piece of epithelium where cells undergo a
partial EMT associated with cancer progression. Such cells acquire increased
migratory properties while retaining some cell-cell adhesiveness. This
scenario leads to enhanced migration, partly coordinated by cell-cell
adhesion, and may result in formation of cell aggregates. In the present
experiment, we assume the extracellular matrix (ECM) to be mildly
heterogeneous and not very dense. This is comparable with ECM found in healthy
epithelium prior to tissue remodelling by cancer cells. Therefore, a one-
dimensional numerical set-up, as discussed below, provides a suitable
description of an early stage of cancer invasion under uncomplicated ECM
topology. Simplifying spatial complexity and thus reducing the number of key
parameters helps elucidate the interplay between myopic diffusion and cell-
cell adhesion components of migration at this stage.

All parameter values used are collected in Subsection 5.1. Their selection has
been made for illustrative purposes and is not guided by any particular
application. The obtained results are discussed in Subsection 5.1 and the
simulation approach is explained in Subsection 5.2.

###  5.1 Results

final time | $25$
---|---
$d$ | $1$
spatial domain | $(0,6)$
$r$ | $1$
$a$ | $1$
$F(x)$ | $1$
$D[q](x)$ | $D_{0}+\delta x$
$\chi(x)$ | $\chi_{0}$
$u_{0}(x)$ | $5$
(a) Fixed model parameters.

$D_{0}$ | $\\{0.15,3,5,8\\}$
---|---
$\delta$ | $\\{0,1\\}$
$\chi_{0}$ | $\\{0.5,1\\}$
$\theta$ | $\\{0,1\\}$
(b) Varying model parameters.

size of time integration interval | $0.01$
---|---
size of time mesh for pdepe | $0.01\cdot\frac{1}{41}$
size of spatial mesh for pdepe | $0.01$
(c) Mesh sizes. Table 5.1: Parameters used in simulations of 3.20.

We singled out the following test parameters: the minimum value of the
diffusion coefficient, $D_{0}$, a diffusion perturbation parameter, $\delta$,
the constant adhesion sensitivity, $\chi_{0}$, and a binary parameter $\theta$
that if set to be zero renders the diffusion non-myopic. These variables allow
us to adjust the magnitudes of three distinct flux components: the canonical
diffusion along the spatial density gradient with diffusion coefficient
$\frac{1}{2\cdot 9}(D_{0}+\delta x)$, the advection to the left with a
constant speed $\theta\frac{1}{2\cdot 9}\delta$, and the non-local adhesion
with strength $\frac{1}{2\cdot 2}\chi_{0}$ and the distance-independent
adhesion force with $F\equiv 1$. The goal of our numerical study here is to
understand how the combination of the three motion effects can unfold.

Since initially the population is homogeneously distributed over the spatial
domain, it would have remained so in the absence of adhesion. In our
simulations, the adhesion strength $\chi_{0}$ is strictly positive, so that
adhesion is always present and, as expected, it promotes the formation of
aggregates.

In our first series of simulations we take $\theta=1$, meaning that diffusion
is myopic. The scenario involving a constant diffusion coefficient, i.e.
$\delta=0$, was previously examined in [23], and we regard it as the control
case. The results are shown in the first two columns of Subsection 5.1. In
addition to a single aggregate or a pair of aggregates displayed in [23] we
also found, for a weak diffusion corresponding to $D_{0}=0.15$ that three very
tight aggregates can be observed for $\chi_{0}=0.5$ in Subsection 5.1. For
$\chi_{0}=1$, they are about twice as dense, see Subsection 5.1. Comparing the
plots in the first column of Subsection 5.1, we see that substantially
increasing the diffusion coefficient diminishes the density of the aggregates,
slows down their formation, and reduces their number. The latter can occur
early on or later in time, as, e.g. Subsections 5.1 and 5.1 respectively
convey. The aggregates are eventually spaced more than one unit apart, as
could be expected for $r=1$. Considering that the adhesion effect increases
with growing density inside the sensing region, the impact of a density-
independent diffusion is most noticeable during the initial accumulation
phase. If diffusion brings density accumulations close enough, adhesion will
cause them to merge, as shown, e.g. in Subsection 5.1.

Next, we introduce a linear perturbation into the diffusion coefficient by
taking $\delta=1$. The results are shown in the last two columns of Subsection
5.1. This symmetry-breaking effect both increases diffusion, most notably on
the right half of the spatial domain, and adds advection to the left. Already
for small $D=0.15$, three aggregates can no longer be sustained, with the
accumulations in the middle and to the right merging into a single aggregate,
see, e.g. Subsection 5.1. Depending on $D_{0}$ and $\chi_{0}$, we observe
either two aggregates with the left one now being tighter than the right one,
compare, e.g. Subsection 5.1 with Subsection 5.1, or a single aggregate
emerging sooner or later in the left half of the domain, as can be observed,
e.g. in Subsection 5.1.

Now we set $\theta=0$ and $\delta=1$. This eliminates the constant-speed
advection to the left, rendering the diffusion non-myopic while still
retaining the spatial heterogeneity in the diffusion coefficient. The
numerical results for this case are presented in the odd columns of Subsection
5.1 next to the corresponding graphs for $\theta=1$ in the even columns of the
same Figure. The comparison of the plots reveals that advection due to
diffusion being myopic primarily serves to preserve the position of the
aggregates. Only for $D_{0}=8$ and $\chi_{0}=0.5$ it is different, see
Subsections 5.1 and 5.1. Here in the absence of the density-independent
advection in the left direction a tight aggregate fails to form. This might be
because the density accumulation fails to reach the left half of the domain
and hence stays under the influence of a strong diffusion that precludes a
substantial aggregation.

$D_{0}=$8 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=8,delta=0,chi0=0.5.jpeg}
\put(20.0,80.0){$(\chi_{0},\delta)=(0.5,0)$} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=8,delta=0,chi0=1.jpeg}
\put(20.0,80.0){$(\chi_{0},\delta)=(1,0)$} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=8,delta=1,chi0=0.5.jpeg}
\put(20.0,80.0){$(\chi_{0},\delta)=(0.5,1)$} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=8,delta=1,chi0=1.jpeg}
\put(20.0,80.0){$(\chi_{0},\delta)=(1,1)$} \end{overpic}

Figure 5.1:  Figure 5.2:  Figure 5.3:  Figure 5.4:  Figure 5.5:  Figure 5.6:
Figure 5.7:  Figure 5.8:  Figure 5.9:  Figure 5.10:  Figure 5.11:  Figure
5.12:  Figure 5.13:  Figure 5.14:  Figure 5.15:  Figure 5.16:  Figure 5.17:
Kymographs of numerical solutions of 5.1 for $\theta=1$ and various values of
$D_{0},\chi_{0}$, and $\delta$.

$D_{0}=$5 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=5,delta=0,chi0=0.5.jpeg}
\put(20.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=5,delta=0,chi0=1.jpeg}
\put(20.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=5,delta=1,chi0=0.5.jpeg}
\put(20.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=5,delta=1,chi0=1.jpeg}
\put(20.0,80.0){} \end{overpic} $D_{0}=$3 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=3,delta=0,chi0=0.5.jpeg}
\put(20.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=3,delta=0,chi0=1.jpeg}
\put(20.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=3,delta=1,chi0=0.5.jpeg}
\put(20.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=3,delta=1,chi0=1.jpeg}
\put(20.0,80.0){} \end{overpic} $D_{0}=$0.15 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=0.15,delta=0,chi0=0.5.jpeg}
\put(20.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=0.15,delta=0,chi0=1.jpeg}
\put(20.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=0.15,delta=1,chi0=0.5.jpeg}
\put(20.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=0.15,delta=1,chi0=1.jpeg}
\put(20.0,80.0){} \end{overpic}

$D_{0}=$8 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=8,delta=1,chi0=0.5.jpeg}
\put(20.0,80.0){$(\chi_{0},\theta)=$(0.5,1) } \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_NonMyopic/T=25,D=8,delta=1,chi0=0.5.jpeg}
\put(10.0,80.0){$(\chi_{0},\theta)=$(0.5,0) } \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=8,delta=1,chi0=1.jpeg}
\put(20.0,80.0){$(\chi_{0},\theta)=$(1,1) } \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_NonMyopic/T=25,D=8,delta=1,chi0=1.jpeg}
\put(10.0,80.0){$(\chi_{0},\theta)=$(1,0) } \end{overpic}

Figure 5.18:  Figure 5.19:  Figure 5.20:  Figure 5.21:  Figure 5.22:  Figure
5.23:  Figure 5.24:  Figure 5.25:  Figure 5.26:  Figure 5.27:  Figure 5.28:
Figure 5.29:  Figure 5.30:  Figure 5.31:  Figure 5.32:  Figure 5.33:  Figure
5.34: Kymographs of numerical solutions to 5.1 for $\theta=1$ (odd columns)
and $\theta=0$ (even columns), $\delta=1$, and various values of $D_{0}$ and
$\chi_{0}$.

$D_{0}=$5 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=5,delta=1,chi0=0.5.jpeg}
\put(20.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_NonMyopic/T=25,D=5,delta=1,chi0=0.5.jpeg}
\put(10.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=5,delta=1,chi0=1.jpeg}
\put(20.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_NonMyopic/T=25,D=5,delta=1,chi0=1.jpeg}
\put(10.0,80.0){} \end{overpic} $D_{0}=$3 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=3,delta=1,chi0=0.5.jpeg}
\put(20.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_NonMyopic/T=25,D=3,delta=1,chi0=0.5.jpeg}
\put(10.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=3,delta=1,chi0=1.jpeg}
\put(20.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_NonMyopic/T=25,D=3,delta=1,chi0=1.jpeg}
\put(10.0,80.0){} \end{overpic} $D_{0}=$0.15 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=0.15,delta=1,chi0=0.5.jpeg}
\put(20.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_NonMyopic/T=25,D=0.15,delta=1,chi0=0.5.jpeg}
\put(10.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_Myopic/T=25,D=0.15,delta=1,chi0=1.jpeg}
\put(20.0,80.0){} \end{overpic} 0.23
\begin{overpic}[width=476.98463pt]{MATLAB/1D_NonMyopic/T=25,D=0.15,delta=1,chi0=1.jpeg}
\put(10.0,80.0){} \end{overpic}

###  5.2 Method

Our solver for 5.1 was written in MATLAB [36]. In this Subsection, we briefly
describe our numerical scheme. To start, we decomposed the time interval
$[0,25]$ into intervals of equal size of $0.01$. Let $t_{k}$,
$k\in\\{0,1,\dots,25\times 100\\}$, be the corresponding mesh points. Starting
with $\tilde{u}_{0}:=u_{0}$, we successively determined approximations
$\tilde{u}_{k}$ of $u(t_{k},\cdot)$ by solving numerically the IBVP

$\displaystyle
2\partial_{t}\tilde{u}(t,x)=\frac{1}{9}\partial_{x}\left((D_{0}+\delta
x)\partial_{x}\tilde{u}(t,x)\right)$
$\displaystyle\phantom{2\partial_{t}\tilde{u}(t,x)}\\!\\!\\!\\!+\partial_{x}\left(\tilde{u}(t,x)\left(\theta\frac{1}{9}\delta-\frac{1}{2}\chi_{0}\int_{(-1,1)\cap(-x,-x+6)}\tilde{u}_{k}(x+\xi)\operatorname{sign}(\xi)\,d\xi\right)\right)$
$\displaystyle\text{ in }\left(0,0.01\right]\times(0,6),$ (5.18a) no-flux
boundary conditions $\displaystyle\text{ on
}\left(0,0.01\right]\\!\times\\!\\{0,6\\},$ (5.18b)
$\displaystyle\tilde{u}=\tilde{u}_{k}$ $\displaystyle\text{ in
}\\{0\\}\times(0,6),$ (5.18c)

and then set

$\displaystyle\tilde{u}_{k+1}:=\tilde{u}(0.01,\cdot).$

To solve the implicit-explicit equation 5.18, we used MATLAB’s function pdepe
[38]. The sizes of the equidistant meshes that we passed to pdepe were
$1/41\cdot 0.01$ for time and $0.01$ for space. To apply pdepe in this case,
one needs first to produce the function described by the non-local term, as it
is part of the expression for the flux. We discretised that integral term for
the entire spatial domain replacing it by the right Riemann sums

$\displaystyle
s_{m}:=0.01\sum_{l=a(m)}^{b(m)}\tilde{u}_{k}\left(0.01(m+l)\right)\operatorname{sign}\left(0.01l\right),$
(5.19)

where

$\displaystyle a(m):=\max\\{-99,-m\\},\quad b(m):=\min\\{100,-m+100\cdot
6\\},\qquad m\in\\{0,\dots,100\cdot 6\\},$

and then interpolated these values using MATLAB’s function griddedInterpolant
[37] to create a function of $x$. To compute the sums 5.19, we applied
MATLAB’s function xcorr [39] as follows. Let

$\displaystyle\omega:=(\underbrace{-1,\dots,-1}_{(100-1)\text{
times}},0,\underbrace{1,\dots,1}_{100\text{ times}}),$ $\displaystyle
u:=(\tilde{u}_{k}(0),\tilde{u}_{k}(0.01),\dots,\tilde{u}_{k}(6)),$
$\displaystyle w:=\texttt{xcorr}(u,\omega),$

then

$\displaystyle(s_{0},\dots,s_{100\cdot 6})=0.01\left(w_{(100\cdot
6+1)-(100-1)},\dots,w_{2(100\cdot 6+1)-100}\right).$

##  6 Discussion

CCA plays a pivotal role in the development and functioning of multicellular
organisms. Notably, it regulates cell migration, either promoting or
inhibiting it. Macroscopic mathematical models can contribute to a better
understanding of adhesion effects because they are amenable to both rigorous
mathematical analysis and in silico studies, and numerical results for these
models can be compared to medical images.

In this paper we aimed at devising a new multiscale approach to modelling cell
migration driven by such effects as CCA and anisotropic diffusion. After
reviewing previously available approaches in Subsection 1.2, we derived in
Section 3 classes of IPDE models containing non-local adhesion and myopic
diffusion, the classical model in [4] being their special case. We further
extended our approach in Section 4 where we derived a novel model 4.34 that
can account for subcellular binding dynamics of CAMs, molecules responsible
for cells sticking to each other.

Our modelling can serve as a starting point for considerably more realistic
models for adhesion-driven motion. One of our simplifying modelling
assumptions, (Ai) in Section 4, was that all junctions are cell-cell junctions
of one and the same type. The main application that we had in mind here were
the adherens junctions such as facilitated by E-cadherins, see Subsection 1.1.
In reality, cells form a variety of junctions and, moreover, subpopulations
with distinct adhesion properties can be involved. Our approach could be
extended to accommodate such complexities.

The usefulness of such models as 4.34 hinges on their solvability. As
previously announced in Remark 4.6, we prove local solvability, yet for a
quasilinear non-myopic degenerate diffusion, in a coming paper. We postpone to
that work a discussion of challenges that arise in connection with treatment
of equation 4.34b. While non-myopic diffusion is often adopted in modelling
cell motion, our simulation results in Section 5 underscore the difference
that a myopic diffusion can make. We are going to settle the solvability for
4.34 in a future work.

## Acknowledgement

  * •

The authors thank Christina Surulescu (RPTU Kaiserslautern-Landau) for
stimulating discussions.

  * •

The authors were supported by the Engineering and Physical Sciences Research
Council [grant number EP/T03131X/1].

  * •

For the purpose of open access, the authors have applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising.

  * •

All data is provided in full in Section 5 of this paper.

## References

  * [1] Bruce Alberts et al.  “Molecular biology of the cell”  Garland Science, 2008
  * [2] Luigi Ambrosio  “Transport equation and Cauchy problem for $BV$ vector fields”  In _Invent. Math._ 158.2, 2004, pp. 227–260  DOI: 10.1007/s00222-004-0367-2
  * [3] Luigi Ambrosio, Nicola Fusco and Diego Pallara  “Functions of bounded variation and free discontinuity problems”, Oxford Mathematical Monographs  The Clarendon Press, Oxford University Press, New York, 2000, pp. xviii+434
  * [4] Nicola J. Armstrong, Kevin J. Painter and Jonathan A. Sherratt  “A continuum approach to modelling cell-cell adhesion”  In _J. Theoret. Biol._ 243.1, 2006, pp. 98–113  DOI: 10.1016/j.jtbi.2006.05.030
  * [5] George I. Bell  “Models for the Specific Adhesion of Cells to Cells”  In _Science_ 200.4342, 1978, pp. 618–627  DOI: 10.1126/science.347575
  * [6] Nicola Bellomo, Abdelghani Bellouquid, L. Gibelli and N. Outada  “A Quest Towards a Mathematical Theory of Living Systems”  Birkhäuser, 2018
  * [7] Andreas Buttenschön, Thomas Hillen, Alf Gerisch and Kevin J. Painter  “A space-jump derivation for non-local models of cell-cell adhesion and non-local chemotaxis”  In _J. Math. Biol._ 76.1-2, 2018, pp. 429–456  DOI: 10.1007/s00285-017-1144-3
  * [8] Yu Cai, Nitesh Shashikanth, Deborah E Leckband and Daniel K Schwartz  “Cadherin diffusion in supported lipid bilayers exhibits calcium-dependent dynamic heterogeneity”  In _Biophysical journal_ 111.12 Elsevier, 2016, pp. 2658–2665  DOI: 10.1016/j.bpj.2016.10.037
  * [9] José A. Carrillo, Raluca Eftimie and Franca Hoffmann  “Non-local kinetic and macroscopic models for self-organised animal aggregations”  In _Kinet. Relat. Models_ 8.3, 2015, pp. 413–441  DOI: 10.3934/krm.2015.8.413
  * [10] A. Chauvière, T. Hillen and L. Preziosi  “Modeling cell movement in anisotropic and heterogeneous network tissues”  In _Networks & Heterogeneous Media_ 2.2 American Institute of Mathematical Sciences (AIMS), 2007, pp. 333–357  DOI: 10.3934/nhm.2007.2.333
  * [11] Li Chen, Kevin Painter, Christina Surulescu and Anna Zhigun  “Mathematical models for cell migration: a non-local perspective”  In _Philosophical Transactions of the Royal Society B_ 375.1807 The Royal Society, 2020, pp. 20190379  DOI: 10.1098/rstb.2019.0379
  * [12] Gregor Corbin et al.  “Modeling glioma invasion with anisotropy- and hypoxia-triggered motility enhancement: from subcellular dynamics to macroscopic PDEs with multiple taxis”  In _Math. Models Methods Appl. Sci._ 31.1, 2021, pp. 177–222  DOI: 10.1142/S0218202521500056
  * [13] Anne Dietrich, Niklas Kolbe, Nikolaos Sfakianakis and Christina Surulescu  “Multiscale modeling of glioma invasion: from receptor binding to flux-limited macroscopic PDEs”  In _Multiscale Model. Simul._ 20.2, 2022, pp. 685–713  DOI: 10.1137/21M1412104
  * [14] R. L. Dobrušin  “Vlasov equations”  In _Funktsional. Anal. i Prilozhen._ 13.2, 1979, pp. 48–58, 96
  * [15] Maria Eckardt, Kevin J. Painter, Christina Surulescu and Anna Zhigun  “Nonlocal and local models for taxis in cell migration: a rigorous limit procedure”  In _J. Math. Biol._ 81.6-7, 2020, pp. 1251–1298  DOI: 10.1007/s00285-020-01536-4
  * [16] Raluca A. Eftimie  “Modeling group formation and activity patterns in self-organizing communities of organisms” Thesis (Ph.D.)–University of Alberta (Canada)  ProQuest LLC, Ann Arbor, MI, 2008, pp. 174  URL: http://gateway.proquest.com/openurl?url_ver=Z39.88-2004&rft_val_fmt=info:ofi/fmt:kev:mtx:dissertation&res_dat=xri:pqdiss&rft_dat=xri:pqdiss:NR46312
  * [17] Christian Engwer, Thomas Hillen, Markus Knappitsch and Christina Surulescu  “Glioma follow white matter tracts: a multiscale DTI-based model”  In _Journal of Mathematical Biology_ 71.3, 2015, pp. 551–582  DOI: 10.1007/s00285-014-0822-7
  * [18] Peter Friedl and Stephanie Alexander  “Cancer invasion and the microenvironment: plasticity and reciprocity”  In _Cell_ 147.5 Elsevier, 2011, pp. 992–1009  DOI: 10.1016/j.cell.2011.11.016
  * [19] A. Gerisch and M. A. J. Chaplain  “Mathematical modelling of cancer cell invasion of tissue: local and non-local models and the effect of adhesion”  In _J. Theoret. Biol._ 250.4, 2008, pp. 684–704  DOI: 10.1016/j.jtbi.2007.10.026
  * [20] F. Golse  “Mean field kinetic equations”, Lecture notes at the École Polytechnique, Univ. Paris-Saclay, 2013
  * [21] Douglas Hanahan and Robert A Weinberg  “Hallmarks of cancer: the next generation”  In _cell_ 144.5 Elsevier, 2011, pp. 646–674  DOI: 10.1016/j.cell.2011.02.013
  * [22] Frederic Heihoff  “Global solutions to a haptotaxis system with a potentially degenerate diffusion tensor in two and three dimensions”  In _Nonlinearity_ 36.2, 2023, pp. 1245–1278
  * [23] Thomas Hillen and Andreas Buttenschön  “Nonlocal adhesion models for microorganisms on bounded domains”  In _SIAM J. Appl. Math._ 80.1, 2020, pp. 382–401  DOI: 10.1137/19M1250315
  * [24] Thomas Hillen and Kevin J. Painter  “Transport and anisotropic diffusion models for movement in oriented habitats”  In _Dispersal, individual movement and spatial ecology_ 2071, Lecture Notes in Math.  Springer, Heidelberg, 2013, pp. 177–222  DOI: 10.1007/978-3-642-35497-7˙7
  * [25] Pierre-Emmanuel Jabin and Zhenfu Wang  “Mean field limit and propagation of chaos for Vlasov systems with bounded forces”  In _J. Funct. Anal._ 271.12, 2016, pp. 3588–3627  DOI: 10.1016/j.jfa.2016.09.014
  * [26] Jan Kelkel and Christina Surulescu  “A multiscale approach to cell migration in tissue networks”  In _Math. Models Methods Appl. Sci._ 22.3, 2012, pp. 1150017, 25  DOI: 10.1142/S0218202511500175
  * [27] Jan Kelkel and Christina Surulescu  “On some models for cancer cell migration through tissue networks”  In _Math. Biosci. Eng._ 8.2, 2011, pp. 575–589  DOI: 10.3934/mbe.2011.8.575
  * [28] Nadia Loy and Luigi Preziosi  “Kinetic models with non-local sensing determining cell polarization and speed according to independent cues”  In _J. Math. Biol._ 80.1-2, 2020, pp. 373–421  DOI: 10.1007/s00285-019-01411-x
  * [29] Nadia Loy and Luigi Preziosi  “Modelling physical limits of migration by a kinetic model with non-local sensing”  In _J. Math. Biol._ 80.6, 2020, pp. 1759–1801  DOI: 10.1007/s00285-020-01479-w
  * [30] Alistair M. Middleton, Christian Fleck and Ramon Grima  “A continuum approximation to an off-lattice individual-cell based model of cell migration and adhesion”  In _J. Theoret. Biol._ 359, 2014, pp. 220–232  DOI: 10.1016/j.jtbi.2014.06.011
  * [31] Daniela Morale, Vincenzo Capasso and Karl Oelschläger  “A rigorous derivation of a nonlinear integro-differential equation from a system of stochastic differential equations for an aggregation model”, Preprint / Interdisziplinäres Zentrum für Wissenschaftliches Rechnen, IWR ARRAY(0x562b2b01dfb8)  Heidelberg: IWR, SFB 359, 1998, pp. 17 pp.
  * [32] Daniela Morale, Vincenzo Capasso and Karl Oelschläger  “An interacting particle system modelling aggregation behavior: from individuals to populations”  In _J. Math. Biol._ 50.1, 2005, pp. 49–66  DOI: 10.1007/s00285-004-0279-1
  * [33] H.G. Othmer and T. Hillen  “The diffusion limit of transport equations II: chemotaxis equations”  In _SIAM J. Appl. Math._ 62, 2002, pp. 1122–1250
  * [34] Benoît Perthame, Weiran Sun and Min Tang  “The fractional diffusion limit of a kinetic model with biochemical pathway”  In _Z. Angew. Math. Phys._ 69.3, 2018, pp. Paper No. 67, 15  DOI: 10.1007/s00033-018-0964-3
  * [35] Benoit Perthame, Weiran Sun, Min Tang and Shugo Yasuda  “Multiple asymptotics of kinetic equations with internal states”  In _Math. Models Methods Appl. Sci._ 30.6, 2020, pp. 1041–1073  DOI: 10.1142/S0218202520400060
  * [36] The MathWorks Inc.  “MATLAB version: 9.14 (R2023a)”  Natick, Massachusetts, United States: The MathWorks Inc., 2023  URL: https://www.mathworks.com
  * [37] The MathWorks Inc.  “griddedInterpolant documentation”  Natick, Massachusetts, United States: The MathWorks Inc., 2023  URL: https://uk.mathworks.com/help/matlab/ref/griddedinterpolant.html
  * [38] The MathWorks Inc.  “pdepe documentation”  Natick, Massachusetts, United States: The MathWorks Inc., 2023  URL: https://uk.mathworks.com/help/matlab/ref/pdepe.html
  * [39] The MathWorks Inc.  “xcorr documentation”  Natick, Massachusetts, United States: The MathWorks Inc., 2023  URL: https://uk.mathworks.com/help/matlab/ref/xcorr.html
  * [40] Michael Winkler and Christina Surulescu  “Global weak solutions to a strongly degenerate haptotaxis model”  In _Commun. Math. Sci._ 15.6, 2017, pp. 1581–1616  DOI: 10.4310/CMS.2017.v15.n6.a5
  * [41] Anna Zhigun and Christina Surulescu  “A novel derivation of rigorous macroscopic limits from a micro-meso description of signal-triggered cell migration in fibrous environments”  In _SIAM J. Appl. Math._ 82.1, 2022, pp. 142–167  DOI: 10.1137/20M1365442

## Appendix A

The following Lemma extends the result that is well-known for skew-symmetric
kernels $K$, see, e.g. [20, Chapter 3 §2 Theorem 3.2.1(b)]. The proof, which
is very similar, is provided for the reader’s convenience.

######  Lemma A.1.

Let $m,N\in\mathbb{N}$. Let $A\in
C([0,\infty)\times\mathbb{R}^{m};\mathbb{R}^{m})$ and $K\in
C([0,\infty)\times\mathbb{R}^{m}\times\mathbb{R}^{m};\mathbb{R}^{m})$. Define

$\displaystyle{\cal
K}\mu(t,Z):=\int_{\mathbb{R}^{m}}K(t,Z,Z^{\prime})\,\mu(dZ^{\prime}),\qquad\operatorname{tr}K(t,Z):=K(t,Z,Z),\qquad
Z\in\mathbb{R}^{m}.$ (A.1)

Consider the ODE system

$\displaystyle\frac{dZ_{i}}{dt}=\frac{1}{N}\sum_{\underset{j\neq
i}{j=1}}^{N}K(\cdot,Z_{i},Z_{j})+A(\cdot,Z_{i}).$ (A.2)

Let $(Z_{1}^{T},\dots,Z_{N}^{T})^{T}\in C^{1}([0,T);\mathbb{R}^{mN})$ be a
solution to A.2 for some $T\in[0,\infty]$. Then, the empirical measure

$\displaystyle\mu_{N}:=\frac{1}{N}\sum_{i=1}^{N}\delta_{Z_{i}}$

satisfies

$\displaystyle\partial_{t}\mu_{N}+\nabla_{z}\cdot\left(\left({\cal
K}\mu_{N}-\frac{1}{N}\operatorname{tr}K+A\right)\mu_{N}\right)=0$ (A.3)

on $(0,T)\times\mathbb{R}^{m}$ in the distributional sense.

###### Proof.

Throughout the proof, equations are understood to be satisfied in the sense of
distributions on $(0,T)\times\mathbb{R}^{m}$.

Let $i\in\\{1,\dots,N\\}$. From the theory of CTEs it is known that (see, e.g.
[20]) that

$\displaystyle\partial_{t}\delta_{Z_{i}}+\nabla_{z}\cdot\left(\left(\frac{1}{N}\sum_{\underset{j\neq
i}{j=1}}^{N}K(\cdot,\cdot,Z_{j})+A\right)\delta_{Z_{i}}\right)=0.$ (A.4)

Using A.1 and the property

$\varphi\delta_{\bar{Z}}=\varphi(\bar{Z})\delta_{\bar{Z}}$

which holds for any $\bar{Z}\in\mathbb{R}^{m}$ and continuous function
$\varphi=\varphi(Z)$, we can rewrite A.4 as follows:

$\displaystyle 0=$
$\displaystyle\partial_{t}\delta_{Z_{i}}+\nabla_{z}\cdot\left(\left(\frac{1}{N}\sum_{\underset{j\neq
i}{j=1}}^{N}K(\cdot,\cdot,Z_{j})+A\right)\delta_{Z_{i}}\right)$
$\displaystyle=$
$\displaystyle\partial_{t}\delta_{Z_{i}}+\nabla_{z}\cdot\left(\left(\frac{1}{N}\sum_{{j=1}}^{N}{\cal
K}\delta_{Z_{j}}-\frac{1}{N}\operatorname{tr}K+A\right)\delta_{Z_{i}}\right)$
$\displaystyle=$
$\displaystyle\partial_{t}\delta_{Z_{i}}+\nabla_{z}\cdot\left(\left({\cal
K}\mu_{N}-\frac{1}{N}\operatorname{tr}K+A\right)\delta_{Z_{i}}\right).$ (A.5)

Computing the arithmetic mean with respect to $i$ on both sides of A.5 yields
A.3. ∎