# Numerical Approaches for non-local Transport-Dominated PDE Models with Applications to Biology

Johan Marguet, Raluca Eftimie, Alexei Lozinski

► To cite this version:

Johan Marguet, Raluca Eftimie, Alexei Lozinski. Numerical Approaches for non-local Transport-Dominated PDE Models with Applications to Biology. Computational and Applied Mathematics Journal, In press. hal-04989162

HAL Id: hal-04989162

https://hal.science/hal-04989162v1

Submitted on 20 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Numerical Approaches for non-local
Transport-Dominated PDE Models with
Applications to Biology

Johan Marguet¹, Raluca Eftimie¹*†, Alexei Lozinski¹*†

¹Laboratoire de Mathématiques de Besançon, 16 Route de Gray,
Besançon, 25000, France.

*Corresponding author(s). E-mail(s): raluca.eftimie@univ-fcomte.fr; alexei.lozinski@univ-fcomte.fr;

Contributing authors: johan.marguet@univ-fcomte.fr;

†These authors contributed equally to this work.

**Abstract**

Transport-dominated partial differential equation models have been used exten-
sively over the past two decades to describe various collective migration phenom-
ena in cell biology and ecology. To understand the behaviour of these models
(and the biological systems they describe) different analytical and numerical
approaches have been used. While the analytical approaches have been discussed
by different recent review studies, the numerical approaches are still facing dif-
ferent open problems, and thus are being employed on a rather ad-hoc basis for
each developed non-local model. The goal of this review is to summarise the basic
ideas behind these transport-dominated non-local models, to discuss the current
numerical approaches used to simulate these models, and finally to discuss some
open problems related to the applications of these numerical methods, in par-
ticular the finite element method. This allows us to emphasize the opportunities
offered by this numerical method to advance the research in this field. In addi-
tion, we present in detail some numerical schemes that we used to discretize these
non-local equations; in particular a new semi-implicit scheme we introduced to
stabilize the oscillations obtained with classical schemes.

**Keywords:** Non-local Models, Biology and Ecology, Numerical Approaches, Finite Elements
# 1 Introduction

Collective behaviours of cells, bacteria, insects, animals and even humans [1] have attracted human attention for at least two millennia [2], due to the spatial and spatio-temporal patterns exhibited by these organisms. Many of these collective behaviours are leading to the directed, coordinated movement of organisms: from the directed migration of cells during health (e.g., tissue morphogenesis [3, 4]) and disease (e.g., cancer [3], abnormal wound healing [5]), to the directed movement of flocks of birds [6] or herds of ungulates [7] during their annual migrations or swarm of bees [8]. The interactions between the organisms inside these biological/ecological aggregations can be either local (when organisms interact only with immediate neighbours), or non-local (when organisms can interact with others further away), see Figure 1. For more

<figure><img src="image_2.png" /><figcaption>**Fig. 1** Difference between local and non-local interactions at the cellular level. In both cases, there are two types of interaction: bio-mechanical (represented by the solid and striped arrows) and bio-chemical (represented by bright pink molecules). On the left, we say that these interactions are local since they appear between directly neighboring cells. Whereas on the right, these interactions are non-local since the mechanical traction of cell 1 towards the left leads to the movement towards the left of cells 2 and 4 which are neighbors but also of cells 3 and 5 which are not directly neighboring. In addition, cell 1 can chemically interact with distant cells 3 and 5 via protrusions carrying a chemical messenger.</figcaption></figure>

than half a century, researchers have acknowledged and further investigated non-local social interactions that are the result of long-distance communication between the different members of these biological/ecological aggregations: from non-local animal-animal interactions due to long-distance auditory, chemical or visual signals [9–11], to non-local cell-cell interactions via bio-chemical signals (e.g., hormones [12] and airinemes [13]) and bio-mechanical forces (e.g., via cellular traction stresses exerted
through the substrate [14]). In the context of these experimental research directions, the last two-three decades have witnessed the development of numerous mathematical models aimed at describing such non-local interactions across various cell and animal aggregations: from discrete (individual-based) models, that focus on changes in the position and velocity of individual organisms [15–18], to continuum (integro-partial differential equation) models that focus on changes on the densities of populations of organisms [19, 20]. While the discrete models have been investigated mainly computationally [1, 17, 21] (due to a lack of methodological approaches for their analytical investigation), the continuum models have been investigated both analytically and computationally [19, 20, 22]. However, the details of the numerical methods used to discretize these continuum models were not always available, since many researchers have been using classical approaches (i.e., finite difference and finite volume schemes combined with standard quadrature methods for the computation of integrals), which can be found in classical numerical analysis books [23]. Nevertheless, this lack of details might impede some (early-career) researchers to further develop and investigate numerically various mathematical models that incorporate long-distance biological/ecological interactions. For this reason, in this study we aim to summarise different numerical approaches that could be used to simulate the spatial and spatio-temporal dynamics exhibited by this class of non-local PDE models in 1D and 2D.

**Terminology.** Before presenting briefly the various types of non-local mathematical models that have been developed over the years, we need to clarify the terminology used throughout this study. Random cell/bacterial/animal movement is described by *diffusion terms*. Directed cell/bacterial/animal movement in response to a stimulus (e.g., food, stiffer substrates) is described by *taxis terms*. In the literature there are different types of taxis: *chemotaxis* (movement aligned with a chemical gradient), *haptotaxis* (cell movement up a gradient of substrate-bound molecules), *durotaxis* (cell movement up a gradient of extracellular matrix stiffness) [24]. All these taxis behaviors are modelled by similar terms: $\nabla \cdot (\chi u \nabla c)$, where $c$ describes the concentration/density of chemicals/molecules/etc. that determine the gradient $\nabla c$, $u$ is the density of individuals/cells/animals that move towards/away the chemical gradient $\nabla c$, and $\chi$ is the chemotactic sensitivity coefficient (that might be a constant or might depend on the concentration of the chemical $c$ or on other factors in the environment). Due to this directional transport of individuals/cells/animals with density $u$, these *taxis models* are also referred to as *advection models* or *transport models*. While this terminology is common for local models [25, 26], it has been used also to describe non-local fluxes that are the result of social interactions (e.g., attractive movement towards higher densities of conspecifics located further away), as those modeled in equation (2) [27]. Nevertheless, some studies, mainly applied to cell biology, refer to these non-local models also as *adhesion models*, since the non-local operator describes cell-cell or cell-matrix adhesion [27, 28]. Other studies, mainly applied to ecology, refer to these non-local models as *aggregation models* [29–33].

The non-local interactions between different members of the same social community (e.g., a cell population or an animal population), which are located at various spatial
positions (some closer, some further away) are described by integrals over different spatial ranges. These integrals can be included into the mathematical equations describing the changes in the cell/animal population densities through various terms that have different biological meanings:

*   **Non-local interactions in the reaction part** [34], to describe the information received by an organism about its possibility of proliferation/elimination within a spatial range around its current position. The general structure of these models is:

$$ \frac{\partial u}{\partial t} = \nabla \cdot (D \nabla u(t, \mathbf{x})) + a[u], \qquad (1) $$

where $u(\mathbf{x}, t)$ is the density of cells/organisms at position $\mathbf{x} \in \mathbb{R}^n$, $n \in \{1, 2, 3\}$, and time $t$, $D$ is the diffusion coefficient (which can be constant or density-dependent), $a[u]$ is the non-local source term often in the form $a[u] = F(u) + I[g(u)]$. Here, $F(u)$ describes local factors that impact proliferation/elimination and $I[g(u)] = \int_{\Omega} g(x, u(y, t)) dy$ describes non-local factors that impact proliferation/elimination (with $g$ a linear or non-linear function for the elimination/proliferation of population $u$, and $\Omega$ the domain of interest for these non-local interactions, which might not be the whole $\mathbb{R}^n$).

*   **Non-local interactions in the advection part** [27, 35, 36], to describe the directional flux generated by interactions with conspecifics positioned further away from a reference individual. The terminology used for this class of non-local models, that have more-or-less the same structure, is quite varied: from *non-local advection models* for non-local taxis (e.g., haptotaxis or chemotaxis) [37], to *aggregation models* [30, 31, 38, 39]. The general form of these models is:

$$ \frac{\partial u}{\partial t} = \nabla \cdot (D \nabla u(t, \mathbf{x}) - \beta u(t, \mathbf{x}) K[u]) + F(u), \qquad (2) $$

where $u(\mathbf{x}, t)$ is the density of cells/organisms at position $\mathbf{x} \in \mathbb{R}^n$, $n \in \{1, 2, 3\}$, and time $t$, $D$ is the diffusion coefficient (which can be constant or density-dependent [37]), $\beta$ is the advection coefficient (which can also be constant or density-dependent), $K[u]$ is the non-local operator and $F(u)$ is the reaction term. In these models, the non-local operator $K[u]$ describes the relative force of interactions exerted by a cell/organism at position $\mathbf{y}$ on a cell/organism at position $\mathbf{x}$. There is also a second class of models that include non-local terms in the advection part: not directly through the convolution operator (as seen in equation (2)), but rather through the gradient of the convolution [29]:

$$ \frac{\partial u}{\partial t} = \nabla \cdot (D \nabla u(t, \mathbf{x}) - \beta u(t, \mathbf{x}) \nabla M[u]) + F(u). \qquad (3) $$

Unlike the previous model, in equation (3) the interaction potential $M[u]$ represents a non-local measure of the cell population, the gradient of which influences the movement of the population [27]. The other terms and parameters in equation
(3) have the same meaning as in equation (2).

Remark 1 Note that for some specific cases, models (2) and (3) can turn out to be equivalent. For example, if in (2) we choose the following 1D kernel

$$K[u] = \hat{\omega} * u \text{ with } \hat{\omega}(x) = \begin{cases} 1, & \text{if } 0 < x < r, \\ -1, & \text{if } -r < x < 0, \\ 0, & \text{else,} \end{cases}$$

and if in (3) we choose the 1D kernel

$$\mathcal{M}[u] = \tilde{\omega} * u \text{ with } \tilde{\omega}(x) = \begin{cases} |x|, & \text{if } -r < x < r, \\ r, & \text{else,} \end{cases}.$$

then we observe that we have $\partial_x \mathcal{M}[u] = K[u]$. Thus, models (2) and (3) are equivalent.

Remark 2 We need to mention that there are also other classes of models that consider non-local diffusion processes that result from particles/individuals jumping in/out of neighbouring positions [40]. An example of such a model is given in [41]:

$$\frac{\partial u(t, \mathbf{x})}{\partial t} = \left( \int_{\mathbb{R}^n} \mathcal{J}(\mathbf{x} - \mathbf{y})u(t, \mathbf{y})d\mathbf{y} - u(t, \mathbf{x}) \right) + F(u), \quad (4)$$

where $\mathcal{J}(\mathbf{x} - \mathbf{y})$ describes the probability distribution of jumping from position $\mathbf{y} \in \mathbb{R}^n$ to position $\mathbf{x} \in \mathbb{R}^n$ [41]. Hence, the convolution $\mathcal{J} * u = \int_{\mathbb{R}^n} \mathcal{J}(\mathbf{x} - \mathbf{y})u(t, \mathbf{y})d\mathbf{y}$ gives the rate at which individuals arrive at position $\mathbf{x}$ from all other positions $\mathbf{y}$, while $-u(t, \mathbf{x})$ gives the rate at which individuals are leaving position $\mathbf{x}$.

Note that these non-local terms that describe jumping between different spatial locations could be interpreted as models with non-local interactions incorporated into the source term, as we can write $a[u] = \mathcal{J} * u + F(u)$.

In this review we focus only on models with non-local terms inside the advection component since, from a numerical perspective, they pose particular challenges due to the advection operator combined with the particularities of biologically-inspired kernels (which can have singularities if one assumes very strong repulsion when $\mathbf{y} = \mathbf{x}$, to avoid cells/organisms occupying the same spatial position). As mentioned before, since the published studies are usually lacking detailed discussions about the various numerical approaches used to discretise and simulate these non-local models, the goal of this review is to fill in this gap and to summarise these different numerical approaches while discussing their advantages and drawbacks, and emphasising the open numerical questions in the field. This way, we hope to provide a more structured approach to choosing appropriate numerical schemes for those non-local models, where the non-local interactions lead to an advective flux.

The rest of this review is structured as follows: in Section 3 we summarise briefly some of the analytical approaches used to investigate these non-local models (we do not focus more on this aspect, since a very detailed review of this topic was recently published [42]). In Section 4 we summarise the numerical approaches (mainly based on finite differences and finite volume schemes) used to discretise these models. Since
the finite element (FE) approaches are almost nonexistent, even if they are well-
suited for 2D and 3D models that are important to describe the dynamics of tissues,
in Section 5 we will summarise the finite element method, and present a series of
simulations for the dynamics of some very simple 1D and 2D models. In Section 6 we
discuss some open problems related to numerical approaches for this class of non-local
models. We conclude this review with Section 7, where we present a final overview of
these numerical methods, and mention some future research perspectives triggered by
the recent advances on non-local models with applications to biology and medicine.

## 2 Brief discussion of non-local transport-dominated models in biology and ecology

In this section we discuss in more detail the biological background behind some simple non-local models developed in the context of cell biology (see Section 2.1) and ecology (see Section 2.2). We also show the general form of the equations derived to describe these biological and ecological interactions. By presenting them in parallel, we aim to emphasise the similarities between the models used in these two different domains. While these aspects were discussed in more detail in other reviews [27, 42, 43], it is important to present them here as well (but very briefly), for the completeness of the discussion. Especially since these details are not always known to researchers in the numerical analysis field.

### 2.1 Non-local models for cell biology

In cell biology, cell-cell and cell-matrix adhesion are crucial processes in the organiza-
tion and maintenance of tissues during development (e.g., morphogenesis) and during
pathological situations (e.g., cancer, wound healing) [4, 44]. While these processes
occur locally via adhesion molecules such as integrins and cadherins [45], their effects
are perceived not only locally but also at a distance (see also Figure 1). In general,
cell-cell interactions can occur over different spatial ranges:

