# Numerical study on a multi-dimensional pressureless Euler-type model with non-local interactions and chemotaxis for collective cell migration

Marta Menci, Roberto Natalini, Tommaso Tenna

► To cite this version:

Marta Menci, Roberto Natalini, Tommaso Tenna. Numerical study on a multi-dimensional pressureless Euler-type model with non-local interactions and chemotaxis for collective cell migration. 2025. hal-05202105

HAL Id: hal-05202105

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

Preprint submitted on 6 Aug 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 study on a multi-dimensional pressureless
Euler-type model with non-local interactions and
chemotaxis for collective cell migration

M. Menci¹, R. Natalini², and T. Tenna³

¹Università Campus Bio-Medico di Roma, Rome, Italy

²Istituto per le Applicazioni del Calcolo, Consiglio Nazionale delle
Ricerche, Rome, Italy

³Université Côte d'Azur, Nice, France - Sapienza Università di Roma,
Rome, Italy

Abstract

In this paper we propose a numerical study of macroscopic models for collective cell migration, focusing on a multi-dimensional pressureless Euler-type model with non-local interactions coupled with chemotaxis, rigorously derived from microscopic dynamics. Different mechanical interactions are investigated, including attraction-repulsion effects. Moreover, the model is extended to the case of different populations of interacting cells. The validity of such macroscopic model and its agreement with the microscopic dynamics is finally assessed through a parameter estimation analysis in a specific setting.

# 1 Introduction

Many different real-world collective phenomena exhibit self-organization behaviors, usually driven by long-range mechanical interactions among agents. In the biological field, collective cell migration is also induced by chemical signals, detected by the cells and guiding their oriented motion. Several processes are regulated by the response of cell populations to the change of chemical concentrations in the environment and different mathematical models were developed for an accurate description of these phenomena [37, 45]. Starting from the seminal Keller-Segel model introduced in [28, 29], where the evolution of the density of cells is described by a parabolic equation and the concentration of a chemoattractant is given by a parabolic (or elliptic) equation, a wide number of variations have been proposed. Among others, in [39] authors present a minimal version of the original model, that has been analytically investigated in different directions. A more complete analysis and some variations can be found in [27, 44]. In more recent years, hyperbolic models have been developed. In particular, the models in [26, 47, 49] have been proposed to describe the early stages of the vasculogenesis process and they consist of a quasi-linear hyperbolic system, directly derived from macroscopic balance laws, coupled
with a parabolic equation for the chemoattractant. The advantage of the hyperbolic approach over parabolic Keller-Segel-type models relies in its ability to reproduce emerging pattern, i.e. the formation of the network of capillaries, as numerically observed in [43]. Moreover, the models account for the presence of a pressure term as a key factor, which authors consider to prevent overcrowding of cells, on a phenomenological basis.

Nevertheless, the behavior of cells at the microscopic scale seems to play a crucial role in modeling collective cell migration, providing a detailed description of the mechanical interactions among cells. In particular, second-order differential models allow to take into account the crucial role of persistence effects. This key feature of cells behavior is not observed in other collective phenomena, for instance in pedestrian dynamics, where accelerations are almost instantaneous and inertia-based effects are negligible. In recent years, a novel class of *hybrid coupled models* has been proposed, gathering the advantages of a multiscale description [20–22]. These models are based on a *discrete in continuous* approach: cells are described as discrete entities moving in the environment and chemotactic signals, influencing the dynamics, as continuous concentrations. Microscopic and macroscopic scales are strongly coupled: the dynamics of cells are influenced by the gradient of the chemical concentration and, at the same time, the evolution of the chemical signal is influenced by cells positions. From the analytical point of view, well-posedness results for generalized versions of hybrid models have been investigated over the last years [32–35].

From the modeling and numerical perspective, the detailed level of description of the microscopic approach can become unsuitable if the phenomena under investigation involve high number of cells. A macroscopic approach, modeling the evolution in time of macroscopic quantities of interest, as cell density and velocity field, better fulfills the requirements. For coupled hybrid systems, the issue concerning the macroscopic counterpart of the microscopic dynamics naturally arises. In this perspective, the recent results in [42] represent the first step toward rigorous derivation of a macroscopic model. Here authors rigorously proved the mean-field limit of the microscopic model towards a kinetic equation coupled with chemotaxis, and its hydrodynamic limit, under strong assumptions. The main feature of the resulting Euler-type model is the absence of pressure terms, in the face of non-local integral terms which keep into account the microscopic interactions, establishing a clear link between the macroscopic model and the microscopic dynamics. As a drawback, the absence of the pressure term in Euler-type systems leads to possible blow-up of the solutions [4, 7, 11]. Non-local pressureless Euler systems have been derived and investigated in different directions [1, 25, 36, 41]. In particular, a wide number of results deals with alignment kind of interactions, inspired by the seminal microscopic alignment model in [18].

**Aim of the paper** In this paper we focus on the macroscopic scale of description, in a multi-dimensional setting. A numerical investigation of different models is performed, comparing seminal models of the literature with the recent pressureless Euler-type model coupled with chemotaxis derived in [42]. The absence of pressure enables the development of a model that does not require phenomenological terms, relying instead on nonlocal terms that retain information about microscopic interactions. In this perspective, we extend the kind of mechanical interactions included in the model, considering attraction-repulsion effects together with alignment. The non-locality of interactions and the coupling with the chemical signal represent the main features of the model with respect to previous works of the literature. Moreover, the original structure of the dynamics
is extended to a multi-population case, accounting for different types of interacting cells. This choice allows the construction of more complex and realistic models, where non-local interactions are characteristic of each individual species. Inspired by successful applications of microscopic hybrid models in simulating Cancer-On-Chip dynamics [8,9], we test the applicability of the macroscopic model in this context through a calibration procedure. Due to the current lack of experimental data on the macroscopic scale, we present an estimation procedure linking the microscopic dynamics, already calibrated on real experimental data, and the macroscopic one. This strategy enables to generate synthetic data from the microscopic model, in order to obtain a greater number of different scenarios to assess the macroscopic one.

**Structure of the paper** The rest of the paper is organized as follows. In Section 2 we analyze in details the coupled models, starting from the microscopic scale to the macroscopic one. In particular, we present the extension of the non-local Euler model to the case of two populations of cells and including different interactions. The numerical study of the models is deferred to Section 3, aiming at investigating and comparing the role of the different terms. Section 4 is devoted to parameter estimation for the model introduced in Section 2.3. All the details about the numerical schemes adopted for the approximation of the solutions to the different models are reported in Appendix A.

# 2 Description of the mathematical models

In this section we delve in details the mathematical models for chemotactic cell migration considered in our study, highlighting the role of the different terms. In the following, we distinguish between microscopic and macroscopic coupled models, referring to the scale of observation considered for cell dynamics. In all the considered models, the evolution of chemoattractant is modeled by a linear diffusion equation with a source term which depends on the positions of cells (for the microscopic one) or the density of cells (for the macroscopic ones). We start by reviewing existing models from the literature, then incorporate additional nonlocal effects, both for the microscopic and the macroscopic scale. Finally, we extend the model to account for multiple interacting populations, focusing on the dynamics between immune and cancer cells at the macroscopic level.

## 2.1 Microscopic models

In all the works of the literature [20–22], cells are modeled as circles, endowed with a radius $R > 0$, whose dynamics is described by time update of positions and velocities of their centers. From a mathematical point of view, the general structure of a hybrid model consists of a system of second order ordinary differential equations, modeling the evolution in time of position $\mathbf{x}_i \in \mathbb{R}^d$ and velocity $\mathbf{v}_i \in \mathbb{R}^d$ of the $i$-th cell, $i = 1, ..., N$, coupled with a reaction-diffusion parabolic partial differential equation which describes the evolution of the concentration of the chemoattractant $\varphi = \varphi(x, t)$. The model reads as:

$$ egin{cases} 
\dot{\mathbf{x}}_i = \mathbf{v}_i, \\
\dot{\mathbf{v}}_i = \mathbf{F}_i(t, \mathbf{X}(t), \mathbf{V}(t), \nabla\varphi(t, \mathbf{x}_i)), \\
\partial_s\varphi = D\Delta\varphi - \kappa\varphi + g(x, \mathbf{X}(s)), \quad s \in [0, t],
\end{cases} 
\quad (1) $$
where $\mathbf{X} = (\mathbf{x}_1, \dots, \mathbf{x}_N) \in \mathbb{R}^{Nd}$, $\mathbf{V} = (\mathbf{v}_1, \dots, \mathbf{v}_N) \in \mathbb{R}^{Nd}$ denote the collections of position and velocity of each cell at any time in the time interval considered, $t \in [0, T]$. The chemotaxis equation accounts for constant diffusion and degradation coefficients, respectively denoted with $D, \kappa > 0$. The source term $g$ models the fact that the chemical signal is produced by the cells themselves through the dependence on the cells positions. In [20, 22, 31, 33] $g$ is defined as the sum of characteristic functions over circular domains of fixed radius $R > 0$ modeling cells. Throughout the rest of the paper, we keep the space variables in math mode style for simplicity of notation. Functions $F_i$ model mechanical and chemotactic interactions influencing the dynamics of each cell $i = 1, \dots, N$. In particular, in [21, 42] $F_i$ are of the following form:

$$ \mathbf{F}_i(t, \mathbf{X}(t), \mathbf{V}(t), \nabla\varphi(t, \mathbf{x}_i)) = \frac{1}{N} \sum_{j=1}^{N} \gamma(\mathbf{v}_i - \mathbf{v}_j, \mathbf{x}_i - \mathbf{x}_j) + \eta \nabla\varphi(t, \mathbf{x}_i) + F_{\text{ext}}(\mathbf{x}_i, \mathbf{v}_i), \quad (2) $$

where $\gamma: \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}^d$ is the interaction function, $\varphi: [0, T] \times \mathbb{R}^d \to \mathbb{R}^d$ is the concentration of the chemoattractant and $F_{\text{ext}}$ models external forces. Initial data are given by initial position and velocity of each particle, and initial concentration of chemoattractant

$$ \mathbf{X}(0) := \mathbf{X}_0, \quad \mathbf{V}(0) := \mathbf{V}_0, \quad \varphi(x, 0) := \varphi_0. \qquad (3) $$

Starting from the hybrid coupled system in (1), different kinds of cell-cell interactions have been considered, depending on the specific applications. In the literature, the great majority of agent-based models accounts for alignment and/or attraction-repulsion kind of interactions. In order to describe different contributions in (2), we rewrite the interaction function $\gamma$ as the sum of two contributions, $\gamma = \gamma_1 + \gamma_2$. Among alignment models, the seminal Cucker-Smale model has been conceived to model consensus in flocks of birds and animal groups [17–19]. The force acting on each individual is a weighted average of the differences of its velocity with the ones of each other member of the group, and it is proportional to the inverse of the distance between agents. The interaction function is hence given in this form:

$$ \gamma_1(\mathbf{v}_i - \mathbf{v}_j, \mathbf{x}_i - \mathbf{x}_j) = \gamma_D(||\mathbf{x}_i - \mathbf{x}_j||) (\mathbf{v}_j - \mathbf{v}_i), \quad (4) $$

where $\gamma_D$ is a decreasing function of the distance $||\mathbf{x}_i - \mathbf{x}_j||$ between two agents. In the following we assume $\gamma_D$ as the classical Cucker-Smale communication rate

$$ \gamma_D(||\mathbf{x}_i - \mathbf{x}_j||) = \frac{\beta}{(1 + ||\mathbf{x}_j - \mathbf{x}_i||^2)^\zeta}, \quad (5) $$

where $\beta \in \mathbb{R}^+$ and $\zeta \in \mathbb{R}^+$ is the Cucker-Smale exponent. Concerning attraction and repulsion effects, interactions are often modeled by interaction potentials [12, 23] or introducing specific functions of the distance, modeling short-range repulsion and long-range attraction. In particular, in [20, 22, 48] authors assume an attractive linear elastic force, and a repulsion force proportional to the inverse of the distance between any two agents, defining a piecewise linear function:

$$ \gamma_2(\mathbf{x}_i - \mathbf{x}_j) = \begin{cases} -\omega_{\text{rep}} \left( \frac{1}{||\mathbf{x}_i - \mathbf{x}_j||} - \frac{1}{R_{\text{rep}}} \right) \frac{\mathbf{x}_j - \mathbf{x}_i}{||\mathbf{x}_i - \mathbf{x}_j||}, & \text{if } ||\mathbf{x}_i - \mathbf{x}_j|| \le R_{\text{rep}}, \\ \omega_{\text{adh}} (||\mathbf{x}_i - \mathbf{x}_j|| - R_{\text{rep}}) \frac{\mathbf{x}_j - \mathbf{x}_i}{||\mathbf{x}_i - \mathbf{x}_j||}, & \text{if } R_{\text{rep}} < ||\mathbf{x}_i - \mathbf{x}_j|| \le R_{\text{adh}}, \end{cases} \quad (6) $$
with $\gamma_2 \equiv 0$ if $||\mathbf{x}_i - \mathbf{x}_j|| > R_{adh}$. Here $\omega_{rep}$, $\omega_{adh}$, $R_{rep}$, $R_{adh}$ are constant positive parameters. In particular repulsion occurs when the distance between two agents is less than $R_{rep} > 0$, whereas adhesion occurs at a distance greater than $R_{rep}$ and less than $R_{adh} > R_{rep}$. Several extensions have been proposed, generalizing (1). In [22], a hybrid system with non-local concentration is introduced, accounting for the fact that the gradient of chemical concentration influences cells not only in their centers, but also in their surroundings. The non-locality of the influence of chemical concentration leads to a more complex structure, which have been rigorously analyzed from a theoretical view point in [34].

## 2.2 Macroscopic models

The quasi-linear hyperbolic model for cell migration proposed in [26] was originally conceived to describe vasculogenesis processes. The structure of the model is the following:

$$ \left\{ \begin{array}{l} \partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0, \\ \partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) = \eta \rho \nabla \varphi - \nabla P(\rho) - \alpha \rho \mathbf{u}, \\ \partial_t \varphi = D \Delta \varphi - \kappa \varphi + g(\rho), \end{array} \right. \qquad (7) $$

where $\eta$, $\alpha$, $D$, $\kappa$ are positive constants. The term $\alpha \rho \mathbf{u}$ represents a damping term, which models dissipative interactions with the external environment. In the non-local macroscopic model (11), this force could be included in $F_{ext}$. The function $g$ models the production rate of the chemoattractant. In [26] a constant production is assumed, i.e. $g(\rho) = a \rho$ with $\rho > 0$. In [49] authors choose $g(\rho) = \alpha_1 \rho/(1 + \alpha_2 \rho^2)$, with $\alpha_1$, $\alpha_2 > 0$, modeling the fact that, in aggregate state where cells are packed and density is high, the chemoattractant production almost vanishes. The first equation is the classical mass conservation, since cells do not undergo mitosis or apoptosis during the considered phenomenon, hence the total number of cells is constant in time. The right-hand side of the second equation, modeling the different effects that influence changes in cell directional motion, is the sum of the following different contributions:

1. a chemotactic body force

$$ f_{chem} = \eta \rho \nabla \varphi, \qquad (8) $$

where $\eta > 0$ measures the strength of the cells response to concentration gradients of the chemical substance. The linear dependence on $\rho$ corresponds to the assumption that each cell experiences a similar chemotactic action. In [46] a saturation effect on the amount of chemoattractant is also included.

2. a dissipative interaction with the substrate

$$ f_{diss} = -\alpha \rho \mathbf{u}. \qquad (9) $$

A linear dependence of $f_{diss}$ on $\rho$ is assumed, and its strength is tuned through the constant parameter $\alpha > 0$.

3. the main feature relies on the pressure term

$$ f_{press} = -\nabla P(\rho), \qquad (10) $$

which has been introduced with the aim of preventing overcrowding of cells. In this sense $f_{press}$ is a phenomenological term modeling no-collision dynamics of cells. The
main issue is indeed the choice of the function $P$, which is not justified by empirical evidence. In [13, 26, 49] the structure is inspired to the pressure law for isentropic gases, i.e. $P(\rho) \propto \rho^{\Gamma}$, $\Gamma > 1$.

Let us now focus on the macroscopic model derived in [42], which describes the evolution of cell density $\rho = \rho(x, t)$ and velocity field $\mathbf{u} = \mathbf{u}(x, t)$ through the following non-local pressureless Euler-type system coupled with chemotaxis

$$ 
\left\{ 
\begin{array}{l}
\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0, \\
\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) = \rho \mathcal{I} + \eta \rho \nabla \varphi + \rho F_{ext}, \\
\partial_t \varphi = D \Delta \varphi - \kappa \varphi + g(\rho).
\end{array}
\right. 
\qquad (11) 
$$

Here $\eta$, $D$, $\kappa$ are positive constants and $\mathcal{I} = \mathcal{I}(x, t)$ denotes a non-local integral term, which in [42] is defined as

$$ \mathcal{I}(x) = \int_{\mathbb{R}^d} \gamma(\cdot - y, \mathbf{u}(\cdot) - \mathbf{u}(y))\rho(y) dy. \qquad (12) $$

The main difference with respect to (7) is the absence of the pressure term, formally replaced by the non-local integral term $\mathcal{I}$ in (12). From a mathematical point of view, it is well-known that the presence of pressure may lead to some a priori estimates, preventing blow up of solutions and resulting as a stabilizing force [30]. The key factor of (11) relies on the non-local integral term (12), which highlights the mechanical kind of interactions acting in the microscopic dynamics through the function $\gamma$. Concerning an Euler-alignment system with the interaction function inspired by the Cucker-Smale model, the integral term can be rewritten as

$$ \mathcal{I}_1(x) = \int_{\mathbb{R}^d} \frac{\beta}{(1 + ||x - y||^2)^\zeta} (\mathbf{u}(y) - \mathbf{u}(x)) \rho(y) dy. \qquad (13) $$