* short-range repulsion: cells repel each other to avoid overcrowding [46];

* medium-range alignment: cells tend to align due to traction forces [47], or due to substrate topography [48], which can align the actin fibres in the cell cytoskeleton;

* long-range attraction: cells can detect and interact with other cells positioned at large distances due to adhesion forces exerted by those cells, or through the generation of cell protrusions [49].

The need to consider long-distance cell-cell interactions to understand collective cell
dynamics is supported by various experimental studies. For example, a 2009 study
[47] showed that traction forces are acting many cell rows behind the leading edge of
the invading cell sheet (as they were measured up to 200µm, equivalent to roughly 10
rows of cells, each with an average diameter of 20µm). In addition, cells can extend
membrane protrusions over distances of up to ten times their diameter or even more,
as is the case with cytonemes (that can reach up to 700µm in length [49]), tunnelling
nanotubes and tumour microtubes (that can reach up to 100µm in length [50]). These
membrane protrusions are used in long-distance cell-cell communication and movement coordination, via physical contact between cells positioned far away from each other [50] (see also Figure 1).

All these (and many other) biological experiments are emphasising the importance of incorporating non-local cell-cell and cell-ECM interactions that generate different force distributions, into the mathematical models for collective cell movement. To exemplify such a simple model, let us now define $u(\mathbf{x}, t)$ to be the density of a cell population at time $t$ and position $\mathbf{x} \in \mathbb{R}^n$ with $n \in \{1, 2, 3\}$. The evolution of this density can be given by the following generic model:

$$ \partial_t u = -\nabla \cdot J(\mathbf{x}, t), $$

where $J$ is the population flux resulting from the movement of cells. This flux can be decomposed in two parts $J = J_d + J_a$, where $J_d$ is the diffusive part and $J_a$ is the advective part:

* The diffusive part describes the local random movement of cells, which causes them to spread out in space following Fick's law [51]. This spreading is described by the gradient of the density $u$, multiplied by a diffusion coefficient $D$ (here assumed to be constant):

$$ J_d = -D\nabla u. \qquad (5) $$

* For the advective part, if we assume that a cell located at position $\mathbf{x}$ interacts with another cell located at position $\mathbf{y} \neq \mathbf{x}$, the movement generated by this interaction is the result of the traction forces caused by cell-cell adhesive bonds. Moreover, these adhesive bonds lead to a sequential transmission of traction forces over larger spatial distances (see also Figure 1), and thus cells can sense the biomechanical effect of other cells at a distance. Therefore, we can assume that the movement of a generic cell located at $\mathbf{x}$ results from the sum of all interactions between the cell positioned at $\mathbf{x}$ with all other cells located within its perception radius $r$. Mathematically, this translates into the following advective flux:

$$ J_a = u(\mathbf{x}, t) \int_{B_r} \alpha(\mathbf{y}) \mathbf{K}_r(\mathbf{x}, \mathbf{y}) f(u(\mathbf{y}, t)) \, d\mathbf{y}, \qquad (6) $$

where $\alpha$ is the adhesion coefficient which could vary depending on the distance separating two cells. The product of $\mathbf{K}_r$ (a vector valued kernel) with $f$ (a scalar valued function) describes the induced direction of the cell movement and its magnitude.

Hence, the generic model can take the following form [19]:

$$ \partial_t u = D\Delta u - \nabla \cdot \left( u(\mathbf{x}, t) \int_{B_r} \alpha(\mathbf{y}) \mathbf{K}_r(\mathbf{x}, \mathbf{y}) f(u(\mathbf{y}, t)) \, d\mathbf{y} \right) + F(u), \qquad (7) $$

where function $F$ was added here to describe the proliferation and decay of cells, to ensure a more general although still simple model. As it was shown in [19, 27], the coefficient $\alpha$ has a crucial importance for cellular aggregation: if this $\alpha$ coefficient is
below a certain threshold, then a dispersed cell population remains dispersed over time; whereas if $\alpha$ is above this threshold, the population could aggregate.

Furthermore, when modelling the interactions between different populations of cells inside the same tissue, one needs to consider systems of coupled equations (where coupling can be done through local and/or non-local terms). Thus, model (7) can be generalised to describe the dynamics of $p \in \mathbb{N}^+$ cell populations:

$$ \partial_t u_i = D_i \Delta u_i - \sum_{j=1}^{p} \nabla \cdot \left( u_i(\mathbf{x}, t) \int_{B_{r_i}} \alpha_{i,j}(\mathbf{y}) \mathbf{K}_{i,j}(\mathbf{x}, \mathbf{y}) f_{i,j}(u(\mathbf{y}, t)) \, d\mathbf{y} \right) + F(u_1, \dots, u_i, \dots, u_p), \quad (8) $$

where $D_i$ is the diffusion coefficient for the $i^{th}$ population, $\alpha_{i,j}$ is the adhesion coefficient between the $i^{th}$ and $j^{th}$ population, $r_i$ the sensing radius of $i^{th}$ cell population, $\mathbf{K}_{i,j}$ and $f_{i,j}$ describe the induced direction of movement and its magnitude. Moreover, function $F(u_1, \dots, u_i, \dots, u_p)$ describes the proliferation and decay of these cells (sometimes as a result of interactions between cells $u_i$ and $u_j$). For different examples of such models (8), that incorporate different functions (and interaction kernels), we refer the reader to [27].

Before concluding this discussion on models for cell movement we need to emphasise that in many cases, these non-local models for cell-cell interactions are coupled with ODEs describing, for example, the remodelling of the extracellular matrix (ECM). The ECM does not move in space; it can only be produced and degraded by different cells. Therefore, these models for cell-ECM interactions take the form:

$$ \frac{\partial u_i}{\partial t} = D_i \Delta u_i - \sum_{j=1}^{p} \nabla \cdot \left( u_i(\mathbf{x}, t) \int_{B_{r_i}} \alpha_{i,j}(\mathbf{y}) \mathbf{K}_{i,j}(\mathbf{x}, \mathbf{y}) f_{i,j}(u(\mathbf{y}, t); e(t)) \, d\mathbf{y} \right) + F(u_1, \dots, u_i, \dots, u_p), \quad (9a) $$

$$ \frac{de}{dt} = p_e(u_i(\mathbf{x}, t)) - d_e(u_i(\mathbf{x}, t))e(t). \quad (9b) $$

In the equations above, the function $f_{i,j}$ can depend also on the density of the ECM $e$; and $p_e$ and $d_e$ functions describe the secretion and degradation of the ECM (by different cells in the system). All other terms are as described for model (8).

## 2.2 Non-local models for ecology

Animals, like cells, can move collectively and form spatial and spatio-temporal social aggregations. The variety of these aggregations (e.g., swarms of insects, flocks of birds, schools of fish or herds of ungulates) has attracted the attention of both biologists and mathematicians who tried to identify common mechanisms behind them. Different empirical and mathematical studies have suggested that the three types of social interactions (i.e., repulsion, attraction and alignment) that are involved in the formation of these aggregations, might act on different spatial ranges [1, 52]:
• at very short distances, animals move away from each other, to avoid physical contact [53];

• at a slightly greater distance, animals tend to align their movements (e.g., walking in single file) [54];

• over large distances, animals tend to approach each other (if they are too far apart) in order to maintain group cohesion [53, 55, 56].

The spatial distances over which these interactions occur depend on the communication ranges of animals (similar to cells' sensing radius discussed in the previous section). For example, migrating birds can coordinate their movement direction and speed across distances of more than 250 meters [55], while swarming insects coordinate their movements within 200 mm [57]. Moreover, such aggregations can persist for short periods of time (as for swarms of mosquitoes [58]) or for longer periods of time, such as a season or an entire lifetime (as for wolves [59]). For example, over a wider time window, we notice that animals do not generally use all the available space but confine themselves to an area called home range, which they keep for a season or for their entire life [60]. The presence of these home ranges for different animal species is associated with species-specific processes such as memory (as some animals tend to return to the same area after their migration or hibernation) [61], or stigmergy (as some animals leave physical traces, such as scent marks, when they move across the landscape) [62]. Therefore, the simple mathematical models that we will review next to describe the spatio-temporal dynamics of various ecological aggregations, include not only local/non-local animal-animal interactions but also local/non-local animal-environment interactions.

To present these models for animal movement, let us define $u(\mathbf{x}, t)$ to be the density of an animal population at time $t$ and position $\mathbf{x} \in \mathbb{R}^n$ with $n \in \{1, 2, 3\}$. The spatio-temporal changes in this density can be described by the following generic non-local reaction-diffusion-advection model [27]:

$$ \partial_t u = D\Delta u - \nabla \cdot [u\nabla\nu(\omega_r * f(u))] + F(u), \quad (10) $$

with $D$ the constant diffusion coefficient, $\nu$ the advection coefficient, the product between $\omega_r$ (scalar function) and $f$ having the same meaning as for the taxis-like model (7), and the function $F$ representing the births and deaths of the animal population.

Furthermore, it is quite common to couple the non-local equation with an ODE representing the process of creating memories and/or traces of passages, and their destruction [43]:

$$ \partial_t u = D\Delta u - \nabla \cdot [u\nabla\nu(\omega_r * f(u; m))] + F(u), \quad (11a) $$

$$ \partial_t m = \gamma u - \delta m, \quad (11b) $$

with $m(\mathbf{x}, t)$ the cognitive map, and the constants $\gamma$ and $\delta$ describing the production/destruction of these memories. Model (11) can also be easily adapted to describe
interactions between multiple animal species (as it was done for model (8)):

$$ \partial_t u_i = D_i \Delta u_i - \nabla \cdot [u_i \nabla \sum_{j=0}^{p} \nu_{i,j}(\omega_r * f(u_i; m_i))] \\ + F(u_1, \dots, u_i, \dots, u_p), \qquad (12a) $$

$$ \partial_t m_i = \gamma u_i - \delta m_i. \qquad (12b) $$

Throughout this Section 2.1 we summarised the most general form of models used to describe nonlocal interactions in cell biology and in ecology, to emphasise the similarities of these models. However, we need to emphasise here also an important difference, which plays a particular role in numerical simulations: the spatial domain. While the ecological domain is relatively simple and can be described by a rectangular domain, the biological domain can take very complex geometrical shapes as cells form tissues and organs. The numerical methods derived to simulate the behaviour of these models should be able to accommodate these complex domains (for more on this see Section 6.5). Next, we summarise briefly some of the analytical approaches used to investigate these models (Section 3) followed by a more in-depth discussion of numerical approaches that can be used to simulate the dynamics of these models (Section 4)

## 3 Analytical approaches for non-local models

Since the overall goal of this article is to summarize the numerical approaches, in this Section we only mention (without going into details) certain existing analytical results for non-local models, where the non-locality appears in the advective part (models (2) and (3)). For a recent detailed review of different analytical approaches applied to study such non-local models (with non-localities either in the reaction, diffusion or advection terms), we refer the reader to [27, 42].

Many analytical studies focus on the existence and uniqueness of bounded solutions displayed by these non-local models (2) and (3): from the existence (under certain hypotheses on the initial condition $u_0$ or even on the source function $f$) of positive weak and/or global solutions, [30, 35, 37, 63], to the convergence of these solutions towards solutions of the associated local models (i.e. where $K[u] = \text{constant}$ in Equation (2)) when the sensing radius $r$ of the cell approaches zero [19, 37]. In addition, various studies [19, 27] focused on pattern formation, as these models can produce very rich spatial and spatio-temporal patterns: from different types of stationary states to traveling fronts and pulses. Analytical approaches to investigate pattern formation start with a stability analysis of spatially homogeneous states (where all individuals/particles are uniformly spread over the domain) and/or heterogeneous states (where individuals/particles are aggregated at particular spatial positions). For example, equation (2) could exhibit spatial aggregations only for parameter $\beta$ above a threshold value [19, 27].
Pattern formation in more complex systems of non-local PDEs (with more than 2-3 population variables) in two and even three spatial dimensions is much more difficult to be performed analytically. In this case, numerical bifurcations approaches could help. Nevertheless, continuation algorithms for these non-local PDEs are very few and not standard (due to complexities of these equations), as we will discuss in Section 6.3.

# 4 Numerical approaches for non-local models

Let us now get to the heart of the matter, and discuss the numerical approaches implemented to solve such aggregation-diffusion problems. Although finite differences and finite volumes methods are often used to numerically solve non-local problems, the details of these approaches are rarely made explicit in the literature, for such models. This is why we will begin this section with a brief on the application of finite difference (Section 4.1) and finite volume (Section 4.2) methods to non-local models. Note that since there are numerous such schemes that one could use to discretize these equations, here we will focus only on some simple commonly-used schemes that we also show how to implement numerically. Then, in the section 5 we will look at resolution of these non-local models by the finite element method, for which there are very few studies in the literature.

*Remark 3* In this article, we have chosen to deal only with “classical” numerical approaches, namely finite difference, finite volume and finite element methods. Nevertheless, there are other numerical approaches that work well for non-local advection-dominated models, such as the spectral methods. We refer the reader to the papers of A. Gerisch [25, 64] and A. Buttenschoen [65] in which this method is well described and applied to non-local models. We prefer to not present spectral methods in detail in this study, since we want to focus mainly on the very few applications of finite element approaches for these non-local models, which further allows us to discuss the current open problems in this field.

To present these different numerical methods, we focus on the following simple and generic non-local model (13), where the non-locality is included in the advective part. This model, from the paper of Armstrong, Painter and Sherratt [19], is a special case of the model (7) described previously in Section 2.1 and of model (10) described previously in Section 2.2:

$$ \frac{\partial u}{\partial t} = D\Delta u - \nabla \cdot (uK[u]), \qquad (13) $$

with the non-local term $K[u]$ taking the following form

$$ K[u](\mathbf{x}, t) = \alpha \int_{B_r(0)} g(u(\mathbf{x} + \mathbf{y}, t))\omega(||\mathbf{y}||_2) \frac{\mathbf{y}}{||\mathbf{y}||_2} d\mathbf{y}, \qquad (14) $$

with $\mathbf{x} \in \mathbb{R}^n$ with $n \in \{1, 2, 3\}$, $r$ the cell/animal sensing radius, $\alpha$ the coefficient describing the strength of interactions between cells/animals, $\omega$ the kernel describing the magnitude of these forces according to the distance of the cell/animal located in
$\mathbf{y} \in \mathbb{R}^n$ relative to $\mathbf{x}$. Finally, $\mathbf{y} \mapsto \frac{\mathbf{y}}{\|\mathbf{y}\|_2}$ is the vector indicating the direction of these forces.

Since the goal of this review is to showcase some numerical approaches for these non-local models, throughout this section we consider the following simple functions (which do not incorporate any of the biological realism discussed in Section 2.1): $g(u) = u$ and $\omega(\cdot) = 1$. Note that usually $\int_{B_r} \omega = 2$, but in our case $\omega$ is a constant so we let $\omega = 1$ and integrate the constant $\frac{1}{2|B_r|}$ in $\alpha$. The simplicity of these choices will allow us to present clearly the numerical schemes (without being sidetracked by complicated mathematical expressions that could arise from the discretisation of more realistic but highly nonlinear functions $g$ and $\omega$). Also to keep the model as simple as possible, we will consider only the non-local PDE model (13)-(14), and ignore the coupled PDE-ODE models discussed above in equations (9) and (12). These simplifications on model equations (which, as mentioned above, might not be very biologically realistic) also allow us to compare, for the purpose of this review, our numerical simulations obtained with different numerical schemes (implemented using an in-house Python code as well as the open-source FEniCS platform) with the numerical simulations in [19].

Finally, the non-local model (13) and (14) has to be completed with boundary conditions, which need to consider the impact of the non-local operator [28]. Examples of biologically-realistic boundary conditions are: (i) periodic boundary conditions (where cells/animals leave the domain through one end and re-enter it through the opposite end), which are usually employed to describe “infinite-like” domains where the boundaries do not affect the dynamics inside the domain; or (ii) zero-flux boundary conditions (where cells/animals cannot cross the domain boundaries). For the simulations performed throughout this Section, we will focus only on periodic boundary conditions. Mathematically, by working on a 1D interval of length 2L (or a 2D square or a 3D cube with sides of lengths 2L), this amounts to looking for a solution $u$ that is 2L-periodic. Hence, the non-local operator $K[u]$ is also 2L-periodic by definition (14).

*Remark 4* For the description of the different numerical schemes in sub-sections 4.1, 4.2 and 5.1 below, we focus only on the following 1D model (although we present also simulations for 2D models):

$$ \partial_t u = -D \partial_{xx}^2 u + \partial_x(uK[u]), \quad (15) $$

with

$$ K[u](x, t) = \alpha \int_{-r}^{r} u(x + y, t) \operatorname{sgn}(y) \, dy, \quad (16) $$

In this $K[u]$ we have considered $\omega(|y|) = 1$, $g(u) = u$, and $\alpha$ constant. In addition, we work with periodic conditions at the domain boundaries.

## 4.1 Finite differences

Before going into the details of the discretization by finite difference of non-local models, let us cite a few articles in which this method is used to carry out simulations: [20, 31]. To describe the finite differences method, we place ourselves in the following framework: the domain is $\Omega = [-L, L]$ discretized in $N$ sub-intervals $[x_i, x_{i+1}]$ with $i = 0, \dots, N-1$. Hence, we note $h = 2L/N$ the mesh step and $\tau$ the
time step such as $t_n = n\tau$.

For $n \in \mathbb{N}$, we assume to have computed $(u_i^n)_{i=0,...,N}$ an approximation of $(u(x_i, t_n))_{i=0,...,N}$, and $(K_i^n)_{i=0,...,N}$ an approximation of $(K[u](x_i, t_n))_{i=0,...,N}$. At the next time step we obtain an approximation of $u(x_i, t_{n+1})$ denoted $u_i^{n+1}$ by using the following discretization of equation (15) : we consider an Euler finite difference scheme for the approximation of the time derivative

$$ \frac{\partial}{\partial t} u(x_i, t_n) \approx \frac{u_i^{n+1} - u_i^n}{\tau}, $$

combined with an implicit centered difference scheme for the second-order spatial derivative (for stability reasons [23])

$$ \frac{\partial^2}{\partial x^2} u(x_i, t_{n+1}) \approx \frac{u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}}{h^2}, $$

and a mix implicit-explicit centered difference scheme for the transport term:

$$ \frac{\partial}{\partial x} [u(x_i, t_{n+1}) K[u](x_i, t_n)] \approx \frac{u_{i+1}^{n+1} K_{i+1}^n - u_{i-1}^{n+1} K_{i-1}^n}{2h}, $$

We could also consider the explicit centered difference scheme for the transport term:

$$ \frac{\partial}{\partial x} [u(x_i, t_n) K[u](x_i, t_n)] \approx \frac{u_{i+1}^n K_{i+1}^n - u_{i-1}^n K_{i-1}^n}{2h}, $$

but as we will see in Section 5.2.1, explicit scheme can leads to numerical oscillations, depending on the value of the time step and the parameters of the equation. In contrast, the mix explicit-implicit scheme for the transport term makes these oscillations disappear.

In the previous equations, $K_i^n$ is an approximation of the non-local operator evaluate at point $x_i$ and time $t_n$. This term could be approximated using the following trapezoid quadrature method (or any other quadrature methods [66]):

$$
\begin{align}
K[u](x_i, t_n) &= \alpha \int_{x_i-r}^{x_i+r} u(y, t_n) \operatorname{sgn}(y - x_i) dy \nonumber \\
&= \alpha \sum_{k=0}^{N_r-1} \left( \int_{x_{i+k}}^{x_{i+k+1}} u(y, t_n) dy - \int_{x_{i-k-1}}^{x_{i-k}} u(y, t_n) dy \right) \nonumber \\
&\approx \frac{\alpha h}{2} \sum_{k=0}^{N_r-1} \left( u_{i+k+1}^n + u_{i+k}^n - u_{i-k}^n - u_{i-k-1}^n \right) \nonumber \\
&= \alpha h \left( \frac{u_{i+N_r}^n - u_{i-N_r}^n}{2} + \sum_{k=1}^{N_r-1} (u_{i+k}^n - u_{i-k}^n) \right) = K_i^n, \tag{17}
\end{align}
$$
with $N_r = \lfloor r/\hbar \rfloor$ the number of sub-intervals $[x_i, x_{i+1}]$ included in “the sensitivity interval” $[0, r]$. Moreover, thanks to the periodicity of $u$, $\forall j < 0$, $u_j = u_{j+N}$, and similarly $\forall j > N$, $u_j = u_{j-N}$.

Thus, $\forall i = 1, ..., N - 2$, we have the following scheme:

$$u_i^{n+1}(1 + 2\beta) + u_{i+1}^{n+1}(-\beta + \gamma K_{i+1}^n) + u_{i-1}^{n+1}(-\beta - \gamma K_{i-1}^n) = u_i^n.$$

By considering also periodic conditions at the domain boundary (i.e., $i = 0$ and $i = N - 1$) leads to the following scheme

$$u_0^{n+1}(1 + 2\beta) + u_1^{n+1}(-\beta + \gamma K_1^n) + u_{N-1}^{n+1}(-\beta - \gamma K_{N-1}^n) = u_0^n, \quad (18a)$$

$$u_i^{n+1}(1 + 2\beta) + u_{i+1}^{n+1}(-\beta + \gamma K_{i+1}^n) + u_{i-1}^{n+1}(-\beta - \gamma K_{i-1}^n) = u_i^n, \quad \forall i = 1, \dots, N - 2, \quad (18b)$$

$$u_{N-1}^{n+1}(1 + 2\beta) + u_0^{n+1}(-\beta + \gamma K_N^n) + u_{N-2}^{n+1}(-\beta - \gamma K_{N-2}^n) = u_{N-1}^n, \quad (18c)$$

$$u_N^{n+1} = u_0^{n+1}, \quad (18d)$$

$$\text{with } \beta = \frac{D\tau}{h^2}, \gamma = \frac{\tau}{2h}.$$

In Section 4.3 we show the numerical implementation of this scheme for model (15)-(16).

## 4.2 Finite volumes

We now present a discretization of the non-local equation (13) using a finite volume approach, as these methods are widely used for such non-local equations [19, 35, 38, 67]. We start by rewriting equation (15) as:

$$\partial_t u + \partial_x J(x) = 0, \quad (19)$$

with the flux $J = J_d + J_a$ defined previously (5)-(6) in Section 2.1. To describe the finite volumes method for the non-local problem (15), we start by considering the same domain as in the previous section: $\Omega = [-L, L]$ that we discretize in $N$ sub-intervals $C_i = [x_{i-1/2}, x_{i+1/2}]$ with $x_{i\pm 1/2} = -L + (i \pm 1/2)h$, for $i = 1, ..., N$. Hence, we note $h = 2L/N$ the mesh step and $\tau$ the time step such that $t_n = n\tau$.

Hence, by integrating equation (19) on the interval $C_i$ (or a similar volume for higher spatial dimensions), we have:

$$\int_{C_i} \partial_t u(x, t) \, dx + \int_{C_i} \partial_x J(x, t) \, dx = 0.$$

We can use the same Euler finite difference scheme as in Section 4.1, for the approximation of the time derivative, to obtain:

$$h \frac{u_i^{n+1} - u_i^n}{\tau} + J(x_{i+1/2}, t_{n+1}) - J(x_{i-1/2}, t_{n+1}) = 0,$$
with $u_i^n$ be the average of $u(x, t_n)$ on $C_i$. As in the previous Section 4.1, do not use a completely implicit scheme. Indeed, the diffusive flux is treated implicitly but for the transport flux we use a mix explicit-implicit discretization approach, to avoid numerical oscillations (see Section 5.2.1). Hence, $J_{i\pm 1/2}^n$ is an approximation of the following flux:

$$J_{i\pm 1/2}^{n+1} \approx -D \frac{\partial u(x_{i\pm 1/2}, t_{n+1})}{\partial x} + u(x_{i\pm 1/2}, t_{n+1}) K[u](x_{i\pm 1/2}, t_n).$$

So we have, $\forall i = 1, ..., N - 1$

$$u_i^{n+1} + \frac{\tau}{h} (J_{i+1/2}^{n+1} - J_{i-1/2}^{n+1}) = u_i^n. \quad (20)$$

Then, we approximate the first spatial derivative with a central difference scheme (centered at $x_{i+1/2}$)

$$\frac{\partial u(x_{i+1/2}, t_{n+1})}{\partial x} \approx \frac{u_{i+1}^{n+1} - u_i^{n+1}}{h},$$

and the value of the non-local term $K[u]$ and of $u$ at point $x_{i\pm 1/2}$ by

$$K[u](x_{i\pm 1/2}, t_n) \approx \frac{K_{i\pm 1}^n + K_i^n}{2}, \quad u(x_{i\pm 1/2}, t_{n+1}) \approx \frac{u_{i\pm 1}^{n+1} + u_i^{n+1}}{2},$$

with $K_i^n$ an approximation of $K[u](x_i, t_n)$ calculated as in the previous Section 4.1. Hence, the scheme (20) becomes

$$u_i^{n+1} + \frac{\tau}{h} \left( -D \frac{u_{i+1}^{n+1} - u_i^{n+1}}{h} + \frac{u_{i+1}^{n+1} + u_i^{n+1}}{2} \frac{K_{i+1}^n + K_i^n}{2} + D \frac{u_i^{n+1} - u_{i-1}^{n+1}}{h} - \frac{u_i^{n+1} + u_{i-1}^{n+1}}{2} \frac{K_i^n + K_{i-1}^n}{2} \right) = u_i^n,$$

and so

$$u_i^{n+1} \left( 1 + 2\beta + \gamma \frac{K_{i+1}^n - K_{i-1}^n}{2} \right) + u_{i+1}^{n+1} \left( -\beta + \gamma \frac{K_{i+1}^n + K_i^n}{2} \right) + u_{i-1}^{n+1} \left( -\beta - \gamma \frac{K_i^n + K_{i-1}^n}{2} \right) = u_i^n, \quad \forall i = 1, ..., N-2,$$

with $\beta = \frac{D\tau}{h^2}$ and $\gamma = \frac{\tau}{2h}$.

Moreover, by adding periodic boundaries conditions, we obtain the following system:

$$u_0^{n+1} \left( 1 + 2\beta + \gamma \frac{K_1^n - K_{N-1}^n}{2} \right) + u_1^{n+1} \left( -\beta + \gamma \frac{K_1^n + K_0^n}{2} \right) + u_{N-1}^{n+1} \left( -\beta - \gamma \frac{K_0^n + K_{N-1}^n}{2} \right) = u_0^n, \quad (21a)$$
$$ u_i^{n+1} \left( 1 + 2\beta + \gamma \frac{K_{i+1}^n - K_{i-1}^n}{2} \right) + u_{i+1}^{n+1} \left( -\beta + \gamma \frac{K_{i+1}^n + K_i^n}{2} \right) \\ + u_{i-1}^{n+1} \left( -\beta - \gamma \frac{K_i^n + K_{i-1}^n}{2} \right) = u_i^n, \quad \forall i = 1, \dots, N-2, \qquad (21b) $$

$$ u_{N-1}^{n+1} \left( 1 + 2\beta + \gamma \frac{K_0^n - K_{N-2}^n}{2} \right) + u_0^{n+1} \left( -\beta + \gamma \frac{K_0^n + K_{N-1}^n}{2} \right) \\ + u_{N-2}^{n+1} \left( -\beta - \gamma \frac{K_{N-1}^n + K_{N-2}^n}{2} \right) = u_{N-1}^n, \qquad (21c) $$

$$ u_N^{n+1} = u_0^{n+1}. \qquad (21d) $$

In Section 4.3 we will see an example of numerical simulations obtained with the implementation of this finite volume scheme for model (15)-(16).

## 4.3 Numerical simulation