Changing the integral function, one can simulate different kind of interactions. In particular, we will consider repulsion-attraction among cells, using the interaction function $\gamma_2$ in (6). In this case, the integral term can be written as

$$ \mathcal{I}_2(x) = \int_{\mathbb{R}^d} \gamma_2(x - y) \rho(y) dy. \qquad (14) $$

In Section 3 we present numerical simulations of the two macroscopic models here considered, in order to compare their specific features. The rigorous derivation of (11) from the microscopic dynamics (1) can be found in [42].

## 2.3 Extension to two-populations models

In a microscopic description of the dynamics, the hybrid model (1) has been extended to the case of two populations of interacting cells. In [9], authors presented a hybrid model describing immune cells migration under the effects of chemicals released by cancer cells, which are fixed in the domain. The cell-cell interactions for the immune system dynamics is described by (4)-(6). The dynamics of immune cells towards cancer cells is described
by a radial repulsion term, depending on cells position. Let $\mathbf{y}_j$ be the position of the j-th cancer cell, then function $\gamma_3$ is introduced as

$$ \gamma_3(\mathbf{x}_i - \mathbf{y}_j) := -\omega_{\text{rep}}^{\text{tum}} \left( \frac{1}{\|\mathbf{x}_i - \mathbf{y}_j\|} - \frac{1}{R_{\text{rep}}^{\text{tum}}} \right) \frac{\mathbf{y}_j - \mathbf{x}_i}{\|\mathbf{x}_i - \mathbf{y}_j\|}, \quad \text{if } \|\mathbf{x}_i - \mathbf{y}_j\| \le R_{\text{rep}}^{\text{tum}}. \quad (15) $$

and $\gamma_3 \equiv 0$ if $\|\mathbf{x}_i - \mathbf{y}_j\| > R_{\text{rep}}^{\text{tum}}$, where $\omega_{\text{rep}}^{\text{tum}}$. Here $R_{\text{rep}}^{\text{tum}}$ are constant positive parameters. With this choice, repulsion occurs when the distance between the i-th immune cell and the j-th cancer cell is less than $R_{\text{rep}}^{\text{tum}}$. The final model for the description of immune cells migration in presence of cancer cells can be written as [9]:

$$ \left\{ \begin{array}{l} \dot{\mathbf{x}}_i = \mathbf{v}_i, \\ \\ \dot{\mathbf{v}}_i = \frac{1}{N} \sum_{j=1}^{N} (\gamma_1(\mathbf{v}_i - \mathbf{v}_j, \mathbf{x}_i - \mathbf{x}_j) + \gamma_2(\mathbf{x}_i - \mathbf{x}_j)) \\ \\ \qquad + \sum_{j=1}^{M} \gamma_3(\mathbf{x}_i - \mathbf{y}_j) + \eta \nabla \varphi(t, \mathbf{x}_i) + F_{\text{ext}}(\mathbf{x}_i), \\ \\ \partial_s \varphi = D \Delta \varphi - \kappa \varphi + \sum_{j=1}^{M} \chi(x - \mathbf{y}_j), \quad s \in [0, t], \end{array} \right. \quad (16) $$

where $x \in \mathbb{R}^d$, $N$ is the total number of immune cells, $M$ is the total number of cancer cells of radius $R_{\text{tum}}$ and $\chi \in C_c^1$ is a generic function, usually chosen as a weighted regularization of the characteristic function

$$ \mathbb{1}(z) = \begin{cases} 1 & \text{if } \|z\| < R_{\text{tum}}, \\ 0 & \text{otherwise.} \end{cases} $$

System (16) is complemented by initial conditions, given by microscopic positions and velocities of each immune cell:

$$ \mathbf{x}_i(0) = \mathbf{x}_i^0, \qquad \mathbf{v}_i(0) = \mathbf{v}_i^0. \quad (17) $$

Since cancer cells do not move in the considered setting, the initial condition for the parabolic equation is given by the chemoattractant $\varphi$ produced at time $t = 0$, namely

$$ \varphi(x, 0) = \sum_{j=1}^{M} \chi(x - \mathbf{y}_j). \quad (18) $$

Inspired by the structure of (11) and its link with the microscopic scale previously described, we here propose a macroscopic Euler-type model keeping memory of the key factor of (16). More specifically, the source term of the parabolic equation for the chemoattractant here depends on the density of cancer cells, as in (16). We denote by $\rho(x, t)$ and $\zeta(x, t)$ respectively the density of immune and cancer cells, respectively, at time $t$ in $x$. To model interactions between immune and cancer cells in (20), we introduce a non-local integral term depending on the microscopic interaction modeled by $\gamma_3$ in (11), namely

$$ \mathcal{I}_3(x) = \int_{\mathbb{R}^d} \gamma_3(x - y)\zeta(y) dy. \quad (19) $$
This term links the microscopic interactions among cells with the macroscopic density evolution. The macroscopic two populations Euler-type model reads as

$$ egin{cases} 
\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0, \\
\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) = \rho \mathcal{I} + \eta \rho \nabla \varphi + \rho F_{ext}, \\
\partial_t \varphi = D \Delta \varphi - \kappa \varphi + \chi * \zeta,
\end{cases} 
\quad (20) $$

where the non-local integral term $\mathcal{I}$ is the sum of the alignment and attraction-repulsion contributions,

$$ \mathcal{I} = \mathcal{I}_1 + \mathcal{I}_2 + \mathcal{I}_3, \quad (21) $$

with $\mathcal{I}_1$, $\mathcal{I}_2$ and $\mathcal{I}_3$ respectively defined in (13), (14) and (19) and cancer cells density $\zeta$ is assumed to be constant in time.

## 3 Macroscopic coupled models: 2D numerical simulations

In this section we present 2D numerical solutions to the models introduced in the previous sections. The choice of the two-dimensional domain lies in the fact that most of the experimental data at our disposal refers to bi-dimensional settings. Indeed, biological experiments on two-dimensional devices offer a cost-effective tool to study cells interactions in a microfluidic environment. To the best of our knowledge, up to now only 1D numerical simulations of the model in (11) can be found in the literature [31], and accounting only for alignment kind of interactions.

We approximate the solutions through finite difference schemes on a square spatial domain $\Omega = [0, 1] \times [0, 1]$, with a mesh size $\Delta x$ in x-direction and $\Delta y$ in y-direction. In all simulations we fix $\Delta x = \Delta y = 0.02$. The time interval is set to $[0, T]$, with $T = 1$ and the numerical scheme adopt a time-variable discretization step $\Delta t$, according to the CFL condition. The numerical scheme is based on discrete kinetic approximations of the models, in which the transport term and the relaxation term are treated using a splitting technique [3, 40]. More precisely, the finite difference scheme for the transport term involves flux limiters in space, improving the numerical scheme conceived in [31] for the 1D case. The numerical approximation of the parabolic equation for the chemotaxis is based on a standard centered scheme. A detailed description of the numerical scheme is given in Appendix A.

The initial macroscopic density is obtained as the sum of N Gaussian bumps whose amplitude is assumed to be of the order of a positive parameter $\sigma$. Heuristically, in the context of cell migration phenomena, $\sigma$ represents the non-dimensional average cell radius: each Gaussian bump, centered in random positions $\mathbf{x}_j$, $j = 1, \dots, N$, approximates a cell endowed with its radius $\sigma$. Thus, the initial density is given by

$$ \rho(x, 0) = \frac{\tau}{2\pi\sigma^2} \sum_{j=1}^{N} \exp\left(-\frac{|x - \mathbf{x}_j|}{2\sigma^2}\right), \quad (22) $$

where $\tau$ is a positive normalization parameter. We set $N = 200$ and $\sigma = 0.015$ to obtain the initial density shown in Fig. 1. In the following tests we assume the same initial condition, setting null initial conditions for the velocity and the chemoattractant

$$ \rho_0(x) = \rho(x, 0) \quad \mathbf{u}(x, 0) = \mathbf{u}_0(x) = 0, \quad \varphi(x, 0) = \varphi_0(x) = 0 \quad \forall x \in \Omega. \quad (23) $$
We assume homogeneous Neumann boundary conditions for the density $\rho$ and for the chemoattractant $\varphi$, and zero boundary conditions for the normal component of the velocity field $\mathbf{u}$:

$$ \nabla \rho \cdot \mathbf{n} \Big|_{\partial\Omega} = 0, \quad \nabla \varphi \cdot \mathbf{n} \Big|_{\partial\Omega} = 0, \quad \mathbf{u} \cdot \mathbf{n} \Big|_{\partial\Omega} = 0. \qquad (24) $$

We present different representative numerical simulations of the models in (11) and (7), showing the role of different effects involved in the dynamics.

In all of the following numerical simulations we fix $\alpha = 0$, $\kappa = 1$, $D = 0.1$, $\alpha_1 = 30$, $\alpha_2 = 0.2$, unless otherwise mentioned.

## 3.1 Pressureless Euler with chemotaxis without non-local interactions

We start our analysis simulating (11) neglecting the non-local interaction term $\mathcal{I}$ together with external forces. The structure of the system hence reduces to a pressureless Euler system, coupled with chemotaxis. In our simulations, the finite-time blow-up which characterizes pressureless dynamics is observed, coherently with previous works of the literature on pressureless Euler (without chemotactic interactions) [1, 5, 10].

Our focus is on the coupling with chemotaxis, hence we perform numerical simulations, varying the value of the parameters $\eta$. We observe that the chemotactic effect acts as an aggregating factor, directing the velocity field towards points of higher density. Fig. 2 refers to the results for three different values of $\eta$. In particular, the snapshots of the density are acquired at the blow-up time instant. We can observe that the time required to detect a blow-up in the solution decreases with respect to $\eta$. This is reasonable, since higher values of the investigated parameter correspond to a strongest influence of the chemotactic gradient on the density dynamics.

## 3.2 Isentropic Euler with chemotaxis

Let now focus on the model in (7). In relevant works of the literature of vasculogenesis [2, 26, 49], the presence of the phenomenological pressure term together with the evolution of the chemical signal drives the formation of a network structure. In particular, the pressure term acts as a repulsion force, meaning that overcrowding regions are prevented. The structure of the function $P(\rho)$ to be considered is widely debated. In this simulation, we assume the pressure function introduced in [13], i.e.

$$ P(\rho) = \begin{cases} (\rho - \rho_0)^3 & \text{if } \rho > \rho_0, \\ 0 & \text{if } \rho \le \rho_0, \end{cases} \qquad (25) $$

with $\rho_0 = 4$. The presence of an activation threshold $\rho_0$ models the fact that the pressure acts as a cell repulsion at sufficiently high cell density, whereas no effect is expected at lower densities. In Fig. 3, we observe the evolution in time of the density $\rho$, starting from the initial condition in Fig. 1. The network formation is strictly dependent to the threshold $\rho_0$ introduced in the expression for the pressure (25): cell migration toward other cells is strongly damped when the density reaches its threshold, since the pressure term to prevent overcrowding of cells becomes non-zero.
### 3.3 Pressureless Euler with chemotaxis and non-local interactions

We numerically investigate the non-local pressureless Euler-type model in (11), with the non-local interaction functions given in (4)-(6). We consider both alignment and attraction-repulsion effects as mechanical interactions. The structure of functions $\gamma_1$ and $\gamma_2$ are modeling choices, but depending on the ones chosen (or experimentally observed) at the microscopic scale. This is the main difference and a key factor of the model with respect to (11), where the choice of the pressure function is pure phenomenological.

Fig. 4 shows different snapshots of a numerical simulation of (11). From a qualitative perspective, numerical evidence shows that the pressureless model is still able to reproduce network formation. Formally, from a mathematical point of view, the pressure term is indeed replaced by the non-local interactions. With respect to the numerical simulation in Fig. 3, it seems that the presence of a mechanical attraction and alignment effect, in addition to chemotaxis, leads to a less dissipative behavior of the network.

<figure><img src="image_2.png" /><figcaption>Figure 1: Initial condition for $\rho$ in all simulations.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Numerical simulations of (11) with $\mathcal{I} = 0$, for two different values of $\eta$. Increased values of $\eta$ result in higher maximum density values, reducing the time at which the blow-up is observed.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Numerical simulation of (7) with $\eta = 0.5$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Numerical simulation of (11) with $\eta = 0.5$.</figcaption></figure>

# 4 Parameter Estimation

In this section we present the procedure followed to estimate the parameters of the model (20), focusing on the key parameters that are not given in the literature. In particular, we validate the model on a reference scenario inspired by Cancer-on-Chip experiments, hence involving two populations of interacting cells. Using parameter values available at the microscopic scale, obtained in [8,9], we present an optimization procedure to estimate the unknown parameters, after converting the information at the microscopic scale to the macroscopic one. In order to reduce the number of parameters to be estimated, we fix those ones known from the literature or already experimentally determined, such as cells radii, see in [9]. In our study, synthetic initial data are used to perform the calibration. We remark that, whereas the hybrid model has been already validated against real data at the microscopic scale, this is the first study concerning the calibration of the macroscopic model.

In our approach we formulate a minimization problem, aiming at finding the optimal values for the parameters of interest, minimizing the “distance” between synthetic (or experimental) data and the numerical approximation of the density solution to (20). In fact, the density of cells represents the quantity of interest for clinicians, as well as the information which can be easily acquired through the experiments.
## 4.1 Linking microscopic data to macroscopic solution

The first step of the calibration procedure consists of generating synthetic data through the hybrid model (16), using the literature parameters available. Synthetic data will play the role of experimental data in the formulation of the parameter estimation problem. The kind of data we are interested in are the trajectories of cells, i.e. positions of cells at each time of the simulation. The crucial point relies in the comparison between the microscopic information and the macroscopic counterpart. Following the idea in [6], we adopt a Kernel Density Estimation (KDE) approach, which is known to provide an efficient tool estimate a probability density function from a sample of observation. In our approach, we will make use of the KDE to create a macroscopic (density) counterpart starting from a discrete set of microscopic positions of cells.

**Definition 4.1** Let $\mathcal{X} = \{\mathbf{x}_1, \dots, \mathbf{x}_m\}$ be a sample set of n-dimensional random vectors $\mathbf{x}_i$, distributed according to an unknown probability density function $f$. Then the multivariate kernel density estimate $\hat{f}_{H,n}$ of $f$ is defined as

$$ \hat{f}_{H,n}(x) = \frac{1}{m} \sum_{i=1}^{m} K_{H,n}(x - \mathbf{x}_i) = \frac{1}{m \det H} \sum_{i=1}^{m} K(H^{-1}(x - \mathbf{x}_i)), \quad (26) $$

where $x, \mathbf{x}_i \in \mathbb{R}^n$, the bandwidth matrix $H \in \mathbb{R}^{n \times n}$ and multivariate kernel function $K$.

The bandwidth matrix $H$ is symmetric and positive definite, and the kernel function $K_{H,n}$ is a normalized non-negative and bounded function for all $x \in \mathbb{R}^n$.

In our study $\mathcal{X}$ is the set of positions of immune cells, dealing with the two-dimensional case. For the kernel function, we choose the following multivariate Gaussian kernel function

$$ K(x) = (2\pi)^{-\frac{n}{2}} e^{-\frac{1}{2}x^T x}. \quad (27) $$

Rigorous results on the convergence of the kernel density estimate $\hat{f}_{H,n}(x)$ towards a probability density function $f$ with a bounded variation kernel function in the one-dimensional setting can be found in [14, 38].

We point out that choice of the coefficients of the bandwidth matrix $H$ in the multi-dimensional case is fundamental to control the smoothness of the estimate. In the one-dimensional case, where the bandwidth matrix $H$ reduces to one parameter $h$, a large value of $h$ yields to an oversmoothed estimate, increasing the bias and losing some important features in the estimation. On the other hand, a small $h$ induces an undersmoothed estimate, increasing the variance and including misleading features introduced by noise. In our 2D setting, we assume independent smoothing along the $x$ and $y$ directions, obtaining $H$ in a diagonal form. In details, we assume

$$ H = \mathrm{diag}(h\sigma, h\sigma), \quad (28) $$

where $\sigma$ is the radius of the cells and $h \in \mathbb{R}$ is one of the parameter that will be estimated in the procedure.

In Fig. 5 we show the macroscopic density obtained from a set of random microscopic positions of $m = 100$ cells through the KDE approach, for two different choices of $h$.
<figure><img src="image_6.png" /><figcaption>Figure 5: KDE starting from a set of random microscopic positions (red dots), for two different values of h in (28).</figcaption></figure>

## 4.2 Minimization Problem

To perform the parameter estimation, we define the objective function to be minimized as

$$J(\Theta) = \frac{1}{N_t} \sum_{i=1}^{N_t} \left( \frac{\|\rho(t_i, \Theta) - \rho_{\text{data}}(t_i)\|_2}{\|\rho_{\text{data}}(t_i)\|_2} \right)^2, \quad (29)$$

where $N_t$ is the number of time steps involved in the minimization process, $\rho_{\text{data}}$ is the reference solution obtained through the KDE described in the previous section and $\Theta \in \mathbb{R}^p$ is the vector of $p$ parameters to be determined. In order to overcome possible ill-posedness, we regularize the problem by including a Tikhonov regularization term as [6, 15]:

$$R_{\lambda}(\Theta) = \lambda^2 \|\Theta - \Theta_0\|_2^2, \quad (30)$$

where $\Theta_0 \in \mathbb{R}^p$ is a priori estimate of the parameters involved in the estimation procedure. Indeed, the presence of the regularization term allows to damp contributions from data errors and rounding errors, keeping the norm of $\lambda \|(\Theta - \Theta_0)\|$ of reasonable size [24]. In the following, we will choose small values of $\lambda$ such that the norm of $\|\Theta - \Theta_0\|$ will not give a significant contribution in the parameter estimation.

The estimation can be seen as an inverse problem, aiming at solving

$$\Theta_{opt} = \arg \min_{\Theta \in \mathcal{A}} K(\Theta), \quad (31)$$