To verify that the two schemes finite differences (18) and finite volumes (21) give the same result, we implement the resolution of these two schemes in Python (see Figure 2). We take as initial condition

$$ u_0(x) = 1 + \text{random}(0, 10^{-2}), $$

and all parameters of the equation (15) and of the discretization are the same in the two different schemes, see Table 1.

<table><caption>**Table 1:** Description of model (15)-(16) parameters with their values used for the numerical simulations presented in Figures 2, 3, and 4.</caption><thead><tr><td>Parameters</td><td>Value</td><td>Description</td></tr></thead><tbody><tr><td>D</td><td>1</td><td>Diffusion coefficient</td></tr><tr><td>α</td><td>10</td><td>Adhesion coefficient</td></tr><tr><td>r</td><td>1</td><td>Sensitivity radius</td></tr><tr><td>L</td><td>10</td><td>Ω = [-L, L]</td></tr><tr><td>N</td><td>1000</td><td>Number of intervals</td></tr><tr><td>τ</td><td>0.01</td><td>Step time</td></tr></tbody></table>

# 5 A new numerical approach for non-local models: finite elements method

Although various numerical approaches have already been considered for different non-local models [19, 20, 31, 35, 38, 67], very few deal with the finite element method due to the mathematical technicalities caused by the non-local terms [5]. This is why we present here in detail the implementation of this method for the non-local aggregation-diffusion equation (13) with periodic boundaries conditions seen previously. We start
<figure><img src="image_3.png" /><figcaption>Fig. 2 Numerical resolution of equation (15) using numerical schemes (18) (red dashed curve) and (21) (blue continuous curve). We show here the solution curves at times $t = 0, t = 0.01, t = 0.5, t = 0.75, t = 1$ and $t = 10$. The parameter values are as in Table 1.</figcaption></figure>

with the technical details of this numerical method; then we see the different numerical simulations carried out in both Python (in-house code) and FEniCS (open access computing platform) in the one-dimensional (1D) and two-dimensional (2D) cases, for one or two cell populations.
## 5.1 Finite Elements discretization

To simplify the discretization by finite elements, we place ourselves in the 1D case (15). Let $\Omega = [-L, L]$ denote the 1D spatial domain and $I = [0, T]$ denote the temporal interval. Thus, the problem in which we are interested is the equation (15) where $K[u]$ is the non-local form (16), coupled with periodic conditions at the edges:

$$ egin{cases} 
\frac{\partial u(x,t)}{\partial t} = D \frac{\partial^2 u(x,t)}{\partial x^2} - \frac{\partial}{\partial x}(u(x,t)K[u](x,t)), \\ 
u(x,t=0) = u_0(x), \\ 
u \text{ is } 2L\text{-periodic.} 
\end{cases} 
\quad (22) $$

For spatial discretization, we use the finite element method. Thus, we multiply the first line of the system (22) by a test function $v$ which belongs to the function space $H(\Omega)$. By defining:

$$ H(\Omega) = \{v \in H^1(\Omega), v \text{ is } 2L\text{-periodic}\}. $$

Next, we integrate over the spatial domain the first line of system (22) multiplied by $v$ and apply Green's formula. Thus, we obtain:

$$ \int_{-L}^{L} \frac{\partial u(x,t)}{\partial t} v(x) dx = \int_{-L}^{L} \left[ D \frac{\partial^2 u(x,t)}{\partial x^2} - \frac{\partial}{\partial x} (u(x,t)K[u](x,t)) \right] v(x) dx, $$

and so,

$$ \int_{-L}^{L} \frac{\partial u(x,t)}{\partial t} v(x) dx = D \left( \left[ \frac{\partial u(x,t)}{\partial x} v(x,t) \right]_{-L}^{L} - \int_{-L}^{L} \frac{\partial u(x,t)}{\partial x} \frac{\partial v(x)}{\partial x} dx \right) - \left[ u(x,t)K[u](x,t) v(x,t) \right]_{-L}^{L} + \int_{-L}^{L} u(x,t)K[u](x,t) \frac{\partial v(x)}{\partial x} dx. $$

But since $u$ and $v$ are $2L$-periodic, the terms in square brackets, in the second equality, are zero. Hence, the variational formulation of the problem (22) is:

Find $u \in H(\Omega)$ such that, for all $v \in H(\Omega)$,

$$ \int_{-L}^{L} \frac{\partial u(x,t)}{\partial t} v(x) dx + D \int_{-L}^{L} \frac{\partial u(x,t)}{\partial x} \frac{\partial v(x)}{\partial x} dx - \int_{-L}^{L} u(x,t) K[u](x,t) \frac{\partial v(x)}{\partial x} dx = 0. \quad (23a) $$

Now, for the time discretization, we use a standard Euler scheme (as in Section 4.1 and 4.2). Furthermore, we choose to treat the diffusion term implicitly, while we treat the advective part using a mixed explicit-implicit discretisation approach. These choices were made in order to avoid stability problems of the numerical scheme. We
will see later (Figure 4 in Section 5.2.1) a numerical illustration of this instability which will justify our choice to not treat the advective part in a purely explicit manner. Thus, the weak form (23) becomes:

Find $u^{n+1} \in H(\Omega)$ such that, for all $v \in H(\Omega)$,

$$a_h(u^{n+1}, v; u^n) = l_h(u^n; v), \qquad (24)$$

with

$$a_h(u^{n+1}, v; u^n) = \int_{\Omega} u^{n+1} v \, dx + D\tau \int_{\Omega} \frac{\partial u^{n+1}}{\partial x} \frac{\partial v}{\partial x} \, dx - \tau \int_{\Omega} u^{n+1} K(u^n) \frac{\partial v}{\partial x} \, dx, \\ l_h(u^n; v) = \int_{\Omega} u^n v \, dx.$$

Before presenting the problem in its discrete form, we discretize the domain $\Omega$ into $N$ intervals, $[x_i, x_{i+1}] = [x_0 + ih, x_0 + (i+1)h]$ for $i = 0, ..., N-1$, of length $h = \frac{2L}{N}$ and $\tau$ the time step such as $t_n = n\tau$. Now, by defining the discrete space of finite elements $\mathbb{P}_1 := \{P(x) = ax + b | a, b \in \mathbb{R}\}$:

$$V_h = \{v \in C^0(\Omega), v|_{[x_j, x_{j+1}]} \in \mathbb{P}_1[X], \forall j = 0, ..., N-1; v \text{ is } 2L\text{-periodic}\},$$

we can therefore write the discrete variational formulation of the problem (22):