where $\mathcal{A}$ is the set of admissible solutions and $K$ the following nonlinear regularized objective function

$$K(\Theta) = J(\Theta) + R_{\lambda}(\Theta). \quad (32)$$

The minimization is performed using a trust region algorithm, which provides as solution the optimal value $\Theta_{opt}$. We assess the goodness of the parameter estimation by computing $E = J(\Theta_{opt})$, namely the relative error between synthetic data and the solution obtained with the estimated parameters.

## 4.3 Numerical Results

In the following we present the simulation results obtained through the parameters estimation procedure for the macroscopic model (20), accounting for two population of cells.
In order to generate synthetic data starting from the microscopic level, we run (16) in the time interval [0, 1] considering $N = 80$ immune cells, whose initial positions are chosen randomly in the domain $[0, 1] \times [0, 1]$ and $M = 3$ tumor cells, located in $\mathbf{y}_1 = [0.5, 0.7]$, $\mathbf{y}_2 = [0.3, 0.3]$, $\mathbf{y}_3 = [0.8, 0.5]$, respectively. Initial null velocities of cells is assumed and initial chemical concentration is chosen as in (18). The value of the parameters considered for the reference test are shown in Table 1, corresponding to adimensional values of the parameters in [9]. The numerical approximation of the solution to the hybrid model, used as synthetic data in this Section, is based on an implicit Euler scheme for the ODE in (16) and a Crank-Nicolson finite difference scheme for the parabolic equation, with spatial discretization performed with $\Delta x = \Delta y = 0.02$. We refer to [8] for a detailed description of the numerical schemes for the hybrid model. Starting from initial positions of immune and cancer cells, namely $\mathbf{x}_i$ and $\mathbf{y}_j$, the KDE is used to build the initial density $\rho_0$ to run the macroscopic model. We recall that cancer cells are fixed in the domain, with $\mathbf{y}_j(t) = \mathbf{y}_j$ for any $t$, hence the KDE for the cancer cell is done offline before starting the optimization procedure. In the following tests we set the bandwidth parameter value $h$ in (28) for the three cancer cells equal to 1.2, while estimating the parameter value for the immune cells. Fig. 6 shows the initial density $\rho_0$ plotted in the $xy$-plane, together with the positions of the immune cells, represented as circles endowed with a radius $R_{imm}$.

<figure><img src="image_7.png" /><figcaption>Figure 6: Plot of initial density $\rho_0$ and initial position of immune cells (black circles).</figcaption></figure>

The parameters to be estimated in (20) are

$$ \Theta = [\eta \quad \omega_{rep} \quad \omega_{adh} \quad \beta \quad \omega_{rep}^{tum} \quad h] \in \mathbb{R}_{+}^{6}. \quad (33) $$

We recall that $\eta$ tunes the effect of the chemotactic gradient, $\omega_{rep}$, $\omega_{adh}$ are respectively the coefficients of the repulsion and the adhesion effect between immune cells (6), $\beta$ is the coefficient of the alignment term, $\omega_{rep}^{tum}$ is the coefficient of the repulsion effect between immune and cancer cells and $h$ is the bandwidth of the kernel density for the immune cells. The other parameters are fixed according to Table 1.

For fixed $N_t$ time steps, the numerical solution to (20) is compared to the reference one when computing the value of $\Theta_{opt}$ minimizing the objective function (32). In the following tests, we set $N_t = 5$ and $\lambda^2 = 10^{-6}$ for the regularization term.

The optimization starts from values in $\Theta_0$ given by the value of the parameters at the microscopic scale. Then the set of admissible value $\mathcal{A}$ is chosen as $[0, 50\Theta_0]$, taking into account physical conditions as the non-negativity of the coefficients and the fact that, by modeling construction, $R_{adh} > R_{rep}$.
<table><caption>Table 1: Table of parameters values considered to generate synthetic data for Test1 - Test3.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td>R<sub>imm</sub></td><td>Radius of ICs</td><td>2 &middot; 10<sup>-2</sup></td></tr><tr><td>R<sub>tum</sub></td><td>Radius of TCs</td><td>5 &middot; 10<sup>-2</sup></td></tr><tr><td>&alpha;</td><td>Damping Coefficient</td><td>1 &middot; 10<sup>2</sup></td></tr><tr><td>&kappa;</td><td>Degradation Rate of &phi;</td><td>2 &middot; 10<sup>-1</sup></td></tr><tr><td>&xi;</td><td>Production Rate of &phi;</td><td>1 &middot; 10<sup>3</sup></td></tr><tr><td>D</td><td>Diffusion Coefficient of &phi;</td><td>4.5 &middot; 10<sup>1</sup></td></tr><tr><td>R<sub>adh</sub></td><td>Radius of action of Adhesion between ICs</td><td>6 &middot; 10<sup>-2</sup></td></tr><tr><td>R<sub>rep</sub></td><td>Radius of action of Repulsion between ICs</td><td>4 &middot; 10<sup>-2</sup></td></tr><tr><td>R<sup>tum</sup><sub>rep</sub></td><td>Radius of action of Repulsion between ICs and TCs</td><td>7 &middot; 10<sup>-2</sup></td></tr><tr><td>R<sub>align</sub></td><td>Radius of action of Alignment between ICs</td><td>6 &middot; 10<sup>-2</sup></td></tr></tbody></table>

**Test 1: neglecting repulsion interaction between immune and cancer cells** In this framework, we do not consider mechanical interactions between cancer and immune cells, setting $\omega_{rep}^{tum} = 0$ in (15) and in (20). We set parameters value as in Table 1 and

$$ \eta = 6, \quad \omega_{rep} = 5 \cdot 10^2, \quad \omega_{adh} = 4, \quad \beta = 2 \cdot 10^3, \quad h = 1.2, \qquad (34) $$

to compute the trajectories of immune cells. Hence we run the optimization procedure, setting the components of $\Theta_0$ as the values in (34). The optimization procedure gives the following optimal solution

$$ [71.62 \quad 6.74 \cdot 10^2 \quad 39.43 \quad 1.68 \cdot 10^3 \quad 0 \quad 2.45], \qquad (35) $$

with estimation error

$$ E = 6.73 \cdot 10^{-2}. \qquad (36) $$

Fig. 7 shows five plots of the numerical simulation of (20) with the estimated parameters at different times. Positions of the immune cells, generating the reference synthetic data, are marked as black disks. We observe that, also from a qualitative point of view, a good level of agreement is achieved between the two scales. In particular, the effect of the chemical signal produced by the three cancer cells plays a crucial role. At the microscopic scale, cells are chemically attracted by the chemical signal, and since $\omega_{rep}^{tum} = 0$, they overlap cancer cells. On the microscopic scale, the density gets higher in the surrounding of $\mathbf{y}_1, \mathbf{y}_2, \mathbf{y}_3$.

**Test 2: the effect of alignment** With respect to attraction and repulsion kind of interactions, the role of alignment in the dynamics is less evident. For this reason, we simulate the same scenario of Test 1 but neglecting alignment, hence assuming $\beta = 0$ both at the microscopic level in (16) and in the macroscopic one (20). Consequently, we set again the components of $\Theta_0$ as the values in (34), with $\beta = 0$. We obtain

$$ \Theta_{opt} = [81.3 \quad 7.52 \cdot 10^2 \quad 39.2 \quad 0 \quad 0 \quad 2.4], \qquad (37) $$

with

$$ E = 6.01 \cdot 10^{-2}. \qquad (38) $$

In Fig. 8 we show the comparison of the macroscopic densities computed respectively with (35) and (37) at final time. In the literature of hybrid models, it has already been
<figure><img src="image_8.png" /><figcaption>Figure 7: Numerical simulation of (20) with estimated parameters. (a)-(e) Plots at different times of the (macroscopic) density $\rho$ versus (microscopic) positions of the cells, represented by black circles.</figcaption></figure>

observed that alignment plays a crucial role in the context of Cancer-on-chip modeling [8, 21] at the microscopic scale. Due to the presence of alignment, the convergence of immune cells towards cancer cells less strong, allowing cells to align with others in the surroundings (Fig. a). Neglecting alignment, immune cells are trapped by cancer cells (Fig. b), almost clustering around them. Black arrows are used to plot the momentum $\rho u$ evaluated in each point of the numerical grid, highlighting the behavior of the system at the macroscopic scale.

**Test 3: repulsion between immune and cancer cells** In the third scenario, we consider mechanical interactions between immune and cancer cells. Despite the lack of experimental data, this could be considered a more realistic setting, in which a repulsion effect among cancer cells and immune cells occurs at a distance between the centers less than $R_{rep}^{tum}$, avoiding non feasible overlapping. In order to generate synthetic data, we run the microscopic model (16) starting with the same initial condition of the previous tests, in the same parameter setting of Table 1, adding a non-null value for $\omega_{rep}^{tum}$. Hence we start the optimization procedure with

$$ \eta = 6, \quad \omega_{rep} = 5 \cdot 10^2, \quad \omega_{adh} = 4, \quad \beta = 2 \cdot 10^3, \quad \omega_{rep}^{tum} = 8.5 \cdot 10^2 \quad h = 1.2, \quad (39) $$

obtaining

$$ \Theta_{opt} = [80.2 \quad 1.32 \cdot 10^3 \quad 33.6 \quad 4.27 \cdot 10^3 \quad 2.31 \cdot 10^2 \quad 2.5], \quad (40) $$

with

$$ E = 9 \cdot 10^{-2}. \quad (41) $$
<figure><img src="image_9.png" /><figcaption>Figure 8: Comparison at final time of numerical simulations of (20) with estimated parameters assuming a) $\mathcal{I} = \mathcal{I}_1 + \mathcal{I}_2 + \mathcal{I}_3$, b) $\mathcal{I} = \mathcal{I}_2 + \mathcal{I}_3$. Black arrows describing the momentum $\rho u$, are plotted against the plots of the density $\rho$.</figcaption></figure>

Fig. 9 shows different plots of the simulations with the estimated parameters. With respect to Fig. 7, we observe the effect of the repulsion between cancer and immune cells: the density plot is less concentrated on tumor cells, plotted as blue disks endowed with radius $R_{tum}$, and immune cells do no overlap to them.

We remark that, in previous works of the literature of the related field [9, 15], the estimation of the parameters involved in the considered models is performed against synthetic data produced by the model itself. Hence, it reasonable to expect a better agreement in term of relative error between the parameters used to generate synthetic data and the optimal parameters obtained through the minimization process. The issue we here address introduces a higher degree of difficulty, since synthetic data are generated through the hybrid model (16), which acts on a different scale compared to the model under investigation (20). The numerical results here obtained show that the non-local pressureless Euler model (20) represents a valid counterpart for (16) on the macroscopic scale, with the advantages of a clear link to the underlying microscopic dynamics.

**Sensitivity Analysis** To conclude our analysis, we evaluate the sensitivity of the model (20), following the approach in [15,16]. Starting with the optimal values found in (40), we evaluate the effects of variations of the model parameters in (33) on the maximal density of immune cells at final time. Denoting with $\mathcal{Y}(\Theta) = \mathcal{Y}(\Theta_1, ..., \Theta_6)$ its value running the model with parameters given in $\Theta$, the sensitivity of the output to any parameter $\Theta_i$ is obtained as

$$S = \frac{|\mathcal{Y}(\Theta_1, ..., \Theta_i \pm \delta, ..., \Theta_6) - \mathcal{Y}(\Theta_1, ..., \Theta_6)|}{\mathcal{Y}(\Theta_1, ..., \Theta_6)} \frac{\Theta_i}{\delta}.$$

Here we consider a 5% perturbation of the parameters choosing $\delta = 0.05\Theta_i$. The results are shown in Table2. As expected, parameters $\eta$, $w_{rep}$ and $w_{rep}^{tum}$ play a significant role in the cells density. Compared to the others, parameter $w_{adh}$ tuning attraction mildly influences the output. This is reasonable, since chemotaxis already plays an adhesive role in the dynamics. The alignment parameter $\beta$ appears to have a small influence on the maximum density value, confirming the results obtained in Test2.
<figure><img src="image_10.png" /><figcaption>Figure 9: Numerical simulation of (20) with estimated parameters. (a)-(e) Plots at different times of the (macroscopic) density $\rho$ versus (microscopic) positions of the immune and tumor cells, represented by black and blue circles, respectively.</figcaption></figure>

<table><caption>Table 2: Sensitivity study with $\delta = 5\%$ in the simulation set of Test3 with reference values given in (40). $S_y$ is the sensitivity considering as model output the maximum density of immune cell at final time.</caption><thead><tr><th>Parameter</th><th>S<sub>y</sub></th></tr></thead><tbody><tr><td>&eta; + &delta;</td><td>1.0978</td></tr><tr><td>&eta; - &delta;</td><td>1.0765</td></tr><tr><td>w<sub>rep</sub> + &delta;</td><td>0.4238</td></tr><tr><td>w<sub>rep</sub> - &delta;</td><td>0.4545</td></tr><tr><td>w<sub>adh</sub> + &delta;</td><td>0.1084</td></tr><tr><td>w<sub>adh</sub> - &delta;</td><td>0.1067</td></tr><tr><td>&beta; + &delta;</td><td>0.0012</td></tr><tr><td>&beta; - &delta;</td><td>0.0012</td></tr><tr><td>w<sup>tum</sup><sub>rep</sub> + &delta;</td><td>0.4310</td></tr><tr><td>w<sup>tum</sup><sub>rep</sub> - &delta;</td><td>0.4843</td></tr><tr><td>h + &delta;</td><td>0.0124</td></tr><tr><td>h - &delta;</td><td>0.0106</td></tr></tbody></table>

# 5 Conclusions

In this paper we numerically investigate the pressureless Euler-type model rigorously derived in [42] as macroscopic counterpart for the class of hybrid coupled models. The derivation requires strong assumptions on the initial data, which are not satisfied in real-world applications. Starting with a comparison with previous models of the literature,
including a phenomenological pressure term, we observe the capability of this new model in reproducing typical patterns of collective cell migration. The investigation is complemented by a quantitative calibration against data collected in Cancer-On-Chip settings. To this aim, we extend the model in different directions: including two populations of cells, with different behaviors, and including attraction-repulsion kind of interactions in addition to alignment. Extending the model to a multi-population framework enables a more realistic representation of the phenomenon, as it captures the distinct roles of various immune cell types in collective dynamics. All the numerical simulations are performed in a multidimensional setting. Starting from synthetic data generated at the microscopic scale, the parameter estimation procedure in three different two-dimensional scenarios shows a good level of agreement between the two scales. The availability of new experimental data, with (macroscopic) information concerning the chemical concentration, will allow to further exploit the potential of the proposed approach. From a theoretical point of view, the non-locality of the model and the absence of a pressure term still represent challenging aspects for theoretical studies, leaving space for further investigations.

## Acknowledgments

M. Menci and T.Tenna are members of GNCS-INDAM research group, R. Natalini is member of GNAMPA-INDAM research group. M. Menci is funded by INdAM - GNCS Project, CUP E53C24001950001, entitled “MODA: Integrating MOdel-based and DAta-Driven Methods for Multiscale Biological Systems”. T. Tenna is funded by the European Union’s Horizon Europe research and innovation programme under the Marie Skłodowska-Curie Doctoral Network DataHyking (Grant No.101072546).

## A Numerical Scheme for macroscopic models

Let us consider the Euler-type systems coupled with chemotaxis studied in this paper in a two-dimensional bounded domain $\Omega \subseteq \mathbb{R}^2$, rewritten in the following compact form

$$ \begin{cases} \partial_t w + \partial_{x_1} A_1(w) + \partial_{x_2} A_2(w) = F(w, \varphi), \\ \partial_t \varphi = D \Delta \varphi + a\rho - b\varphi, \end{cases} \qquad (42) $$

where $w = (\rho, \rho u_1, \rho u_2)$ and

$$ A_1(w) = \begin{pmatrix} \rho u_1 \\ \rho u_1^2 + \varepsilon P(\rho) \\ \rho u_1 u_2 \end{pmatrix}, \quad A_2(w) = \begin{pmatrix} \rho u_2 \\ \rho u_1 u_2 \\ \rho u_2^2 + \varepsilon P(\rho) \end{pmatrix}. \qquad (43) $$

and

$$ F(w, \varphi) = \begin{pmatrix} 0 \\ (1 - \varepsilon)\rho \int \gamma(\cdot - y, u_1(\cdot) - u_1(y))\rho(y)dy + \eta\rho\partial_{x_1}\varphi - \alpha\rho u_1 \\ (1 - \varepsilon)\rho \int \gamma(\cdot - y, u_2(\cdot) - u_2(y))\rho(y)dy + \eta\rho\partial_{x_2}\varphi - \alpha\rho u_2 \end{pmatrix}. \qquad (44) $$

The choice $\varepsilon = 1$ leads to (7), whereas the choice $\varepsilon = 0$ leads to the non-local pressureless Euler type model (11).
**Numerical Approximation of the Euler-type system** In order to numerically solve (7) and (11), we consider a class of numerical schemes based on a discrete kinetic approximation for multidimensional systems of conservation laws, following [3]. This approach enables us to approximate a general quasi-linear conservation laws of the form

$$ \partial_t w + \nabla_x A(w) = 0, \qquad (45) $$

by a sequence of semi-linear hyperbolic systems, through a BGK relaxation equation, leading to

$$ \partial_t f^\varepsilon + \sum_{j=1}^d \Lambda_j \partial_{x_j} f^\varepsilon = \frac{1}{\varepsilon} (M(w^\varepsilon) - f^\varepsilon), \qquad (46) $$

with Cauchy initial condition $f^\varepsilon(x, 0) = f_0^\varepsilon(x)$ and by setting

$$ u^\varepsilon := \sum_{i=1}^l f_i^\varepsilon. \qquad (47) $$