$$\left\{ \begin{array}{l} \text{Find } u_h^{n+1} \in V_h \text{ such that, } \forall v_h \in V_h, \\ a(u_h^{n+1}, v_h; u_h^n) = l(u_h^n; v_h). \end{array} \right. \qquad (25)$$

But solving this problem (25) amounts to solve the following problem (see Appendix A for calculation details):

$$u_0^{n+1} \left( \frac{2}{3} + 2\beta + 2\gamma(K_1^n - K_{N-1}^n) \right) + u_1^{n+1} \left( \frac{1}{6} - \beta + 2\gamma(2K_0^n + K_1^n) \right) \\ + u_{N-1}^{n+1} \left( \frac{1}{6} - \beta - 2\gamma(2K_0^n + K_{N-1}^n) \right) = 4u_0^n + u_1^n + u_{N-1}^n, \qquad (26a)$$

$$u_i^{n+1} \left( \frac{2}{3} + 2\beta + 2\gamma(K_{i+1}^n - K_{i-1}^n) \right) + u_{i+1}^{n+1} \left( \frac{1}{6} - \beta + 2\gamma(2K_i^n + K_{i+1}^n) \right) \\ + u_{i-1}^{n+1} \left( \frac{1}{6} - \beta - 2\gamma(2K_i^n + K_{i-1}^n) \right) = 4u_i^n + u_{i+1}^n + u_{i-1}^n, \forall i = 1, \dots, N-2, \qquad (26b)$$

$$u_{N-1}^{n+1} \left( \frac{2}{3} + 2\beta + 2\gamma(K_0^n - K_{N-2}^n) \right) + u_0^{n+1} \left( \frac{1}{6} - \beta + 2\gamma(2K_{N-1}^n + K_0^n) \right) \\ + u_{N-2}^{n+1} \left( \frac{1}{6} - \beta - 2\gamma(2K_{N-1}^n + K_{N-2}^n) \right) = 4u_{N-1}^n + u_0^n + u_{N-2}^n, \qquad (26c)$$
$$u_{N}^{n+1} = u_{0}^{n+1}, \qquad (26d)$$

with $\beta = \frac{D\tau}{2h^2}$ and $\gamma = \frac{\tau}{2h}$. In the previous system, we introduce $K_j^n$ which is the following approximation of the non-local operator $K[u](x_j, t_n)$ (see also Appendix A for calculation details):

$$K_j^n := K(u_h)(x_j, t_n) \approx -\alpha h \mathcal{F}^{-1}(\mathcal{F}(U_h^n)\mathcal{F}(\omega))(x_j, t_n),$$

where $\mathcal{F}$ (resp. $\mathcal{F}^{-1}$) is the Fourier transform (resp. inverse Fourier transform).

## 5.2 Numerical simulations

Now that we have described the discretisation with finite element method for the non-local problem (13), we can move on to the different numerical simulations carried out with an in-house Python code (which does not use finite element libraries) and with a FEniCS code (which uses already-implemented finite element libraries).

### 5.2.1 One population model

*Uni-dimensional case.*

For the 1D case (15), the adhesive flow $K[u]$ take the form (16). In Figure 3 we use the FEniCS software to numerically solve the discrete system(25). We also solved (25) using an in-house Python code ran with the same parameter values (see Table 1) and the results were the same (not shown here). For the initial conditions, see Section 4.3. In the simulation 3, we can see that amplitude of the solution become stable for large time.

*Comparison between explicit scheme and mix explicit-implicit scheme, for the advective part.*

This paragraph aims to justify our choice of choosing a mix explicit-implicit scheme to treat the advection part of the equation (15). Indeed, we have carried out two numerical simulations with these two schemes, for the same parameters values (see Table 1). We notice that for the explicit scheme for the advective part, i.e. the following scheme:

$$\int_{\Omega} u^{n+1} v \, dx + D\tau \int_{\Omega} \frac{\partial u^{n+1}}{\partial x} \frac{\partial v}{\partial x} \, dx = \tau \int_{\Omega} u^n K(u^n) \frac{\partial v}{\partial x} \, dx + \int_{\Omega} u^n v \, dx, \quad (27)$$

some numerical oscillations appear if $\tau$ is not very small, whereas those oscillations disappear when we use the mix explicit-implicit scheme (24) for the advective part, see Figure 4. For the initial conditions we have chosen:

$$u_0(x) = e^{-x^2} + e^{-(x-L)^2} + e^{-(x+L)^2} + e^{-(x-L/2)^2} + e^{-(x+L/2)^2}.$$

Thus, instead of using the explicit scheme and having to refine the time step $\tau$ to make these oscillations disappear, we choose to use the explicit-implicit mix scheme.
<figure><img src="image_4.png" /><figcaption>Fig. 3 Numerical solution of the equation (15) using FEM implemented in FEniCS, at different time. Parameter values are as in Table 1.</figcaption></figure>

These simulations illustrate our choice of scheme for the finite element part, but these results would have been similar for the finite difference and finite volume schemes.
<figure><img src="image_5.png" /><figcaption>Fig. 4 Numerical simulations of model (15)-(16), using the parameter values in Table 1 (except $\alpha = 15$) obtained with : (left column) the explicit scheme (27), (right column) the mix explicit-implicit scheme.</figcaption></figure>

### Bi-dimensional case.

Next we perform simulations with the finite element method for the 2D version of model (13) using FEniCS software (see Table 2 for parameters values). As for the unidimensional case, we also perform simulations using an in-house Python code ran with the same parameter values (see Table 2) and the results were the same (not shown here). For the initial conditions we have chosen

$$u_0(x, y) = 0.2 + \text{random}(0, 10^{-2}).$$

Moreover, in the non-linear term $K[u]$, the linear form $g$ and the kernel $\omega$ are the same as for the 1D case (i.e. $g(u) = u$, $\omega(||y||_2) = \frac{1}{2|B_r|}$ and $f(u) = 0$). The results of these simulations are presented in Figure 5. We see that the small 2D perturbations give
rise to spatial aggregations with finite densities, which are stable in time. Moreover, the solution is positive. Therefore, the numerical scheme behaves well.

<table><caption>**Table 2:** Description of model parameters with their values used for the 2D numerical simulations presented in Figure 5.</caption><thead><tr><th>Parameters</th><th>Value</th><th>Description</th></tr></thead><tbody><tr><td><i>D</i></td><td>1</td><td>Diffusion coefficient</td></tr><tr><td><i>α</i></td><td>10</td><td>Adhesion coefficient</td></tr><tr><td><i>r</i></td><td>1</td><td>Sensitivity radius</td></tr><tr><td><i>L</i></td><td>2.5</td><td>Ω = [-<i>L</i>, <i>L</i>] × [-<i>L</i>, <i>L</i>]</td></tr><tr><td><i>N</i></td><td>80</td><td>Number of intervals</td></tr><tr><td><i>τ</i></td><td>0.1</td><td>Step time</td></tr></tbody></table>

### 5.2.2 Two populations model

We now consider the modelling of two cell populations, *u* and *v*, interacting through adhesion. In this new model, we must take into account the different attraction and adhesion relationships between cells. For each cell type, we assume that there are two adhesion forces: the first representing the adhesion force between cells of the same population, and the second representing the adhesion force between two cells of different populations. Thus, the non-local model for two populations of cells is as follows:

$$ egin{cases} \frac{\partial u}{\partial t} = D\Delta u - \nabla \cdot (uK_u(u, v)), \\ \frac{\partial v}{\partial t} = D\Delta v - \nabla \cdot (vK_v(u, v)), \end{cases} \quad (28) $$

where *D* is the diffusion and, $K_u$ and $K_v$ are the following non-local forms:

$$ K_u(\mathbf{x}, t) = S_u \int_{B_r(0)} g_{uu}(u(\mathbf{x} + \mathbf{y}, t), v(\mathbf{x} + \mathbf{y}, t))\omega(\|\mathbf{y}\|_2) \frac{\mathbf{y}}{\|\mathbf{y}\|_2} dy \\ + C \int_{B_r(0)} g_{uv}(u(\mathbf{x} + \mathbf{y}, t), v(\mathbf{x} + \mathbf{y}, t))\omega(\|\mathbf{y}\|_2) \frac{\mathbf{y}}{\|\mathbf{y}\|_2} dy, \quad (29) $$

and,

$$ K_v(\mathbf{x}, t) = S_v \int_{B_r(0)} g_{vv}(u(\mathbf{x} + \mathbf{y}, t), v(\mathbf{x} + \mathbf{y}, t))\omega(\|\mathbf{y}\|_2) \frac{\mathbf{y}}{\|\mathbf{y}\|_2} dy \\ + C \int_{B_r(0)} g_{vu}(u(\mathbf{x} + \mathbf{y}, t), v(\mathbf{x} + \mathbf{y}, t))\omega(\|\mathbf{y}\|_2) \frac{\mathbf{y}}{\|\mathbf{y}\|_2} dy, \quad (30) $$

with $S_u$ (respectively $S_v$) the self-adhesive strength of population *u* (respectively *v*), *C* the cross-adhesive strength between the populations *u* and *v*, $\omega$ describe the magnitude of these forces and $g_i$ for $i \in [uu, uv, vv, vu]$ functions defining the dependence of the adhesive force on cell density. This model, very common in the literature, has already been solved numerically in numerous papers [19, 27].
<figure><img src="image_6.png" /><figcaption>Fig. 5 Numerical solution of the bi-dimensional equation (13) using FEM implemented in FEniCS. We show the solution at six different times: t = 0, t = 1, t = 1.5, t = 2, t = 10, t = 100. The parameter values used for these simulations are listed in Table 2.</figcaption></figure>

### Uni-dimensional case.

We focused on the one dimensional case of the model (28)-(30) coupled with the following initial condition

$$u_0(x) = 0.2 + \text{rand}(0, 10^{-2}),$$
$$v_0(x) = 0.2 + \text{rand}(0, 10^{-2}).$$

Moreover, unlike the simulations of the model for one population of cells, where we took a very simple form for the function $g$, we take here a more complicated form but more realistic from a biological point of view, for the functions $g_i$ (for $i \in [uu, uv, vv, vu]$)

$$g_{uu} = g_{vu} = \begin{cases} u(1 - u - v), & \text{if } u + v < 1, \\ 0, & \text{else,} \end{cases} \qquad (31)$$

and

$$g_{vv} = g_{uv} = \begin{cases} v(1 - u - v), & \text{if } u + v < 1, \\ 0, & \text{else.} \end{cases}, \qquad (32)$$

These two forms of $g_i$ functions, with, $i \in [uu, uv, vv, vu]$, make it possible to include a “population pressure” which reflects the fact that cells are attracted only towards regions where cell density is lower than 1 [19].

<table><caption>**Table 3:** Description of 1D case model (28)- (30) parameters with their values used for the numerical simulations presented in Figure 6.</caption><thead><tr><th>Parameters</th><th>Value</th><th>Description</th></tr></thead><tbody><tr><td><i>D</i></td><td>1</td><td>Diffusion coefficient</td></tr><tr><td><i>S</i><sub>u</sub></td><td>25 or 200</td><td>Self-adhesive strength of population <i>u</i></td></tr><tr><td><i>S</i><sub>v</sub></td><td>7.5 or 15 or 25</td><td>Self-adhesive strength of population <i>v</i></td></tr><tr><td><i>C</i></td><td>0 or 5 or 30 or 50</td><td>Cross-adhesive strength between population <i>u</i> and <i>v</i></td></tr><tr><td><i>r</i></td><td>1</td><td>Sensitivity radius</td></tr><tr><td><i>L</i></td><td>10</td><td>Ω = [-<i>L</i>, <i>L</i>]</td></tr><tr><td><i>N</i></td><td>1000</td><td>Number of intervals</td></tr><tr><td><i>τ</i></td><td>0.01</td><td>Step time</td></tr></tbody></table>

Now that the model for two populations is well defined, we can discretize it using the same approach as for the one-population model (15). Hence, we solved this model (28) using finite elements. Then we implemented its resolution with Fenics. We tested the resolution of this finite element scheme with different parameter values $S_u$, $S_v$ and $C$. As expected, changing the values of the parameters $S_u$, $S_v$ and $C$ strongly impacts the behavior of the two populations over time. Furthermore, these results are in agreement with Steinberg's hypotheses on cell sorting, the Differential Adhesion Hypothesis [68, 69], which states that:

$C = 0 \rightarrow$ Complet sorting;

$C < S_v$ and $C < S_v \rightarrow$ Partial engulfment of one population by the other;

$S_v < C < S_u \rightarrow$ Engulfment of population $u$ by population $v$;
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Numerical solutions of the uni-dimensional equation (28), with FEM implemented in FEniCS, at different time. From left to right: First column ⇒ $S_u = 25$, $S_v = 7.5$ and $C = 0$, second column ⇒ $S_u = 25$, $S_v = 25$ and $C = 5$, third column ⇒ $S_u = 200$, $S_v = 25$ and $C = 50$, and last column ⇒ $S_u = 25$, $S_v = 15$ and $C = 30$. The rest of parameters are listed in Table 3.</figcaption></figure>

$$ \frac{S_u + S_v}{2} \leq C \rightarrow \text{Mixing.} $$

Our simulations (Figure 6) exactly reproduce these hypotheses: first column ⇒ complete sorting, second column ⇒ partial engulfment, third column ⇒ engulfment of population u by population v, and last column ⇒ column mixing.

### Bi-dimensional case.

Analogous to what we have just done for the one-dimensional case, we implemented the finite element resolution, with FEniCS, of the two-dimensional problem with two populations of cells (28) coupled with the following initial condition:

$$ u_0(x, y) = 0.1 + \text{rand}(0, 10^{-2}), $$

$$ v_0(x, y) = 0.1 + \text{rand}(0, 10^{-2}). $$

Moreover, functions $g_i$ for $i \in [uu, uv, vu, vv]$ have the same form as in the previous uni-dimensional case (31) and (32).

<table><caption>**Table 4:** Description of 2D case model (28)- (30) parameters with their values used for the numerical simulations presented in Figure 7.</caption><thead><tr><td>Parameters</td><td>Value</td><td>Description</td></tr></thead></table>
<table><tr><td>D</td><td>1</td><td>Diffusion coefficient</td></tr><tr><td>S<sub>u</sub></td><td>100 or 200 or 250 or 500</td><td>Self-adhesive strength of population u</td></tr><tr><td>S<sub>v</sub></td><td>50 or 100 or 200</td><td>Self-adhesive strength of population v</td></tr><tr><td>C</td><td>0 or 50 or 75 or 150</td><td>Cross-adhesive strength between population u and v</td></tr><tr><td>L</td><td>2.5</td><td>Ω = [-L, L] × [-L, L]</td></tr><tr><td>N</td><td>80</td><td>Number of intervals</td></tr><tr><td>τ</td><td>0.01</td><td>Step time</td></tr></table>

<figure><img src="image_8.png" /><figcaption>Fig. 7 Numerical solutions of the bi-dimensional equation (28), with FEM implemented in FEiNCS, at different time. From left to right: First column ⇒ $S_u = 250$, $S_v = 200$ and $C = 0$, second column ⇒ $S_u = 200$, $S_v = 200$ and $C = 50$, third column ⇒ $S_u = 500$, $S_v = 50$ and $C = 75$, and last column ⇒ $S_u = 100$, $S_v = 100$ and $C = 150$. The rest of parameters are listed in Table 4.</figcaption></figure>

Here are the results obtained for different values of the parameters $S_u$, $S_v$ and $C$. These simulations (Figure 7) once again confirm Steinberg's hypotheses [68, 69]:

* First column: complet sorting because $C = 0$

* Second column: partial engulfment because $C < S_u$ and $C < S_v$

* Third column: engulfment of population u by population v because $S_v < C < S_u$

* Last column: mixing because $\frac{S_u+S_v}{2} < C$

# 6 Open Numerical Problems

There are various open problems associated with the numerical implementation of 1D/2D/3D non-local models. In the following we discuss some of these issues.
## 6.1 In-house codes vs. open-source computing platforms

There are various open-source computational platforms that can be used to implement faster these non-local mathematical models: FEniCS [70], DUNE [71], FreeFEM [72], Chaste [73], etc. Each of these platforms have their own particularities. For this review study, we chose to focus on FEniCS (open source, that already has function implementing the finite element method) and on Python (in-house code, for which we implemented ourselves the finite element method – see Appendix A). Below we summarise some of the advantages and challenges we encountered when discretising and simulating numerically models (13) and (28).

Open source software (FEniCS in our case) generally works well for classic local problems. They are practical and powerful for this type of models. Indeed, these problems have already been, for the most part, implemented and it is therefore possible to find in the documentation [74] or the online forums (https://fenicsproject.discourse.group/), the answers to the questions that one may ask. However, when we want to solve a more complicated problem, we may encounter implementation difficulties. This was the case for the numerical resolution of non-local problems (13) and (28). Indeed, on several points, the implementation of these non-local models posed some problems with the FEniCS software, because we needed to understand what is behind various functions of the main library `dolfin`. Since FEniCS is coded in Python with libraries from C++, understanding the implementation of the objects we are manipulating is not very intuitive. Let us present a concrete problem that we faced: when working with periodic conditions at the boundaries of the domain, the numbering of the degrees of freedom (DOFs) is no longer the same as when we use more classical conditions at the boundaries (e.g. Dirichlet, Neumann). In the one-dimensional case with classic conditions, the DOFs are numbered from left to right (node $x_0$ corresponds to DOF number 0, node $x_1$ to node 1, etc...). However, with periodic conditions, this order is no longer respected, the DOFs are reorganized in such a way as to make the solver more efficient. Thus, when we wanted to assign specific values to a FEniCS function (as was the case for the initial condition $u_0$ or the non-local term $K[u]$ at each time iteration), this posed a problem because the order of the vector coefficients no longer corresponded to the order of numbering of the DOFs. The wrong coefficient was then assigned to the DOF. We therefore had to find a function that reorganized the DOFs in the right order and then perform the assignment. But it was not easy to find as we had difficulty understanding in what order and why FEniCS did not organize the DOFs in logical order.

Unlike these open-source platforms, in-house codes (Python in our case) allow for more freedom and control of the discretization method chosen, and on the actual implementation of these methods. However, this requires having a fairly solid background in the different numerical methods because it is up to the person who implements to calculate by hand the coefficients of the matrices and vectors to be implemented to solve the problem. In the case of finite elements, the programmer must calculate the integrals himself and then implement them efficiently in matrices. These manual calculations can be a source of errors and the choice of implementation
can be crucial if you do not want the program to take too long to run (e.g. use libraries that facilitate calculation with sparse matrices like `scipy sparse`). These problems specific to in-house codes do not concern open-source platforms. Indeed, when using these platforms, it is generally not necessary to have in-depth knowledge of the finite element method. We just need to know the variational formulation of our problem and then implement it. This therefore reduces the risk of errors in calculations and implementation. In addition, the speed of code compilation no longer depends on the programmer's skills, since the main calculations (e.g. calculations of integrals, matrix inversion) are already implemented in the software and the programmer does not have access to them.

Nevertheless, a problem that is common to open-source software and in-house codes is the efficient and inexpensive implementation of the non-local term at each mesh point and each iteration time. This is why throughout this article, we have seen two possible choices to approach the non-local term $K[u]$; see equations (17) and (A2). The first method (17), which uses the trapezoid quadrature rule, is time-consuming. While with the second approximation (A2), is much faster since it uses the discrete Fourier transform, and more precisely an algorithm (Fast Fourier Transform) allowing the discrete Fourier transform of a vector to be calculated very efficiently.

Finally, we managed to exactly reproduce the results obtained with FEniCS (Figures 3 and 5), by implementing the in-house Python code. This therefore shows that FEniCS software works correctly for non-local models.

## 6.2 Numerical oscillations

It is well known that for the following local advection-diffusion equation:

$$ \partial_t u = D \Delta u - \nabla \cdot (\beta u), $$

numerical oscillations appear when we solve this equation numerically with a mesh step $h$ too large. Indeed, if $1 < \frac{|\beta|h}{2D}$ oscillations appear and it is enough to refine the mesh to make them disappear [23]. Thus, we wondered if such oscillations could appear with the non-local model (13). Indeed, by playing with the parameter $\alpha$ of the non-local form $K[u]$, oscillations appear (as we can see in Figure 8) when the advective part dominates the diffusive part (i.e. when $D \ll \alpha$). However, refining the mesh can become very time-consuming (even more for the non-local model than for the local model). Thus, stabilization methods were created for the local model to make these oscillations disappear without reducing the mesh step $h$: Discontinuous Galerkin [23], Streamline Upwind Petrov-Galerkin [23], Continuous Interior Penalty [75], etc. To our knowledge, there is no such methods for non-local models. Thus, we could try, in the near future, to adapt the existing stabilization methods for the local model to the non-local problem. Or even to completely create new stabilizing methods for such non-local models.
<figure><img src="image_9.png" /><figcaption>**Fig. 8** Finite element solution of the equation (15) at time $t = 1.0$ (left) and zoom on one of the interesting parts showing oscillations (right). Equation and mesh parameters are the same as in Table 1, except $\alpha = 13$.</figcaption></figure>

## 6.3 Numerical continuation for non-local advection-dominated models

To have a better understanding of how these aggregation patterns (shown in Figures 3-7) form, persist and sometimes disappear to give rise to other patterns, one can focus on numerical bifurcation theory, and in particular on numerical continuation approaches, which can trace solution branches as model parameters are varied [76, 77]. The general idea behind these algorithms is to use a predictor-corrector numerical scheme combined with singularity detection approaches to track the parametrised solution branches [78]. While there are various software packages for continuation problems arising in Ordinary Differential Equations (e.g., AUTO, XPP, MatCont) [78], this is not the case for Partial Differential Equations (PDEs), where most of the algorithms developed are applied to the specific problems investigated. One of the reasons for this lack of a general continuation software for PDEs is the fact that there are different spatial discretization methods for each class of PDEs (with specific spatial domain characteristics and boundary conditions) [78]. Also, the large number of degrees of freedom for spatial 2D and 3D PDE models also requires specific continuation algorithms [79].

In regard to PDEs, we note that for classical (local) PDE models one can use the MATLAB package `pde2path`, based on the finite element method (FEM) with $P_1$ Lagrangian elements [79]. (We note that the FEM approaches are relatively common for these continuation algorithms due to the algebraic system that results from the discretization of the FEM model: $F(U, \lambda) = 0$, where $\lambda$ is a model parameter and $F$: $\mathbb{R}^n \times \mathbb{R} \to \mathbb{R}$ a nonlinear function [78, 79].) However, more complex PDEs, such as these non-local advection/transport-dominated models that we discussed here, pose various convergence problems [80]. In general, for these non-local PDEs, specific in-house continuation algorithms have been developed for specific parameter ranges [81].
For the future, we expect that the fast development of these non-local aggregation-diffusion models, and the various 2D and 3D patterns that they can generate, will lead to the development of new continuation algorithms for numerical bifurcation studies.

## 6.4 Data availability and Parameter identification

There are numerous studies that look at parameter and function identification for local models using various approaches: from Tikhonov regularisation [82] to Bayesian approaches [83, 84], and references therein. However, there are very few studies that try to identify parameters for non-local models of the type discussed here; among them we mention [85–87]. For example, the study in [86], focused on non-local models for ecological movement of animal populations inside home ranges. The authors first used synthetic data (the equilibrium solutions of the model) to explore parameter inference via maximum likelihood estimation and Stochastic Gradient Descent (SGD) [86]. Then, they used location data for thirteen meerkat populations (from the Kalahari Meerkat Project), to perform a kernel density estimation to use initial population distributions based on this dataset. The study in [87] focused on local and non-local models for collective movement of cells during cancer invasions, and used Tikhonov regularisation approaches to identify the unknown mutation function (describing the transition from a “normal” less invasive cancer cell sub-population to a “mutated” more invasive sub-population). The authors in [87] used synthetic data (to which different levels of noise was added). A very recent study [85] proposed a Partial Inversion algorithm to estimate the non-local kernel (interaction potential) for a simple aggregation-diffusion equation, from synthetic noisy discrete trajectory data.

We emphasise that most of these studies use synthetic data for parameter/function identification. This is because spatial (and mainly spatio-temporal) data is not always available. However, recent advances in imaging techniques will improve this aspect through the generation of numerous datasets quantifying spatial interactions in both ecology [88] and cell biology [89].

## 6.5 Complex domain geometries

The simulations presented in the previous sections were performed on simple rectangular domains. However, the power of finite element method relies in how it deals with complex geometries. And biological and medical applications of advection-dominated equations have considered for now only simple domains, due to the numerical difficulties associated with simulations on more complex geometries. But such complex geometries are very common in medical applications: the structure of tissues, the shape of organs, etc.. For this reason we focused this review on the finite element method: to emphasise the numerical possibilities of simulating such non-local models on more realistic 2D and 3D geometries.

# 7 Summary and discussion

The goal of this review paper was to summarise the current numerical approaches used to discretize and simulate numerically non-local aggregation PDEs, that are
becoming more and more common in the literature as a result of more and more
experimental studies that emphasise the importance of non-local interactions in cell
biology and ecology.

Before discussing the different numerical approaches to solving non-local equations, we started this review by explaining the importance and origin of non-local models in biology and ecology. Then, we briefly reviewed the different existing analytical approaches for these non-local models: existence results, stability analysis, bifurcations, etc... We did not focus more on these analytical aspects as they were recently reviewed in [42].

Afterwards, we started the presentation of different numerical approaches that
could be used to discretise such non-local models. For this part, we focused on a very
simple and generic non-local model, that did not incorporate any of the biological
details discussed in Section 2.1. The reason for using this simple model (for numerical
illustration purposes) was two folded: (i) it allowed us to present clearly the various
discretisation schemes (see equations (18), (21), (26)); (ii) it allowed us to compare
our numerical results (obtained with these different schemes) with the simulations
presented in [19], for which we did not have any numerical details. In regard to these
numerical approaches, first we discussed two classical finite differences and finite vol-
umes schemes, then a less common finite element scheme. We then presented a series
of numerical simulations for different non-local models solved using the finite element
method. We concluded this review with a discussion on open numerical problems
related to these non-local models. In this part we first compared two ways of imple-
menting the resolution of these non-local models via two finite element schemes, one
leading to numerical oscillations. Then, we discussed open problems related to numer-
ical continuation algorithms for the numerical investigation of different bifurcation
patterns exhibited by these non-local models. Finally we discussed briefly open prob-
lems related to numerical approaches that can be used to for parameter identification
for such non-local models.

Through this review of numerical approaches for non-local advection-dominated problems, we aimed to emphasise the opportunities to further develop the finite element approaches to discretise and simulate these non-local models not only in 1D, but mainly in 2D and 3D, and on various geometries that are important in biological and medical applications (e.g., complex domains simulating various tissues and organs [90] on which non-local interactions among cells, molecules, etc., seem to gain interest [91]). By discussing (in Section 6) various open problems related to the numerical implementation of the finite element schemes, stabilisation approaches for these schemes, as well as the need to develop new schemes for numerical continuation methods (to investigate the bifurcation of the complex spatio-temporal patterns that can be exhibited by these models), we aimed to suggest some new numerical research directions that could advance the current research in the field. This review complemented other recent reviews in the literature related to non-local models: from more modeling-focused and biologically-focused reviews [27, 43, 91] to more analytically-focused reviews [28], to reviews focused on numerical approaches for non-local models (spectral methods [28], finite volume methods [31] and finite difference methods [38]).
Our overall goals here were (i) to show that finite element methods could also be used to discretise and simulate non-local advection-dominated models; (ii) to open this field of non-local models in biology to researchers focused on numerical analysis, to generate new advances in the field.

## Acknowledgements

This work (JM, AL, RE) has been supported by EIPHI Graduate School (contract ANR-17-EURE-0002) and by the Bourgogne-Franche-Comté Region. RE and AL work was also partially supported by the ANR grant ANR-21-CE45-0025-01. RE and JM work was also partially supported by the Région Bourgogne-Franche-Comté "Accueil de Nouvelle Équipe de Recherche (ANER) 2022" grant number 2022-7-13828.

## Conflict of interest

All authors have no conflicts of interest.

## Appendix A Calculation details of FEM

Given the discrete variational formulation (25), we can compute the matrix $A_h^n$ and the right-hand side $B_h^n$, such that:

$$A_{i,j}^n = a(\phi_i, \phi_j; u_h^n),$$

$$B_i^n = l(u_h^n; \phi_i),$$

with $\phi_i$ for $i = 1, ..., N - 1$ are the classical basis functions (see Figure A1):

$$\phi_i(x) = \begin{cases} \frac{x-x_{i-1}}{h}, & \text{if } x \in [x_{i-1}, x_i], \\ \frac{x_{i+1}-x}{h}, & \text{if } x \in [x_i, x_{i+1}], \\ 0, & \text{otherwise,} \end{cases} \quad \text{and} \quad \phi_0(x) = \begin{cases} \frac{x-x_{N-1}}{h}, & \text{if } x \in [x_{N-1}, x_N], \\ \frac{x_1-x}{h}, & \text{if } x \in [x_0, x_1], \\ 0, & \text{otherwise.} \end{cases}$$

Thus, to find the solution to the problem (25) at each iteration in time, it suffices to solve the matrix system $A_h U_h^{n+1} = B_h^n$.

To solve the previous system, we have to make several approximations: first we apply some quadrature rules to calculate the integral

$$\int_{\Omega} \phi_i(x) K[u^n](x, t) \frac{\partial \phi_j(x)}{\partial x} dx,$$

because we cannot calculate its exact value. Then, we need to calculate the value of $K(u^n)$ at each point of the mesh to be able to apply the quadrature rule. To approximate the value of the non-local term, we could have used the same approximation as in Section 4.1, but we chose to try to find another method, since this non-local
<figure><img src="image_10.png" /><figcaption>**Fig. A1** The basis functions $\phi_i$ associated to node $x_i$ with $i \in \{0, ..., N-1\}$. The basic function $\phi_0$ is drawn in red since it is this which translates the periodic conditions at domain boundaries.</figcaption></figure>

term seems to be similar to a convolution. Hence, we first rewrite it in the form of a convolution:

$$K[u](x, t) = -\alpha \int_{-L}^{L} u(x - y, t) \tilde{\omega}(y) = -\alpha(u * \tilde{\omega})(x),$$

with

$$\tilde{\omega}(x) = \begin{cases} \operatorname{sgn}(x), & \text{if } x \in [-r, r] \\ 0, & \text{else.} \end{cases} \qquad (A1)$$

Under this form, we can therefore use the Discrete Fourier Transform to calculate the non-local term $K[u]$. First, we can discretize the convolution with the trapezoidal rule:

$$\begin{aligned} K[u](x, t) &\approx -\alpha h \sum_{k=0}^{N-1} \frac{u(x - x_{k+1}, t)\tilde{\omega}(x_{k+1}) + u(x - x_k, t)\tilde{\omega}(x_k)}{2} \\ &= -\alpha h \left( \frac{u(x - x_0, t)\tilde{\omega}(x_0) + u(x - x_N, t)\tilde{\omega}(x_N)}{2} + \sum_{k=1}^{N-1} u(x - x_k, t)\tilde{\omega}(x_k) \right). \end{aligned}$$

Next, thanks to the periodicity of $u$ and $\tilde{\omega}$, we have:

$$\begin{aligned} K[u](x, t) &= -\alpha h \left( u(x - x_0, t)\tilde{\omega}(x_0) + \sum_{k=1}^{N-1} u(x - x_k, t)\tilde{\omega}(x_k) \right) \\ &= -\alpha h \sum_{k=0}^{N-1} u(x - x_k, t)\tilde{\omega}(x_k). \end{aligned}$$

So, at any point of the mesh,

$$K[u](x_j, t) = -\alpha h \sum_{k=0}^{N-1} u(x_j - x_k, t) \tilde{\omega}(x_k)$$
$$
\begin{align*}
&= -\alpha h \sum_{k=0}^{N-1} u(-L + jh - (-L + kh), t) \tilde{\omega}(x_k) \\
&= -\alpha h \sum_{k=0}^{N-1} u(L + x_{j-k}, t) \tilde{\omega}(x_k).
\end{align*}
$$

Hence, to obtain the value of $K[u]$ at any point of the mesh, we have to calculate $u$
at the points $x_i + L$. In practice, thanks to the periodicity of $u$ and $K[u]$, we calculate
$K[u](x_n + L, t) = \alpha h \sum_{k=0}^{N-1} u(x_{n-k}, t)\tilde{\omega}(x_k)$, then we translate all the values from L.
So, $\forall j = 0, \dots, N-1$ we calculate:

$$
(K[u])_j^n = -\alpha h \sum_{k=0}^{N-1} u_{j-k}^n \tilde{\omega}_k,
$$

where $\tilde{\omega}_k = \tilde{\omega}(x_k)$ for $k = 0, \dots, N-1$.

Knowing that $u$ and $\tilde{\omega}$ are $N$-periodic, we use the following DFT:

$$
(\mathcal{F}(u^n))_j = \sum_{k=0}^{N-1} u_k^n e^{-\frac{2\pi i}{N}kj}.
$$

Hence,

$$
\begin{align*}
(\mathcal{F}(K[u]^n))_m &= \sum_{s=0}^{N-1} (K[u])_s^n e^{-\frac{2\pi i}{N}sm} \\
&= -\alpha h \sum_{s=0}^{N-1} \sum_{k=0}^{N-1} u_{s-k}^n \tilde{\omega}_k e^{-\frac{2\pi i}{N}sm} \\
&= -\alpha h \sum_{k=0}^{N-1} e^{-\frac{2\pi i}{N}km} \tilde{\omega}_k \sum_{s=0}^{N-1} e^{-\frac{2\pi i}{N}(s-k)m} u_{s-k}^n \\
&= -\alpha h \mathcal{F}(u^n)_m \mathcal{F}(\tilde{\omega})_m.
\end{align*}
$$

And so,

$$
K_i^n := K(u_h)(x_i, t_n) \approx -\alpha h \mathcal{F}^{-1}(\mathcal{F}(U_h^n)\mathcal{F}(\tilde{\omega})). \quad (A2)
$$

We can now calculate the coefficients of the matrix $A_h^n$, by using Simpson's rule when is necessary, at each iteration in time. Thanks to the form of the basis functions $\phi_i$, we already know which coefficients of the matrix $A_h^n$ will be non-zeros. In fact, we have (see Figure A1):

$$
\left\{
\begin{array}{l}
\text{supp}(\phi_0) = [x_0, x_1] \cup [x_{N-1}, x_0], \\
\text{supp}(\phi_i) = [x_{i-1}, x_{i+1}], \quad \forall i \in \{1, \dots, N-2\}, \\
\text{supp}(\phi_{N-1}) = [x_{N-2}, x_{N-1}] \cup [x_{N-1}, x_0],
\end{array}
\right.
$$
thus the non-zero coefficients are:

$$
\begin{cases}
A_{0,1}^n, & A_{0,0}^n, & A_{0,N-1}^n, \\
A_{i,i+1}^n, & A_{i,i}^n, & A_{i,i-1}^n, \quad \forall i \in \{1, \dots, N-2\}, \\
A_{N-1,0}^n, & A_{N-1,N-1}^n, & A_{N-1,0}^n,
\end{cases}
$$

so the matrix $A_h^n$ is almost tridiagonal. Hence, let's start by calculating the coefficients
of the diagonal. For $i = 1, \dots, N-2$:

$$
\begin{align*}
A_{i,i}^n &= a(\phi_i, \phi_i, u_h^n) \\
&= \int_{x_{i-1}}^{x_i} \left[ (\phi_i(x))^2 + D\tau \left(\frac{\partial \phi_i(x)}{\partial x}\right)^2 - \tau \phi_i(x) K(u_h^n)(x) \frac{\partial \phi_i(x)}{\partial x} \right] dx \\
&\quad + \int_{x_i}^{x_{i+1}} \left[ (\phi_i(x))^2 + D\tau \left(\frac{\partial \phi_i(x)}{\partial x}\right)^2 - \tau \phi_i(x) K(u_h^n)(x) \frac{\partial \phi_i(x)}{\partial x} \right] dx \\
&\approx \frac{1}{h^2} \left( \left[ \frac{x^3}{3} - x^2 x_{i-1} + xx_{i-1}^2 \right]_{x_{i-1}}^{x_i} + D\tau h - \frac{\tau h^2}{6} (2K_i^n + K_{i-1}^n) \right. \\
&\qquad \left. + \left[ \frac{x^3}{3} - x^2 x_{i+1} + xx_{i+1}^2 \right]_{x_i}^{x_{i+1}} + D\tau h + \frac{\tau h^2}{6} (2K_i^n + K_{i+1}^n) \right) \\
&= \frac{2h}{3} + \frac{2D\tau}{h} + \frac{\tau}{6} (K_{i+1}^n - K_{i-1}^n).
\end{align*}
$$

And with a similar calculation,

$$
\begin{align*}
A_{0,0}^n &= \frac{2h}{3} + \frac{2D\tau}{h} + \frac{\tau}{6}(K_1^n - K_{N-1}^n), \\
A_{N-1,N-1}^n &= \frac{2h}{3} + \frac{2D\tau}{h} + \frac{\tau}{6}(K_0^n - K_{N-2}^n).
\end{align*}
$$

Now, we calculate the coefficients of the superior diagonal, for $i = 0, \dots, N-2$:

$$
\begin{align*}
A_{i,i+1}^n &= a(\phi_i, \phi_{i+1}, u_h^n) \\
&= \int_{x_i}^{x_{i+1}} \left[ \phi_i(x)\phi_{i+1}(x) + D\tau \frac{\partial \phi_i(x)}{\partial x} \frac{\partial \phi_{i+1}(x)}{\partial x} - \tau \phi_i(x) K(u_h^n)(x) \frac{\partial \phi_{i+1}(x)}{\partial x} \right] dx \\
&\approx \frac{1}{h^2} \left( \left[ \frac{x^2}{2}(x_{i+1} + x_i) - xx_i x_{i+1} - \frac{x^3}{3} \right]_{x_{i-1}}^{x_i} - D\tau h + \frac{\tau h^2}{6} (2K_i^n + K_{i+1}^n) \right) \\
&= \frac{h}{6} - \frac{D\tau}{h} + \frac{\tau}{6} (2K_i^n + K_{i+1}^n),
\end{align*}
$$

and

$$
A_{N-1,0}^{n} = \frac{h}{6} - \frac{D\tau}{h} + \frac{\tau}{6}(2K_{N-1}^{n} + K_{0}^{n}).
$$

The same for the lower diagonal, for $i = 1, \dots, N-1$:

$$
A_{i,i-1}^n = a(\phi_i, \phi_{i-1}, u_h^n)
$$
$$
\begin{align*}
&= \int_{x_{i-1}}^{x_i} [\phi_i(x)\phi_{i-1}(x) + D\tau \frac{\partial\phi_i(x)}{\partial x} \frac{\partial\phi_{i-1}(x)}{\partial x} - \tau\phi_i(x)K(u_h^n)(x)\frac{\partial\phi_{i-1}(x)}{\partial x}] dx \\
&\approx \frac{1}{h^2} \left( \left[ \frac{x^2}{2}(x_{i-1} + x_i) - xx_i x_{i-1} - \frac{x^3}{3} \right]_{x_i}^{x_{i+1}} - D\tau h - \frac{\tau h^2}{6}(2K_i^n + K_{i-1}^n) \right) \\
&= \frac{h}{6} - \frac{D\tau}{h} - \frac{\tau}{6}(2K_i^n + K_{i-1}^n),
\end{align*}
$$

and

$$
A_{0,N-1}^{n} = \frac{h}{6} - \frac{D\tau}{h} - \frac{\tau}{6}(2K_{0}^{n} + K_{N-1}^{n}).
$$

To complete our systeme, we have to calculate the right hand side vector $B_h^n$. We use,
for any interior node of the mesh, Simpson’s quadrature rule (but we could also have
used the trapezoid method):

$$
\begin{align*}
B_i^n &\approx \int_{\Omega} u_h^n(x) \phi_i(x) \, dx \\
&= \frac{1}{h} \left[ \int_{x_i}^{x_{i+1}} u_h^n(x) (x_{i+1} - x) \, dx + \int_{x_{i-1}}^{x_i} u_h^n(x) (x - x_{i-1}) \, dx \right] \\
&\approx \frac{h}{6} \left[ u_{i-1}^n + 4u_i^n + u_{i+1}^n \right],
\end{align*}
$$

and

$$
\begin{align*}
B_0^n &\approx \frac{h}{6} \left[ u_{N-1}^n + 4u_0^n + u_1^n \right], \\
B_{N-1}^n &\approx \frac{h}{6} \left[ u_{N-2}^n + 4u_{N-1}^n + u_0^n \right].
\end{align*}
$$

Now that we know the form of the matrix $A_h^n$ and the vector $B_h^n$, we rewrite the
matrix system $A_h^n U_h = B_h^n$ as the system (26).

References

[1] Vicsek, T., Zafeiris, A.: Collective motion. Physics Reports **517**(3-4), 71–140 (2012) https://doi.org/10.1016/j.physrep.2012.03.004

[2] Pliny, the Elder: Natural History, Trans. H. Rackham. Cambridge, MA: Harvard University Press (1961)

[3] Friedl, P., Gilmour, D.: Collective cell migration in morphogenesis, regeneration and cancer. Nature Reviews Molecular Cell Biology, 445–457 (2009) https://doi.org/10.1038/nrm2720

[4] Scarpa, E., Mayor, R.: Collective cell migration in development. Journal of Cell Biology **212**(2), 143–155 (2016) https://doi.org/10.1083/jcb.201508047
[5] Adebayo, O.E., Urcun, S., Rolin, G., Bordas, S.P.A., Trucu, D., Eftimie, R.: Mathematical investigation of normal and abnormal wound healing dynamics: local and non-local model. Mathematical Biosciences and Engineering **20**(9), 17446–17498 (2023) https://doi.org/10.3934/mbe.2023776

[6] Hemelrijk, C.K., Hildenbrandt, H.: Schools of fish and flocks of birds: their shape and internal structure by self-organization. Interface focus **2**(6), 726–737 (2012) https://doi.org/10.1098/rsfs.2012.0025

[7] Yan, X., Wang, X., Zhao, Y., Zhu, Q., Yang, L., Li, Z.: Collective decision-making and spatial patterns in orientation of an endemic ungulate on the tibetan plateau. Current Zoology, 004 (2023) https://doi.org/10.1093/cz/zoad004

[8] Bhagavan, H., Muthmann, O., Brockmann, A.: Structural and temporal dynamics of the bee curtain in the open-nesting honey bee species, apis florea. Apidologie **47**, 749–758 (2016) https://doi.org/10.1007/s13592-016-0428-8

[9] McComb, K., Reby, D., Baker, L., Moss, C., Sayialel, S.: Long-distance communication of acoustic cues to social identity in African elephants. Animal Behaviour **65**(2), 317–329 (2003) https://doi.org/10.1006/anbe.2003.2047

[10] Payne, R., Webb, D.: Orientation by means of long range acoustic signaling in baleen whales. Annals of the New York Academy of Sciences **188**, 110–141 (1971) https://doi.org/10.1111/j.1749-6632.1971.tb13093.x

[11] Wilson, E.O.: Chemical communication in the social insects. Science **149**, 1064–1071 (1965) https://doi.org/10.1126/science.149.3688.1064

[12] Mattaini, K.: Introduction to Molecular and Cell Biology. Creative Commons Attribution-NonCommercial 4.0 International License (2020). Chap. 9

[13] Eom, D.S.: Airinemes: thin cellular protrusions mediate long-distance signalling guided by macrophages. Open Biology **10**(8), 200039 (2020) https://doi.org/10.1098/rsob.200039

[14] Reinhart-King, C.A., Dembo, M., Hammer, D.A.: Cell-cell mechanical communication through compliant substrates. Biophysical Journal **95**(12), 6044–6051 (2008) https://doi.org/10.1529/biophysj.107.127662

[15] Judson, O.P.: The rise of the individual-based model in ecology. Trends in Ecology & Evolution **9**(1), 9–14 (1994) https://doi.org/10.1016/0169-5347(94)90225-9

[16] Mogilner, A., Edelstein-Keshet, L., Bent, L., Spiros, A.: Mutual interactions, potentials, and individual distances in a social aggregation. Journal of Mathematical Biology **47**, 353–389 (2003) https://doi.org/10.1007/s00285-003-0209-7
[17] Couzin, I.D., Krause, J., James, R., Ruxton, G.D., Franks, N.R.: Collective memory and spatial sorting in animal groups. Journal of Theoretical Biology **218**, 1–11 (2002) https://doi.org/10.1006/jtbi.2002.3065

[18] Reynolds, C.W.: Flocks, herds and schools: a distributed behavioural model. In: Proc. of the 14th Ann. Conf. on Computer Graphics and Interactive Techniques, 27-21 July, Anaheim, CA, pp. 27–21 (1987). https://doi.org/10.1145/37402.37406

[19] Armstrong, N.J., Painter, K.J., Sherratt, J.A.: A continuum approach to modelling cell–cell adhesion. Journal of Theoretical Biology **243**(1), 98–113 (2006) https://doi.org/10.1016/j.jtbi.2006.05.030

[20] Mogilner, A., Edelstein-Keshet, L.: A non-local model for a swarm. Journal of Mathematical Biology **38**, 534–570 (1999) https://doi.org/10.1007/s002850050158

[21] Pineda, M., Weijer, C.J., Eftimie, R.: Modelling cell movement, cell differentiation, cell sorting and proportion regulation in dictyostelium discoideum aggregations. Journal of Theoretical Biology **370**, 135–150 (2015) https://doi.org/10.1016/j.jtbi.2015.01.042

[22] Topaz, C.M., Bertozzi, A.L., Lewis, M.A.: A nonlocal continuum model for biological aggregation. Bulletin of Mathematical Biology **68**, 1601–1623 (2006) https://doi.org/10.1007/s11538-006-9088-6

[23] Quarteroni, A., Quarteroni, S.: Numerical Models for Differential Problems. Springer (2009). https://doi.org/10.1007/978-88-470-5522-3

[24] Fortunato, I.C., Sunyer, R.: The forces behind directed cell migration. Biophysica **2**, 548–563 (2022) https://doi.org/10.3390/biophysica2040046

[25] Gerisch, A., Chaplain, M.A.J.: Robust numerical methods for taxis–diffusion–reaction systems: applications to biomedical problems. Mathematical and Computer Modelling **43**(1-2), 49–75 (2006) https://doi.org/10.1016/j.mcm.2004.05.016

[26] Berezovskaya, F.S., Karev, G.P.: Bifurcations of travelling waves in population taxis models. Physics-Uspekhi **42**(9), 917 (1999) https://doi.org/10.1070/PU1999v042n09ABEH000564

[27] Painter, K.J., Hillen, T., Potts, J.R.: Biological modelling with nonlocal advection–diffusion equations. Mathematical Models and Methods in Applied Sciences **34**(1), 57–107 (2024) https://doi.org/10.1142/S0218202524400025

[28] Hillen, T., Buttenschön, A.: Non-Local Cell Adhesion Models: Symmetries and Bifurcations in 1-D. Springer (2021). https://doi.org/10.1007/978-3-030-67111-2
[29] Fabrèges, B., Hivert, H., Balc’h, K.L., Martel, S., Delarue, F., Lagoutière, F., Vauchelet, N.: Numerical schemes for the aggregation equation with pointy potentials. ESAIM: Proceedings and Surveys **65**, 384–400 (2019) https://doi.org/10.1051/proc/201965384

[30] Fellner, K., Hughes, B.D.: Solutions of a non-local aggregation equation: Universal bounds, concavity changes, and efficient numerical solutions. Mathematical Methods in the Applied Sciences **43**(8), 5398–5429 (2020) https://doi.org/10.1002/mma.6281

[31] James, F., Vauchelet, N.: Numerical methods for one-dimensional aggregation equations. SIAM Journal on Numerical Analysis **53**(2), 895–916 (2015) https://doi.org/10.1137/140959997

[32] Carrillo, J.A., Fjordholm, U., Solem, S.: A second-order numerical method for the aggregation equations. Mathematics of Computation **90**(327), 103–139 (2021) https://doi.org/10.48550/arXiv.1804.07796

[33] Bernoff, A.J., Topaz, C.M.: Nonlocal aggregation models: A primer of swarm equilibria. SIAM Review **55**(4), 709–747 (2013) https://doi.org/10.1137/130925669

[34] Lorz, A., Lorenzi, T., Hochberg, M.E., Clairambault, J., Perthame, B.: Populational adaptive evolution, chemotherapeutic resistance and multiple anticancer therapies. ESAIM: Mathematical Modelling and Numerical Analysis-Modélisation Mathématique et Analyse Numérique **47**(2), 377–399 (2013) https://doi.org/10.1051/m2an/2012031

[35] Bitsouni, V., Chaplain, M.A.J., Eftimie, R.: Mathematical modelling of cancer invasion: the multiple roles of TGF-β pathway on tumour proliferation and cell adhesion. Mathematical Models and Methods in Applied Sciences **27**(10), 1929–1962 (2017) https://doi.org/10.1142/S021820251750035X

[36] Bitsouni, V., Eftimie, R.: Non-local parabolic and hyperbolic models for cell polarisation in heterogeneous cancer cell populations. Bulletin of Mathematical Biology **80**(10), 2600–2632 (2018) https://doi.org/10.1007/s11538-018-0477-4

[37] Krasnianski, M., Painter, K., Surulescu, C., Zhigun, A.: Nonlocal and local models for taxis in cell migration: a rigorous limit procedure. Journal of Mathematical Biology **81**, 1251–1298 (2020) https://doi.org/10.1007/s00285-020-01536-4

[38] Delarue, F., Lagoutière, F., Vauchelet, N.: Convergence analysis of upwind type schemes for the aggregation equation with pointy potential. Annales Henri Lebesgue **3**, 217–260 (2020) https://doi.org/10.5802/ahl.30

[39] Bailo, R., Carrillo, J.A., Murakawa, H., Schmidtchen, M.: Convergence of a fully discrete and energy-dissipating finite-volume scheme for aggregation-diffusion equations. Mathematical Models and Methods in Applied Sciences **30**(13),
2487-2522 (2020) https://doi.org/10.1142/S0218202520500487

[40] Othmer, H.G., Dunbar, S.R., Alt, W.: Models of dispersal in biological systems. Journal of mathematical biology **26**(3), 263-298 (1988) https://doi.org/10.1007/BF00277392

[41] Chang, K., Zhang, Z., Liang, G.: Dynamics analysis of a nonlocal diffusion dengue model. Scientific Reports **13**, 15239 (2023) https://doi.org/10.1038/s41598-023-42440-3

[42] Chen, L., Painter, K., Surulescu, C., Zhigun, A.: Mathematical models for cell migration: a non-local perspective. Philosophical Transactions of the Royal Society B **375**(1807), 20190379 (2020) https://doi.org/10.1098/rstb.2019.0379

[43] Wang, H., Salmaniw, Y.: Open problems in pde models for knowledge-based animal movement via nonlocal perception and cognitive mapping. Journal of Mathematical Biology **86**(5), 71 (2023) https://doi.org/10.1007/s00285-023-01905-9

[44] Friedl, P., Locker, J., Sahai, E., Segall, J.E.: Classifying collective cancer cell invasion. Nature Cell Biology **14**(8), 777-783 (2012) https://doi.org/10.1038/ncb2548

[45] Ren, G., Roberts, A.I., Shi, Y.: Adhesion molecules: key players in mesenchymal stem cell-mediated immunosuppression. Cell Adhesion & Migration **5**(1), 20-22 (2011) https://doi.org/10.4161/cam.5.1.13491

[46] Burger, M., Capasso, V., Morale, D.: On an aggregation model with long and short range interactions. Nonlinear Analysis: Real World Applications **8**(3), 939-958 (2007) https://doi.org/10.1016/j.nonrwa.2006.04.002

[47] Trepat, X., Wasserman, M.R., Angelini, T.E., Millet, E., Weitz, D.A., Butler, J.P., Fredberg, J.J.: Physical forces during collective cell migration. Nature Physics **5**(6), 426-430 (2009) https://doi.org/10.1038/nphys1269

[48] Leclech, C., Barakat, A.I.: Is there a universal mechanism of cell alignment in response to substrate topography? Cytoskeleton **78**(6), 284-292 (2021) https://doi.org/10.1002/cm.21661

[49] Yamashita, Y.M., Inaba, M., Buszczak, M.: Specialized intercellular communications via cytonemes and nanotubes. Annual Review of Cell and Developmental Biology **34**, 59-84 (2018) https://doi.org/10.1146/annurev-cellbio-100617-062932

[50] Roehlecke, C., Schmidt, M.H.H.: Tunneling nanotubes and tumor microtubes in cancer. Cancers **12**(4) (2020) https://doi.org/10.3390/cancers12040857

[51] Fick, A.: Ueber diffusion. Annalen der physik (1855) https://doi.org/10.1002/
[52] Berdahl, A., Kao, A., Flack, A., Westley, P., Codling, E., Couzin, I., Dell, A., Biro, D.: Collective animal navigation and migratory culture: From theoretical models to empirical evidence. Philosophical Transactions of the Royal Society B: Biological Sciences **373**, 20170009 (2018) https://doi.org/10.1098/rstb.2017.0009

[53] Cavagna, A., Castello, L., Dey, S., Giardina, I., Melillo, S., Parisi, L., Viale, M.: Short-range interaction vs long-range correlation in bird flocks. Physical Review E **92** (2014) https://doi.org/10.1103/PhysRevE.92.012705

[54] Baltzersen, R.K.: Human Stigmergic Problem Solving. Cambridge University Press (2022). https://doi.org/10.1017/9781108981361.006

[55] Larkin, R., Szafoni, R.: Evidence for widely dispersed birds migrating together at night. Integrative and Comparative Biology **48**, 40–9 (2008) https://doi.org/10.1093/icb/icn038

[56] Cavagna, A., Cimarelli, A., Giardina, I., Parisi, G., Santagati, R., Stefanini, F., Viale, M.: Scale-free correlations in starling flocks. Proceedings of the National Academy of Sciences **107**(26), 11865–11870 (2010) https://doi.org/10.1073/pnas.1005766107

[57] Gorbonos, D., Ianconesu, R., Puckett, J.G., Ni, R., Ouellette, N.T., Gov, N.S.: Long range acoustic interactions in insect swarms: an adaptive gravity model. New Journal of Physics **18**, 073042 (2016) https://doi.org/10.1088/1367-2630/18/7/073042

[58] Bengtsson, J.: Aggregation in non-social insects. Växtskyddsbiologi, Sveriges lantbruksuniversitet. Alnarp, Sweden, 1–18 (2008)

[59] Mech, L.D., Boitani, L.: Wolves: Behavior, Ecology, and Conservation. University of Chicago Press (2019). https://doi.org/10.14430/arctic540

[60] Burt, W.H.: Territoriality and home range concepts as applied to mammals. Journal of Mammalogy **24**(3), 346–352 (1943) https://doi.org/10.2307/1374834

[61] Riotte-Lambert, L., Benhamou, S., Chamaillé-Jammes, S.: How memory-based movement leads to nonterritorial spatial segregation. The American Naturalist **185**(4), 103–116 (2015) https://doi.org/10.1086/680009

[62] Theraulaz, G., Bonabeau, E.: A brief history of stigmergy. Artificial Life **5**(2), 97–116 (1999) https://doi.org/10.1162/106454699568700

[63] Dyson, J., Gourley, S.A., Webb, G.F.: A non-local evolution equation model of cell–cell adhesion in higher dimensional space. Journal of Biological Dynamics
7(sup1), 68-87 (2013) https://doi.org/10.1080/17513758.2012.755572

[64] Gerisch, A., Painter, K.J.: Mathematical modeling of cell adhesion and its applications to developmental biology and cancer invasion, 337-368 (2010) https://doi.org/10.1201/9781420094558-c12

[65] Buttenschön, A., Hillen, T., Gerisch, A., Painter, K.J.: A space-jump derivation for non-local models of cell-cell adhesion and non-local chemotaxis. Journal of mathematical biology **76**, 429-456 (2018) https://doi.org/10.1101/093617

[66] Kress, R.: Numerical Analysis. Springer (2012). https://doi.org/10.1007/978-1-4612-0599-9

[67] Bitsouni, V., Trucu, D., Chaplain, M.A.J., Eftimie, R.: Aggregation and travelling wave dynamics in a two-population model of cancer cell growth and invasion. Mathematical Medicine and Biology: a Journal of the IMA **35**(4), 541-577 (2018) https://doi.org/10.1093/imammb/dqx019

[68] Steinberg, M.S.: Differential adhesion in morphogenesis: a modern view. Current Opinion in Genetics & Development **17**(4), 281-286 (2007) https://doi.org/10.1016/j.gde.2007.05.002. Pattern formation and developmental mechanisms

[69] Foty, R.A., Steinberg, M.S.: Cadherin-mediated cell-cell adhesion and tissue segregation in relation to malignancy. International Journal of Developmental Biology **48**(5-6), 397-409 (2004) https://doi.org/10.1387/ijdb.041810rf

[70] Logg, A., Mardal, K.-A., Wells, G.: Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book. Springer (2012). https://doi.org/10.1007/978-3-642-23099-8

[71] Bastian, P., Blatt, M., Dedner, A., Engwer, C., Klöfkorn, R., Ohlberger, M., Sander, O.: A generic grid interface for parallel and adaptive scientific computing. part i: abstract framework. Computing **82**, 103-119 (2008) https://doi.org/10.1007/s00607-008-0003-x

[72] Hecht, F.: New development in freefem++. Journal of Numerical Mathematics **20**(3-4), 251-265 (2012) https://doi.org/10.1515/jnum-2012-0013

[73] Cooper, F.R., Baker, R.E., Bernabeu, M.O., Bordas, R., Bowler, L., Bueno-Orovio, A., Byrne, H.M., Carapella, V., Cardone-Noott, L., Jonatha, C., et al.: Chaste: cancer, heart and soft tissue environment. Journal of Open Source Software **5**(47) (2020) https://doi.org/10.21105/joss.01848

[74] Petter, L.H., Logg, A.: Solving PDEs in Python. Springer (2017). https://doi.org/10.1007/978-3-319-52462-7
[75] Burman, E.: Robust error estimates in weak norms for advection dominated transport problems with rough data. Mathematical Models and Methods in Applied Sciences **24**(13), 2663-2684 (2014) https://doi.org/10.1142/S021820251450033X

[76] Allgower, E.L., Georg, K.: Introduction to Numerical Continuation Methods. SIAM (1990). https://doi.org/10.1137/1.9780898719154

[77] Rheinboldt, W.C.: Numerical continuation methods: a perspective. Journal of Computational and Applied Mathematics **13**, 103-113 (2000) https://doi.org/10.1016/S0377-0427(00)00428-3

[78] Kuehn, C.: Efficient gluing of numerical continuation and a multiple solution method for elliptic pdes. Applied Mathematics and Computation **266**, 656-674 (2015) https://doi.org/10.1016/j.amc.2015.05.120

[79] Uecker, H.: Continuation and bifurcation in nonlinear pdes – algorithms, applications, and experiments. Jahresbericht der Deutschen Mathematiker-Vereinigung **124**, 43–80 (2022) https://doi.org/10.1365/s13291-021-00241-5

[80] Le, T.T., Eftimie, R.: Numerical challenges for the understanding of localised solutions with different symmetries in non-local hyperbolic systems. Computers & Mathematics with Applications **169**, 112-131 (2024) https://doi.org/10.1016/j.camwa.2024.06.007

[81] Buono, P.-L., Eftimie, R., Kovacic, M., Veen, L.: Kinetic models for pattern formation in animal aggregations: A symmetry and bifurcation approach. Active Particles, Volume 2: Advances in Theory, Models, and Applications, 39–64 (2019) https://doi.org/10.1007/978-3-030-20297-2_2

[82] Braun, E.C., Bretti, G., Natalini, R.: Parameter estimation techniques for a chemotaxis model inspired by cancer-on-chip (coc) experiments. International Journal of Non-Linear Mechanics **140**, 103895 (2022) https://doi.org/10.1016/j.ijnonlinmec.2021.103895

[83] Falco, C., Cohen, D.J., Carrillo, J.A., Baker, R.E.: Quantifying tissue growth, shape and collision via continuum models and bayesian inference. Journal of the Royal Society Interface **20**, 20230184 (2023) https://doi.org/10.1098/rsif.2023.0184

[84] Campillo-Funollet, E., Madzvamuse, A., Venkataraman, C.: Bayesian parameter identification for turing systems on stationary and evolving domains. Bulletin of Mathematical Biology **81**, 81-104 (2019) https://doi.org/10.1007/s11538-018-0518-z

[85] Carrillo, J.A., Estrada-Rodriguez, G., Mikolas, L., Tang, S.: Sparse identification
of nonlocal interaction kernels in nonlinear gradient flow equations via par-
tial inversion. arXiv preprint arXiv:2402.06355 (2024) https://doi.org/10.48550/
arXiv.2402.06355

[86] Ellefsen, E., Rodríguez, N.: Nonlocal mechanistic models in ecology: Numerical methods and parameter inferences. Applied Sciences **13**(19) (2023) https://doi.org/10.3390/app131910598

[87] Alwuthaynani, M., Eftimie, R., Trucu, D.: Inverse problem approaches for mutation laws in heterogeneous tumours with local and nonlocal dynamics. Mathematical Biosciences and Engineering **19**(4), 3720–3747 (2022) https://doi.org/10.3934/mbe.2022171

[88] Greer, A.T., Chiaverano, L.M., Treible, L.M., Briseño-Avena, C., Hernandez, F.J.: From spatial pattern to ecological process through imaging zooplankton interactions. ICES Journal of Marine Science **78**(8), 2664–2674 (2021) https://doi.org/10.1093/icesjms/fsab149

[89] Lewis, S.M., Asselin-Labat, M.-L., Nguyen, Q., Berthelet, J., Tan, X., Wimmer, V.C., Merino, D., Rogers, K.L., Naik, S.H.: Spatial omics and multiplexed imaging to explore cancer biology. Nature Methods **18**, 997–1012 (2021) https://doi.org/10.1038/s41592-021-01203-6

[90] Magomedov, I., Khaliev, M.S., Elmurzaev, A.: Application of finite element analysis in medicine. In: Journal of Physics: Conference Series, vol. 1679, p. 022057 (2020). https://doi.org/10.1088/1742-6596/1679/2/022057. IOP Publishing

[91] Lee, C., Hoopes, M., Diehl, J., Gilliland, W., Huxel, G., Leaver, E., McCann, K., Umbanhowar, J., Mogilner, A.: Non-local concepts and models in biology. Journal of theoretical biology **210**(2), 201–219 (2001) https://doi.org/10.1006/jtbi.2000.228