Here $\varepsilon$ is a positive number, $\Lambda_j$ are real diagonal matrices $l \times l$, $f^\varepsilon = (f_1^\varepsilon, \dots, f_l^\varepsilon)$ and $\mathcal{M}: \mathbb{R}^k \to \mathbb{R}^l$ is a Lipschitz continuous function.

The numerical schemes presented in [3] are constructed by splitting (46) into a homogeneous linear part and an ordinary differential system, which is exactly solved thanks to the particular structure of the source term.

The main advantage of this approximation, especially in the multidimensional case, is the possibility of avoiding the resolution of local Riemann problems in the design of numerical schemes. In this framework, the scalar and the system case are treated in the same way at numerical level and, since all the approximating problems are in diagonal form, the numerical scheme is easily implemented.

Here we detail the construction of the scheme for our model (42). Let us take $\Omega = [0, L] \times [0, L]$ and let use denote $\Delta x$ the space step. We consider the discretization points $\mathbf{x}_\alpha = (\alpha_1 \Delta x, \alpha_2 \Delta x), 0 \le \alpha_i \le N_x + 1$. We denote the time step by $\Delta t$ and the approximation of a function $f$ at time $t_n = n\Delta t$ as $f^n$. Let us introduce a 5-velocities relaxation scheme with

$$ \lambda_1 = \lambda(1, 0), \quad \lambda_2 = \lambda(0, 1), \quad \lambda_3 = \lambda(-1, 0), \quad \lambda_4 = \lambda(0, -1), \quad \lambda_5 = (0, 0), \qquad (48) $$

for some $\lambda > 0$. The time and space steps have to satisfy the CFL stability condition, setting $\Delta t = 0.9\Delta x/\lambda$.

Now we introduce the corresponding Maxwellians $M_i(w) \in \mathbb{R}^3, i = 1, \dots, 5$

$$ M_i(w) = a_i w + b_{i1} A_1(w) + b_{i2} A_2(w), \qquad (49) $$

where

$$ \begin{gathered} a_1 = \dots = a_4 = a, a_5 = 1 - 4a, \quad a < \frac{1}{4}, \qquad (50) \\ b_{11} = b_{22} = -b_{31} = -b_{42} = \frac{1}{2\lambda}, \quad b_{ij} = 0 \quad \text{otherwise,} \end{gathered} $$

in order to obtain consistent conditions for the convergence of the kinetic solution to the macroscopic one. Let us denote by $w^{n,\alpha}$ the approximation of $w$ at the point $\mathbf{x}_\alpha \in \mathbb{R}^2$ and at time $t_n$. We set the discretization of Maxwellians as

$$ f_i^{n,\alpha} = M_i(w^{n,\alpha}), \quad \text{for } i = 1, \dots, 5. \qquad (51) $$
We evolve each of the functions $f_i$, $1 \le i \le 5$, in time through an upwind scheme by following the velocity $\lambda_i$, as

$$f^{n+\frac{1}{2},\alpha} = H^{up}(\Delta t)(f_i^{n,\alpha}), \qquad (52)$$

where

$$H^{up}(\Delta t)(f_i^{n,\alpha}) = f_i^{n,\alpha} - \mu \sum_{j=1}^{2} \lambda_{ij}(f_i^{n,\alpha+e_j} - f_i^{n,\alpha-e_j}) \\ + \mu \sum_{j=1}^{2} |\lambda_{ij}|(f_i^{n,\alpha+e_j} - 2f_i^{n,\alpha} + f_i^{n,\alpha-e_j}), \quad (53)$$

and $\mu = \Delta t/(2\Delta x)$. This upwind scheme could be improved through a flux limiter approach, in order to reduce the diffusion and obtain a second order approximation far from oscillations. Therefore, the new scheme reads as (52), where $H^{up}$ is replaced by

$$H^S(\Delta t) = H^{up}(\Delta t) - \frac{1}{2}\mu \sum_{j=1}^{2} \lambda_{ij} (\mathrm{sign}(\lambda_{ij}) - \mu\lambda_{ij}) \\ \left( \phi_{i,j+\frac{1}{2}} \Delta_{j+\frac{1}{2}} f_i^{n,\alpha} - \phi_{i,j-\frac{1}{2}} \Delta_{j-\frac{1}{2}} f_i^{n,\alpha} \right), \quad (54)$$

where we have introduced the following quantities

$$\phi_{i,j+\frac{1}{2}} := \phi \left( \frac{\Delta_{j-\frac{1}{2}} f_i^{n,\alpha}}{\Delta_{j+\frac{1}{2}} f_i^{n,\alpha}} \right), \quad \Delta_{j+\frac{1}{2}} f_i^{n,\alpha} = f_i^{n,\alpha+e_j} - f_i^{n,\alpha}, \qquad (55)$$

in which $\phi$ is the minmod function defined as $\phi(r) = \max\{0, \min\{r, 1\}\}$. Finally, we just conclude by setting

$$w^{n+1} = \sum_{i=1}^{5} f_i^{n+\frac{1}{2}} + F(w^n). \qquad (56)$$

The $F$ is the source term approximated in explicit as

$$F(w^n) = \begin{pmatrix} 0 \\ \rho^n \mathcal{I}^n + \eta \rho^n D_{x_1} \varphi^n - \alpha \rho^n u_1^n \\ \rho^n \mathcal{I}^n + \eta \rho^n D_{x_2} \varphi^n - \alpha \rho^n u_2^n \end{pmatrix}, \qquad (57)$$

where $\mathcal{I}^n$ is the approximated value of the integral term obtained through a quadrature formula, while $D_{x_i}$ is a three points discretization of the operator $\partial_{x_i} \varphi$.

**Numerical Approximation of the Parabolic Equation** For the discretization of the parabolic equation of the chemoattractant in (42), we apply a classical central finite difference scheme with 5-point stencil for the Laplacian and the $\theta$-method scheme in time. The discretization of the Laplacian for a function $\varphi$ becomes:

$$\Delta\varphi(x_k, y_m) = \frac{\varphi_{k+1,m} - 2\varphi_{k,m} + \varphi_{k-1,m}}{\Delta x^2} + \frac{\varphi_{k,m+1} - 2\varphi_{k,m} + \varphi_{k,m-1}}{\Delta y^2}, \quad (58)$$
while the θ-method for the discretization in time leads to the final scheme for the parabolic equation:

$$ \frac{\varphi_{k,m}^{n+1} - \varphi_{k,m}^{n}}{\Delta t} = D\theta \left( \frac{\varphi_{k+1,m}^{n+1} - 2\varphi_{k,m}^{n+1} + \varphi_{k-1,m}^{n+1}}{\Delta x^2} + \frac{\varphi_{k,m+1}^{n+1} - 2\varphi_{k,m}^{n+1} + \varphi_{k,m-1}^{n+1}}{\Delta y^2} \right) \\ + D(1-\theta) \left( \frac{\varphi_{k+1,m}^{n} - 2\varphi_{k,m}^{n} + \varphi_{k-1,m}^{n}}{\Delta x^2} + \frac{\varphi_{k,m+1}^{n} - 2\varphi_{k,m}^{n} + \varphi_{k,m-1}^{n}}{\Delta y^2} \right) \\ + \theta s_{k,m}^{n+1} + (1-\theta)s_{k,m}^{n}, $$

where $s_{k,m}^{n} = s(x_k, y_m, t_n)$ is an external heat source. The numerical scheme could be complemented with suitable boundary conditions.

The advantage is achieved by choosing $\theta = \frac{1}{2}$, usually known as *Crank-Nicolson scheme*, which increases the approximation order in time to 2, benefiting from unconditional stability.

## References

[1] Giacomo Albi, Young-Pil Choi, and Axel-Stefan Haeck. Pressureless Euler alignment system with control. *Mathematical Models and Methods in Applied Sciences*, 28(09):1635-1664, 2018.

[2] D Ambrosi, F Bussolino, and L Preziosi. A review of vasculogenesis models. *Journal of Theoretical Medicine*, 6(1):1-19, 2005.

[3] Denise Aregba-Driollet and Roberto Natalini. Discrete kinetic schemes for multi-dimensional systems of conservation laws. *SIAM Journal on Numerical Analysis*, 37(6):1973-2004, 2000.

[4] François Bouchut and François James. Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness. *Commun. Partial Differ. Equations*, 24(11-12):2173-2189, 1999.

[5] François Bouchut, Shi Jin, and Xiantao Li. Numerical approximations of pressureless and isothermal gas dynamics. *SIAM J. Numer. Anal.*, 41(1):135-158, 2003.

[6] Elishan Christian Braun. *Organs-On-Chips: mathematical modelling and parameter estimation*. PhD Thesis, Università degli studi Roma Tre, 2021.

[7] Yann Brenier and Emmanuel Grenier. Sticky particles and scalar conservation laws. *SIAM J. Numer. Anal.*, 35(6):2317-2328, 1998.

[8] Gabriella Bretti, Elio Campanile, Marta Menci, and Roberto Natalini. A scenario-based study on hybrid PDE-ODE model for cancer-on-chip experiment. In *Problems in Mathematical Biophysics: A Volume in Memory of Alberto Gandolfi*, pages 37-57. Springer, 2024.

[9] Gabriella Bretti, Adele De Ninno, Roberto Natalini, Daniele Peri, and Nicole Roselli. Estimation Algorithm for a Hybrid PDE-ODE Model Inspired by Immunocompetent Cancer-on-Chip Experiment. *Axioms*, 10(4):243, 2021.
[10] José A Carrillo, Young-Pil Choi, and Sergio P Perez. A review on attractive-repulsive hydrodynamics for consensus in collective behavior. *Active Particles, Volume 1: Advances in Theory, Models, and Applications*, pages 259-298, 2017.

[11] José A Carrillo, Young-Pil Choi, Eitan Tadmor, and Changhui Tan. Critical thresholds in 1D Euler equations with non-local forces. *Mathematical Models and Methods in Applied Sciences*, 26(01):185-206, 2016.

[12] José A Carrillo, Stephan Martin, and Vladislav Panferov. A new interaction potential for swarming models. *Physica D: Nonlinear Phenomena*, 260:112-126, 2013.

[13] Fausto Cavalli, Andrea Gamba, Giovanni Naldi, and Matteo Semplice. Approximation of 2D and 3D models of chemotactic cell movement in vasculogenesis. *Math Everywhere: Deterministic and Stochastic Modelling in Biomedicine, Economics and Industry. Dedicated to the 60th Birthday of Vincenzo Capasso*, pages 179-191, 2007.

[14] Yen-Chi Chen. A tutorial on kernel density estimation and recent advances. *Bio-statistics & Epidemiology*, 1(1):161-187, 2017.

[15] Giorgia Ciavolella, Nathalie Ferrand, Michèle Sabbah, Benoît Perthame, and Roberto Natalini. A model for membrane degradation using a gelatin invadopodia assay. *Bulletin of Mathematical Biology*, 86(3):30, 2024.

[16] Fabrizio Clarelli, Cristiana Di Russo, Roberto Natalini, and Magali Ribot. A fluid dynamics multidimensional model of biofilm growth: stability, influence of environment and sensitivity. *Mathematical medicine and biology: a journal of the IMA*, 33(4):371-395, 2016.

[17] Iain D Couzin, Jens Krause, Nigel R Franks, and Simon A Levin. Effective leadership and decision-making in animal groups on the move. *Nature*, 433(7025):513-516, 2005.

[18] Felipe Cucker and Steve Smale. Emergent behavior in flocks. *IEEE Transactions on Automatic Control*, 52(5):852-862, 2007.

[19] Felipe Cucker and Steve Smale. On the mathematics of emergence. *Japanese Journal of Mathematics*, 2:197-227, 2007.

[20] Ezio Di Costanzo, Alessandro Giacomello, Elisa Messina, Roberto Natalini, Giuseppe Pontrelli, Fabrizio Rossi, Robert Smits, and Monika Twarogowska. A discrete in continuous mathematical model of cardiac progenitor cells formation and growth as spheroid clusters (Cardiospheres). *Mathematical Medicine and Biology: a Journal of the IMA*, 35(1):121-144, 2018.

[21] Ezio Di Costanzo, Marta Menci, Eleonora Messina, Roberto Natalini, and Antonia Vecchio. A hybrid model of collective motion of discrete particles under alignment and continuum chemotaxis. *Discrete & Continuous Dynamical Systems-Series B*, 25(1):443-472, 2020.

[22] Ezio Di Costanzo, Roberto Natalini, and Luigi Preziosi. A hybrid mathematical model for self-organizing cell migration in the zebrafish lateral line. *Journal of Mathematical Biology*, 71(1):171-214, 2015.
[23] Maria R D'Orsogna, Yao-Li Chuang, Andrea L Bertozzi, and Lincoln S Chayes. Self-propelled particles with soft-core interactions: patterns, stability, and collapse. *Physical Review Letters*, 96(10):104302, 2006.

[24] Heinz W. Engl, Martin Hanke, and Andreas Neubauer. *Regularization of inverse problems*, volume 375 of *Math. Appl., Dordr*. Dordrecht: Kluwer Academic Publishers, 1996.

[25] Alessio Figalli and Moon-Jin Kang. A rigorous derivation from the kinetic Cucker-Smale model to the pressureless Euler system with nonlocal alignment. *Analysis & PDE*, 12(3):843–866, 2018.

[26] A. Gamba, D. Ambrosi, A. Coniglio, A. de Candia, S. Di Talia, E. Giraudo, G. Serini, L. Preziosi, and F. Bussolino. Percolation, Morphogenesis, and Burgers Dynamics in Blood Vessels Formation. *Phys. Rev. Lett.*, 90:118101, Mar 2003.

[27] Thomas Hillen and Kevin J Painter. A user's guide to PDE models for chemotaxis. *Journal of Mathematical Biology*, 58(1-2):183, 2009.

[28] Evelyn F Keller and Lee A Segel. Initiation of slime mold aggregation viewed as an instability. *Journal of Theoretical Biology*, 26(3):399-415, 1970.

[29] Evelyn F Keller and Lee A Segel. Model for chemotaxis. *Journal of theoretical biology*, 30(2):225-234, 1971.

[30] R Kowalczyk, A Gamba, and L Preziosi. On the stability of homogeneous solutions to some aggregation models. *Discrete and Continuous Dynamical Systems Series B*, 4(1):203–220, 2004.

[31] Marta Menci, Roberto Natalini, and Thierry Paul. Microscopic, kinetic and hydrodynamic hybrid models of collective motions with chemotaxis: a numerical study. *Mathematics and Mechanics of Complex Systems*, 12(1):47–83, 2023.

[32] Marta Menci and Marco Papi. Global solutions for a path-dependent hybrid system of differential equations under parabolic signal. *Nonlinear Analysis*, 184:172–192, 2019.

[33] Marta Menci and Marco Papi. Existence of solutions for hybrid systems of differential equations under exogenous information with discontinuous source term. *Nonlinear Analysis*, 221:112885, 2022.

[34] Marta Menci, Marco Papi, Maria Michaela Porzio, and Flavia Smarrazzo. On a coupled hybrid system of nonlinear differential equations with a nonlocal concentration. *Journal of Differential Equations*, 361:288–338, 2023.

[35] Marta Menci, Marco Papi, and Flavia Smarrazzo. Existence of solutions for coupled hybrid systems of differential equations for microscopic dynamics and local concentrations. *Communications on Pure and Applied Analysis*, 22(7):2146–2168, 2023.

[36] Sebastien Motsch and Eitan Tadmor. A new model for self-organized dynamics and its flocking behavior. *Journal of Statistical Physics*, 144:923–947, 2011.
[37] James D Murray. *Mathematical biology: I. An introduction*, volume 17. Springer Science & Business Media, 2007.

[38] EA Nadaraya. On non-parametric estimates of density functions and regression curves. *Theory of Probability & Its Applications*, 10(1):186–190, 1965.

[39] Vidyanand Nanjundiah. Chemotaxis, signal relaying and aggregation morphology. *Journal of Theoretical Biology*, 42(1):63–105, 1973.

[40] Roberto Natalini. A discrete kinetic approximation of entropy solutions to multidimensional scalar conservation laws. *Journal of Differential Equations*, 148(2):292–317, 1998.

[41] Roberto Natalini and Thierry Paul. On the mean field limit for Cucker-Smale models. *Discrete and Continuous Dynamical Systems-B*, 27(5):2873–2889, 2022.

[42] Roberto Natalini and Thierry Paul. The mean-field limit for hybrid models of collective motions with chemotaxis. *SIAM Journal on Mathematical Analysis*, 55(2):900–928, 2023.

[43] Roberto Natalini, Magali Ribot, and Monika Twarogowska. A numerical comparison between degenerate parabolic and quasilinear hyperbolic models of cell movements under chemotaxis. *Journal of Scientific Computing*, 63:654–677, 2015.

[44] Kevin J Painter. Mathematical models for chemotaxis and their applications in self-organisation phenomena. *Journal of Theoretical Biology*, 481:162–182, 2019.

[45] Benoît Perthame. *Transport equations in biology*. Springer Science & Business Media, 2006.

[46] M. Scianna, C. G. Bell, and Luigi Preziosi. A review of mathematical models for the formation of vascular networks. *J. Theor. Biol.*, 333:174–209, 2013.

[47] Guido Serini, Davide Ambrosi, Enrico Giraudo, Andrea Gamba, Luigi Preziosi, and Federico Bussolino. Modeling the early stages of vascular network assembly. *The EMBO journal*, 22(8):1771–1779, 2003.

[48] Balint Szabo, GJ Szöllösi, B Gönci, Zs Jurányi, David Selmeczi, and Tamás Vicsek. Phase transition in the collective migration of tissue cells: experiment and model. *Physical Review E*, 74(6):061908, 2006.

[49] Andrea Tosin, Davide Ambrosi, and Luigi Preziosi. Mechanics and chemotaxis in the morphogenesis of vascular networks. *Bulletin of Mathematical Biology*, 68:1819–1836, 2006.