Phenotype-structuring of non-local kinetic models of cell
migration driven by environmental sensing

Tommaso Lorenzi*

Nadia Loy†

Chiara Villa‡§

April 25, 2025

Abstract

The capability of cells to form surface extensions to non-locally probe the surrounding en-
vironment plays a key role in cell migration. The existing mathematical models for migration
of cell populations driven by this non-local form of environmental sensing rely on the simpli-
fying assumption that cells in the population share the same cytoskeletal properties, and thus
form surface extensions of the same size. To overcome this simplification, we develop a kinetic
modelling framework wherein a population of migrating cells is structured by a continuous
phenotypic variable that captures variability in structural properties of the cytoskeleton. This
framework provides a multiscale representation of cell migration, from single-cell dynamics to
population-level behaviours, as we start with a microscopic model that describes the dynam-
ics of single cells in terms of stochastic processes. Next, we formally derive the mesoscopic
counterpart of this model, which consists of a phenotype-structured kinetic equation that
features a phenotype-dependent non-locality. Then, considering an appropriately rescaled
version of this kinetic equation, we formally derive the corresponding macroscopic model,
which takes the form of a partial differential equation for the cell number density. To validate
the formal procedures employed to derive the macroscopic model from the microscopic model,
through the mesoscopic one, we first compare the results of numerical simulations of the two
models. We then compare numerical solutions of the macroscopic model with the results of
cell locomotion assays, to test the ability of the model to recapitulate qualitative features of
experimental observations.

# 1 Introduction

Biological background

Cell migration is at the heart of morphogenesis, is essential to ensure successful wound healing, immune response, and tissue homeostasis in adult organisms, and is also central to the progression of different pathologies, such as the metastatic cascade associated with cancer malignancy (Trepat et al., 2012).

<footnote>*Department of Mathematical Sciences “G. L. Lagrange”, Politecnico di Torino, 10129 Torino, Italy (tommaso.lorenzi@polito.it)</footnote>

<footnote>†Department of Mathematical Sciences “G. L. Lagrange”, Politecnico di Torino, 10129 Torino, Italy (na-dia.loy@polito.it)</footnote>

<footnote>‡Sorbonne Université, CNRS, Université de Paris, Inria, Laboratoire Jacques-Louis Lions UMR7598, F-75005 Paris, France</footnote>

<footnote>§Université Paris-Saclay, Inria, Centre Inria de Saclay, 91120, Palaiseau, France (chiara.villa@inria.fr)</footnote>
A key role in cell migration is played by the capability of cells to form surface extensions, which enable them to probe the surrounding environment even multiple cell diameters away (i.e. to perform non-local environmental sensing) (Friedl and Gilmour, 2009; Friedl et al., 1998; Mitchison and Cramer, 1996). The size and shape of such surface extensions depend on the cytoskeletal properties of the cells, which have been shown to vary not only between cell populations but also amongst cells within the same population, leading to intra-population heterogeneity in cell locomotive phenotype. For instance, fibroblasts in wound healing have been observed to move following the extension of short or long lamellipodia (Trepat et al., 2012); invading tumour cells have been reported to display different morphologies and migratory abilities depending on whether they are more epithelial-like or mesenchymal-like (Jolly et al., 2018); and it has been demonstrated that myoblasts migrate via the extension of long pseudopods or short leading lamellae (Goodman et al., 1989).

Changes in the cell locomotive phenotype can occur following external stimuli, being these chemical or mechanical. For instance: fibroblasts involved in wound healing may increase their lamellipodia size and speed of migration in response to different growth factors (Trepat et al., 2012); epithelial-to-mesenchymal transition in cancer can be induced by both growth factors and increased stiffness of the extracellular matrix (ECM) (Gkretsi and Stylianopoulos, 2018); and the locomotive phenotype of myoblasts has been shown to depend on the adhesive glycoproteins present on the substrate on which they move, likely due to mechanotransduction (Givant-Horwitz et al., 2005; Goodman et al., 1989; Miyamoto et al., 1998; Rousselle and Scoazec, 2020). Independently of the specific pathways responsible for such changes, it is clear that the cytoskeletal structure is not a binary state and this translates into a range of different locomotive phenotypes, which are linked to the characteristics of the surface extensions that are formed by the cells.

In the light of mounting evidence indicating that myoblast cytodifferentiation is deeply affected by the molecular composition of the surrounding ECM, myoblasts locomotion over different ECM components was examined in Goodman et al. (1989) by means of *in vitro* stripe assays. In each assay, parallel stripes of two different glycoproteins (i.e. laminin and fibronectin) were absorbed onto the substrates. Myoblasts (i.e. murine myoblasts, cell line MM14dy) were plated at similar densities at one end of the stripes and were allowed to migrate. Cell migration was then tracked via videomicroscopy. Over laminin stripes, cells migrated generally more rapidly undergoing consecutive cycles of long-pseudopod extension and rapid translocation following release from the substrate. On the other hand, cells moved more slowly on fibronectin stripes (about three times slower than over laminin) and mainly through the formation of short leading lamellae. In Goodman et al. (1989) it was noted that, along with different motile behaviours, cells adhering to different ECM components displayed different cytoskeletal characteristics. In myoblasts migrating over fibronectin, actin<fnref n="1" /> and α-actinin<fnref n="2" /> quickly organised into stress fibres and vinculin<fnref n="3" /> organised into focal contacts, while these proteins remain more sparsely distributed within the cytoskeleton of myoblasts migrating over laminin. Whilst the specific pathways through which cell mechanosensing translates into cytoskeletal changes and different locomotive strategies are not discussed in Goodman et al. (1989), the experimental results therein presented highlight the important role that intra-population heterogeneity and environment-induced changes in locomotive phenotype play in cell migration driven by non-local environmental sensing.

## Mathematical modelling background and content of the article

A variety of mathematical models for cell migration driven by non-local environmental sensing have been proposed in the last twenty years – e.g. see the review (Chen et al., 2020) and references therein. The majority of these models rely on the assumption that cells in the population share the same cytoskeletal properties, which do not evolve in time, and thus sense the surrounding environment by means of surface extensions of the same fixed size. In view of the above biological

<footnote n="1">¹Actin is a family of globular multi-functional proteins that form microfilaments in the cytoskeleton.</footnote>

<footnote n="2">²α-actinin is an actin-binding protein.</footnote>

<footnote n="3">³Vinculin is a membrane-cytoskeletal protein in focal adhesion plaques, which is involved in linkage of integrin adhesion molecules to the actin cytoskeleton.</footnote>
background, this is clearly a simplification that limits the domain of application of these models and their outputs.

To overcome such a simplification, building upon previous work on phenotype-structured models of cell movement reviewed in Lorenzi et al. (2024b), in this article we generalise the non-local kinetic modelling approach of Loy and Preziosi (2020a,b) by developing a modelling framework wherein a population of migrating cells is structured by a continuous phenotypic variable that captures variability in structural properties of the cytoskeleton. Cells with different values of the phenotypic variable form surface projections of different length and can then probe the surrounding environment over regions of different sizes. As a result, in contrast to Loy and Preziosi (2020a,b), the integration domain for the non-locality in the model considered here is non-constant, as it depends on the phenotypic variable.

The framework provides a multiscale representation of migration of cell populations, from single-cell dynamics to population-level behaviours. In fact, the starting point is a microscopic model that describes the dynamics of single cells in terms of stochastic processes. This model takes into account both cell movement, wherein cell reorientation is driven by non-local sensing of the surrounding environment, and environment-induced changes in the cytoskeletal structure. The mesoscopic counterpart of this model, which is formally derived by extending the limiting procedure employed in Conte and Loy (2024) to the phenotypically heterogeneous scenario considered here, is a phenotype-structured kinetic equation for the distribution of cells in each of phase and phenotype space. Then, considering an appropriately rescaled version of this kinetic equation, we formally derive the corresponding macroscopic model, which takes the form of a partial differential equation (PDE) for the cell number density. Being derived from an underlying model of single-cell dynamics, such a macroscopic equation constitutes a sound population-level model of cell migration.

To validate the formal procedures employed to derive the macroscopic model from the microscopic model, through the mesoscopic one, we first compare the results of numerical simulations of the two models. We then compare numerical solutions of the macroscopic model with the results of cell locomotion assays from Goodman et al. (1989), to test the ability of the model to recapitulate qualitative features of experimental observations. We finally report on numerical solutions of a reduced macroscopic model, and compare them with numerical simulations of the original model, in order to assess the robustness of the observed patterns of cell migration. The modelling framework presented here provides a mean to investigate the movement of phenotypically heterogeneous cell populations driven by environmental sensing, in scenarios where the extracellular environment may induce changes in direction, speed, and sensing radius of the cells. In particular, in this article we explore how the interplay between environment-induced changes in speed and in sensing radius may impact on the migration of cell populations.

## Outline of the article

In section 2, we formulate the microscopic model, the mesoscopic and macroscopic counterparts of which are then presented in section 3. In section 4, we report on Monte Carlo simulations of the microscopic model and numerical solutions of the PDE that defines the macroscopic model. In section 5, we conclude with a discussion and propose some future research directions.

## 2 A microscopic model for the migration of phenotypically heterogeneous cell populations

We consider a phenotypically heterogeneous cell population, where cells may differ in structural properties of their cytoskeleton which are implicated in cell migration. We focus on the case where the environment surrounding the cells can affect both reorientation processes, which underlie cell movement, and mechanisms of phenotypic variation, which drive changes in the cytoskeletal structure.
## 2.1 Mathematical modelling of the cell microscopic state

In order to provide a microscopic-scale description of the dynamics of the cell population, at time $t \in \mathbb{R}_0^+$, where $\mathbb{R}_0^+$ is the set of non-negative real numbers (i.e. $\mathbb{R}_0^+ := \mathbb{R}^+ \cup \{0\}$), we describe the microscopic state of the cells by the quadruplet $(\mathbf{X}_t, V_t, \hat{\mathbf{V}}_t, Y_t)$. The random vector $\mathbf{X}_t \in \mathbb{R}^d$, with $d = 1, 2, 3$ depending on the biological problem under study, represents the spatial position of a cell. Moreover, denoting by $\mathbf{V}_t = V_t \hat{\mathbf{V}}_t$ the random vector that represents the cell velocity, $\hat{\mathbf{V}}_t \in \mathbb{S}^{d-1}$ is the unit vector in the direction of $\mathbf{V}_t$ (i.e. the direction of cell polarity and thus of cell motion) and $V_t \in [0, V_{\max}]$ is the corresponding speed. Here, $\mathbb{S}^{d-1}$ is the unit sphere boundary in $\mathbb{R}^d$ and $V_{\max} \in \mathbb{R}^+$ is the maximal cell speed. In the remainder of the article, we will also use the compact notation $\mathcal{V} := [0, V_{\max}] \times \mathbb{S}^{d-1}$. Finally, the random variable $Y_t \in [0, 1]$ models the cell phenotypic state and captures inter-cellular variability in structural properties of the cytoskeleton. In particular, motivated by the experimental results presented in Goodman et al. (1989), we make the following assumptions:

**Assumption 1.** Smaller values of $Y_t$ correlate with a stiffer cytoskeleton (i.e. a cytoskeleton characterised by a higher level of organisation of actin and $\alpha$-actinin into stress fibres and vinculin into focal contacts), which leads to the formation of shorter cell surface extensions (i.e. lamellae).

**Assumption 2.** Larger values of $Y_t$ correlate with a more flexible cytoskeleton (i.e. a cytoskeleton characterised by more sparsely distributed actin, $\alpha$-actinin, and vinculin), which leads to the formation of longer cell surface extensions (i.e. pseudopods).

These assumptions imply that the movement of cells in phenotypic states close to 0 is driven by the extension of short lamellae, whereas the movement of cells in phenotypic states close to 1 is driven by the extension of long pseudopods (cf. the schematics in Figure 1).

<figure><img src="image_1.png" /><figcaption>Figure 1: Schematics illustrating the relationships between the phenotypic state, $y$, and cell characteristics.</figcaption></figure>

The microscopic state of the cells evolves in time due to cell movement, consisting of an alternation of runs over straight lines and reorientations, and phenotypic changes, which we assume to be independent processes. In particular, we model the evolution of the microscopic state $(\mathbf{X}_t, V_t, \hat{\mathbf{V}}_t, Y_t)$ between time $t$ and $t + \Delta t$, with $\Delta t \in \mathbb{R}^+$, through the following system (Conte
and Loy, 2024; Pareschi and Toscani, 2013)

$$
\left\{
\begin{aligned}
Y_{t+\Delta t} &= (1 - \Lambda)Y_t + \Lambda Y'_t \\
\hat{\mathbf{V}}_{t+\Delta t} &= (1 - M)\hat{\mathbf{V}}_t + M\hat{\mathbf{V}}'_t \\
V_{t+\Delta t} &= (1 - M)V_t + MV'_t \\
\mathbf{X}_{t+\Delta t} &= \mathbf{X}_t + \Delta t V_t \hat{\mathbf{V}}_t.
\end{aligned}
\right.
\qquad (1)
$$

In the system (1), $\Lambda$ and $M$ are independent Bernoulli random variables with parameters $\lambda\Delta t$ and $\mu\Delta t$, respectively, where, as per the modelling strategies introduced in the next subsections, the parameter $\lambda \in \mathbb{R}^+$ is the rate of phenotypic changes and the parameter $\mu \in \mathbb{R}^+$ is the rate at which cells change their velocity. Note that we are implicitly assuming $\Delta t$ to be small enough so that $\lambda\Delta t \le 1$ and $\mu\Delta t \le 1$. The equations in system (1) are such that:

* if $\Lambda = 1$ then a phenotypic change occurs and the cell transitions from the original phenotypic state $Y_t$ to a new phenotypic state represented by the random variable $Y'_t \in [0, 1]$, whereas if $\Lambda = 0$ then no phenotypic change occurs and the cell remains in the phenotypic state $Y_t$;

* similarly, if $M = 1$ then a velocity change occurs and the cell velocity switches from $(V_t, \hat{\mathbf{V}}_t)$ to $(V'_t, \hat{\mathbf{V}}'_t) \in \mathcal{V}$, whilst if $M = 0$ then no velocity change occurs and the cell keeps its velocity $(V_t, \mathbf{V}_t)$;

* the cell position evolves according to a free-particle drift.

We assume $Y'_t$, $\hat{\mathbf{V}}'_t$, and $V'_t$ to be distributed according to the following distribution functions

$$Y'_t \sim K[\mathcal{S}^{\ddagger}](y'|y), \quad \hat{\mathbf{V}}'_t \sim \mathcal{B}[\mathcal{S}](\hat{\mathbf{v}}'; y), \quad V'_t \sim \Psi[\mathcal{S}^{\dagger}](v'; \hat{\mathbf{v}}', y). \qquad (2)$$

As detailed in the next subsections, in (2):

* the vector functions

$$\mathcal{S}^{\ddagger} : (t, \mathbf{x}) \in \mathbb{R}_0^+ \times \mathbb{R}^d \mapsto \mathcal{S}^{\ddagger}(t, \mathbf{x}) \in \mathbb{R}^{m^{\ddagger}}, \quad \mathcal{S} : (t, \mathbf{x}) \in \mathbb{R}_0^+ \times \mathbb{R}^d \mapsto \mathcal{S}(t, \mathbf{x}) \in \mathbb{R}^m,$$

$$\mathcal{S}^{\dagger} : (t, \mathbf{x}) \in \mathbb{R}_0^+ \times \mathbb{R}^d \mapsto \mathcal{S}^{\dagger}(t, \mathbf{x}) \in \mathbb{R}^{m^{\dagger}},$$

which are assumed to be given, are introduced in order to incorporate the influence of the surrounding environment on cell dynamics into the model. Specifically: the vector function $\mathcal{S}^{\ddagger}(t, \mathbf{x})$ describes the spatial distributions at time $t$ of $m^{\ddagger}$ environmental factors that may affect the phenotypic state of the cells; the vector function $\mathcal{S}(t, \mathbf{x})$ describes the spatial distributions at time $t$ of $m$ environmental factors that may affect the direction of polarity of the cells; and the vector function $\mathcal{S}^{\dagger}(t, \mathbf{x})$ describes the spatial distributions at time $t$ of $m^{\dagger}$ environmental factors that may affect the cell speed. These could be, for instance, the concentrations of abiotic components of the cell microenvironment, such as nutrients, growth factors, adhesion sites on the ECM, and adhesive glycoproteins like laminin and fibronectin.

* The kernel $K[\mathcal{S}^{\ddagger}](y'|y)$ is a conditional probability distribution that describes the probability for a cell to transition into the phenotypic state $y'$ given that the cell is in the phenotypic state $y$.

* The probability distribution $\mathcal{B}[\mathcal{S}](\hat{\mathbf{v}}'; y)$ describes the probability for a cell in the phenotypic state $y$ to acquire the post-reorientation velocity $\hat{\mathbf{v}}'$.

* The probability distribution $\Psi[\mathcal{S}^{\dagger}](v'; \hat{\mathbf{v}}', y)$ describes the probability for a cell in the phenotypic state $y$ with post-reorientation velocity $\hat{\mathbf{v}}'$ to acquire the post-reorientation speed $v'$.

Note that the probability distributions $K[\mathcal{S}^{\ddagger}]$, $\mathcal{B}[\mathcal{S}]$, and $\Psi[\mathcal{S}^{\dagger}]$ depend on $(t, \mathbf{x})$ through the distributions of the environmental factors $\mathcal{S}^{\ddagger}(t, \mathbf{x})$, $\mathcal{S}(t, \mathbf{x})$, and $\mathcal{S}^{\dagger}(t, \mathbf{x})$.
## 2.2 Mathematical modelling of phenotypic changes

Due to phenotypic changes, which occur at rate $\lambda$, cells at position $\mathbf{x}$ in the phenotypic state $y$ can transition into a new phenotypic state $y'$ with a probability that is given by the kernel $K[\mathcal{S}^{\ddagger}](y'|y)$. The dependence of the phenotypic transition kernel $K$ on the distributions of environmental factors $\mathcal{S}^{\ddagger}(t, \mathbf{x})$ captures the fact that phenotypic changes undergone at time $t$ by cells at position $\mathbf{x}$ may be affected by the local cellular microenvironment (Gkretsi and Stylianopoulos, 2018; Goodman et al., 1989; Trepat et al., 2012). We let the kernel $K$ satisfy the following assumptions

$$K[\mathcal{S}^{\ddagger}](\cdot|y) > 0, K[\mathcal{S}^{\ddagger}](\cdot|y) \in L^1([0, 1]), \quad \int_0^1 K[\mathcal{S}^{\ddagger}](y'|y) dy' = 1 \text{ for a.e. } y \in [0, 1], \quad (3)$$

so that $K[\mathcal{S}^{\ddagger}](\cdot|y)$ is a probability distribution on $[0, 1]$. We then define its average as

$$y_K[\mathcal{S}^{\ddagger}](y) := \int_0^1 y' K[\mathcal{S}^{\ddagger}](y'|y) dy', \quad (4)$$

that is, $y_K[\mathcal{S}^{\ddagger}](y)$ is the mean phenotypic state of $K[\mathcal{S}^{\ddagger}](\cdot|y)$ – i.e. the phenotypic state into which, on average, cells in the phenotypic state $y$ enter as a result of phenotypic changes. This depends on $(t, \mathbf{x})$ through the distributions of the environmental factors $\mathcal{S}^{\ddagger}(t, \mathbf{x})$, and it is such that $y_K : (t, \mathbf{x}, y) \in \mathbb{R}_0^+ \times \mathbb{R}^d \times [0, 1] \to [0, 1]$. We also assume that

$$K[\mathcal{S}^{\ddagger}](y'|\cdot) \in L^{\infty}([0, 1]) \text{ for a.e. } y' \in [0, 1]. \quad (5)$$

## 2.3 Mathematical modelling of cell reorientation

We model cell reorientation as a velocity-jump process (Othmer et al., 1988; Stroock, 1974), whereby cells at position $\mathbf{x}$ in the phenotypic state $y$ may change their velocity at rate $\mu$ and acquire a new velocity (i.e. the post-reorientation velocity), which is prescribed by a turning kernel $T$. In analogy with previous works on the mathematical modelling of cell movement (Chauviere et al., 2007; Loy and Preziosi, 2020a,b), we assume the cell velocity acquired upon reorientation to be independent from the previous one. Moreover, we let the post-reorientation velocity be affected by cell sensing of the surrounding environment. This is mediated by membrane receptors located along surface extensions, which enable the cells to detect environmental factors within a finite sensing radius and then adjust the direction of their polarity (i.e. their direction of motion) and their speed accordingly. In particular, here we generalise the modelling strategies presented in (Loy and Preziosi, 2020a,b) so as to capture cell phenotypic variability by introducing a function

$$R : y \in [0, 1] \mapsto R(y) \in [R_{\min}, R_{\max}] \subset \mathbb{R}^{+}, \quad 0 < R_{\min} < R_{\max}$$

that represents the sensing radius of cells in the phenotypic state $y$. Since, under Assumptions 1 and 2, cells in phenotypic states modelled by larger (smaller) values of $y$ exhibit longer (shorter) surface extensions, and will thus sense environmental factors over a longer (shorter) distance, we assume the function $R$ to be such that

$$R(0) = R_{\min}, \quad R(1) = R_{\max}, \quad \frac{d}{dy} R(y) > 0 \quad \forall y \in [0, 1]. \quad (6)$$

In this framework, we let the turning kernel $T \equiv T[\mathcal{S}, \mathcal{S}^{\dagger}](\boldsymbol{\nu}'; y)$, which prescribes the post-reorientation velocity $\boldsymbol{\nu}' = (\nu', \hat{\nu}')$ acquired by cells at position $\mathbf{x}$ in the phenotypic state $y$ upon sensing of the distributions of environmental factors $\mathcal{S}$ and $\mathcal{S}^{\dagger}$ over a neighbourhood of $\mathbf{x}$ of maximal radius $R(y)$ (i.e. the sensing region), satisfy the following assumptions

$$T[\mathcal{S}, \mathcal{S}^{\dagger}](\cdot; y) > 0, T[\mathcal{S}, \mathcal{S}^{\dagger}](\cdot; y) \in L^1(\mathcal{V}), \quad \int_{\mathcal{V}} T[\mathcal{S}, \mathcal{S}^{\dagger}](\boldsymbol{\nu}'; y) d\boldsymbol{\nu}' = 1 \text{ for a.e. } y \in [0, 1], \quad (7)$$
so that $T[S, S^\dagger](\cdot; y)$ is a probability distribution on $\mathcal{V}$. Setting $\mathbf{v}' = v'\hat{\mathbf{v}}'$, we define

$$ \mathbf{u}_T[S, S^\dagger](y) := \int_{\mathcal{V}} \mathbf{v}' T[S, S^\dagger](\mathbf{v}'; y) \mathrm{d}\mathbf{v}', \quad (8) $$

that is, $\mathbf{u}_T[S, S^\dagger](y)$ is the mean velocity of $T[S, S^\dagger](\cdot; y)$ – i.e. the post-reorientation velocity which is acquired, on average, by cells in the phenotypic state $y$. This depends on $(t, \mathbf{x})$ through the distributions of the environmental factors $S(t, \mathbf{x})$ and $S^\dagger(t, \mathbf{x})$, and it is such that $\mathbf{u}_T: (t, \mathbf{x}, y) \in \mathbb{R}_0^+ \times \mathbb{R}^d \times [0, 1] \to \mathcal{V}$.

In particular, building on the modelling strategies presented in (Loy and Preziosi, 2020a,b), we use the definition

$$ T[S, S^\dagger](\mathbf{v}'; y) := \mathcal{B}[S](\hat{\mathbf{v}}'; y) \Psi[S^\dagger](\mathbf{v}'; \hat{\mathbf{v}}', y) \quad (9) $$

where:

* The probability distribution $\mathcal{B}[S](\hat{\mathbf{v}}'; y)$, which prescribes the direction of the post-reorientation velocity $\hat{\mathbf{v}}'$ of a cell at position $\mathbf{x}$ in the phenotypic state $y$, based on the distributions of environmental factors $S$ in a neighbourhood of $\mathbf{x}$ of radius $r \in [0, R(y)]$, is defined as

$$ \mathcal{B}[S](\hat{\mathbf{v}}'; y) := \frac{\int_{0}^{R(y)} \gamma_S(r, y) b[S](r, \hat{\mathbf{v}}') \mathrm{d}r}{\int_{S^{d-1}}^{R(y)} \gamma_S(r, y) b[S](r, \hat{\mathbf{v}}') \mathrm{d}r \, \mathrm{d}\hat{\mathbf{v}}'} \quad (10) $$

Details on $\gamma_S(r, y)$ and $b[S](r, \hat{\mathbf{v}}')$ are provided below.

* The probability distribution $\Psi[S^\dagger](\mathbf{v}'; \hat{\mathbf{v}}', y)$, which prescribes the post-reorientation speed $\mathbf{v}'$ of a cell at position $\mathbf{x}$ in the phenotypic state $y$, based on the distributions of environmental factors $S^\dagger$ in a neighbourhood of $\mathbf{x}$ of radius $r^\dagger \in [0, R(y)]$ along the post-reorientation direction of motion $\hat{\mathbf{v}}'$, is defined as

$$ \Psi[S^\dagger](\mathbf{v}'; \hat{\mathbf{v}}', y) := \frac{\int_{0}^{R(y)} \gamma_{S^\dagger}(r^\dagger, y) \psi[S^\dagger(t, \mathbf{x} + r^\dagger \hat{\mathbf{v}}')] (\mathbf{v}') \mathrm{d}r^\dagger}{\int_{0}^{V_{\max}} \int_{0}^{R(y)} \gamma_{S^\dagger}(r^\dagger, y) \psi[S^\dagger(t, \mathbf{x} + r^\dagger \hat{\mathbf{v}}')] (\mathbf{v}') \mathrm{d}r^\dagger \mathrm{d}\mathbf{v}'} \quad (11) $$

Details on $\gamma_{S^\dagger}(r^\dagger, y)$ and $\psi[S^\dagger](\mathbf{v}')$ are provided below.

In the definitions given by (10) and (11):

* $\gamma_S : \mathbb{R}_0^+ \times [0, 1] \to \mathbb{R}_0^+$ and $\gamma_{S^\dagger} : \mathbb{R}_0^+ \times [0, 1] \to \mathbb{R}_0^+$ are weights modelling how sensing of the distributions of environmental factors $S$ and $S^\dagger$ over neighbourhoods of $\mathbf{x}$ of radii $r \in [0, R(y)]$ and $r^\dagger \in [0, R(y)]$ affects the direction of polarity and the speed of cells in the phenotypic state $y$. For instance, if the direction of polarity and the speed of cells in the phenotypic state $y$ are affected only by the environment sensed at the edge of the sensing region, then the functions $\gamma_S$ and $\gamma_{S^\dagger}$ can be defined as

$$ \gamma_S(r, y) := \delta(r - R(y)) \quad \text{and} \quad \gamma_{S^\dagger}(r^\dagger, y) := \delta(r^\dagger - R(y)), \quad (12) $$

where $\delta(\cdot)$ is the Dirac delta. Alternatively, these functions may be defined: as Heaviside step functions of $r$ and $r^\dagger$, when the direction of polarity and the speed of the cells are affected by the environmental conditions within the sensing region in a uniform manner; or as monotonically decreasing functions of $r$ and $r^\dagger$, so as to account for a reduced influence of environmental cues closer to the edge of the sensing region (Loy and Preziosi, 2020a,b).

* $b[S](r, \hat{\mathbf{v}}')$ is a non-negative function that describes how reorientation of cell polarity, driven by the distributions of environmental factors $S$ in a neighbourhood of $\mathbf{x}$ of radius $r \in [0, R(y)]$, leads the cells to move in direction $\hat{\mathbf{v}}'$. For additional details about this function, we refer the reader to (Loy and Preziosi, 2020a).
* $\psi[\mathcal{S}^\dagger(t, \mathbf{x}+r^\dagger\hat{\mathbf{v}}')](v')$ is a probability distribution that prescribes the post-reorientation speed $v'$, based on the distribution of the environmental factors $\mathcal{S}^\dagger$ in a neighbourhood of $\mathbf{x}$ of radius $r^\dagger \in [0, R(y)]$ along the post-reorientation direction of motion $\hat{\mathbf{v}}'$, which is such that

$$ \int_{0}^{V_{\max}} \psi[\mathcal{S}^{\dagger}](v') \, dv' = 1, \quad \int_{0}^{V_{\max}} v' \psi[\mathcal{S}^{\dagger}](v') \, dv' = u_{\psi}[\mathcal{S}^{\dagger}]. \qquad (13) $$

In (13), $u_\psi[\mathcal{S}^\dagger]$ is the mean speed of $\psi[\mathcal{S}^\dagger](v')$, that is, the post-reorientation speed which is acquired, on average, by the cells. This depends on $(t, \mathbf{x})$ through the distributions of the environmental factors $\mathcal{S}^\dagger(t, \mathbf{x})$, and it is such that $u_\psi : (t, \mathbf{x}) \in \mathbb{R}_0^+ \times \mathbb{R}^d \to [0, V_{\max}]$.

Overall, the modelling choice corresponding to (9)-(11) implies that at each reorientation there is a simultaneous change of direction and speed: the polarisation direction changes into a new one $\hat{\mathbf{v}}'$ and the new speed $v'$ is chosen by sensing the external factor $\mathcal{S}^\dagger$ along $\hat{\mathbf{v}}'$. Note that here we assume $b[\mathcal{S}]$ and $\psi[\mathcal{S}^\dagger]$ to be such

$$ \mathcal{B}[\mathcal{S}](\cdot; y) \in L^1(\mathbb{S}^{d-1}), \quad \Psi[\mathcal{S}^\dagger](\cdot; \hat{\mathbf{v}}', y) \in L^1([0, V_{\max}]), $$

and we also assume $\gamma_S$ and $\gamma_{S^\dagger}$ to be such that $\mathcal{B}[\mathcal{S}]$ and $\Psi[\mathcal{S}^\dagger]$ are measurable functions of $y$. Under these regularity assumptions, $T[\mathcal{S}, \mathcal{S}^\dagger]$ defined via (9)-(11) is a measurable function of $y$. Moreover, using the Fubini-Tonelli theorem, one can see that $T[\mathcal{S}, \mathcal{S}^\dagger]$ does satisfy the normalisation assumption in (7).

# 3 Corresponding mesoscopic and macroscopic models

In this section, we formally derive first the mesoscopic counterpart of the microscopic model presented in the previous section and then the corresponding macroscopic model.

## 3.1 Preliminaries and notation

We let

$$ (\mathbf{X}_t, V_t, \hat{\mathbf{V}}_t, Y_t) \sim f(t, \mathbf{x}, v, \hat{\mathbf{v}}, y), \qquad (14) $$

where the distribution function $f: \mathbb{R}_0^+ \times \mathbb{R}^d \times [0, V_{\max}] \times \mathbb{S}^{d-1} \times [0, 1] \to \mathbb{R}_0^+$, which represents the cell distribution in each of phase and phenotype space, is a probability density function and is thus such that

$$ \int_{\mathbb{R}^d} \int_0^{V_{\max}} \int_{\mathbb{S}^{d-1}} \int_0^1 f(t, \mathbf{x}, v, \hat{\mathbf{v}}, y) \, dy \, d\hat{\mathbf{v}} \, dv \, d\mathbf{x} = 1. \qquad (15) $$

In the remainder of the article, we will use the more compact notation

$$ f(t, \mathbf{x}, v, \hat{\mathbf{v}}, y) \equiv f(t, \mathbf{x}, \boldsymbol{\nu}, y), \quad \boldsymbol{\nu} \equiv (v, \hat{\mathbf{v}}) \in \mathcal{V}. $$

The macroscopic quantities are defined as the statistical moments of $f$ at $(t, \mathbf{x})$ fixed:

* the number density (per unit volume) of the cell population (i.e. the cell density)

$$ \rho(t, \mathbf{x}) := \int_{\mathcal{V}} \int_{0}^{1} f(t, \mathbf{x}, \boldsymbol{\nu}, y) \, dy \, d\boldsymbol{\nu}; \qquad (16) $$

* the normalised distribution of the cell population in the phase space

$$ p(t, \mathbf{x}, \boldsymbol{\nu}) := \frac{1}{\rho(t, \mathbf{x})} \int_{0}^{1} f(t, \mathbf{x}, \boldsymbol{\nu}, y) \, dy; \qquad (17) $$
• the mean velocity of the cell population (i.e. the ensemble average velocity)

$$ \mathbf{U}(t, \mathbf{x}) := \int_{\mathcal{V}} \mathbf{v} p(t, \mathbf{x}, \boldsymbol{\nu}) \, d\boldsymbol{\nu} \qquad (18) $$

and the corresponding momentum

$$ \rho(t, \mathbf{x}) \mathbf{U}(t, \mathbf{x}) = \rho(t, \mathbf{x}) \int_{\mathcal{V}} \mathbf{v} p(t, \mathbf{x}, \boldsymbol{\nu}) \, d\boldsymbol{\nu}; \qquad (19) $$

• the second order tensor of the cell population

$$ \rho(t, \mathbf{x}) \mathcal{D}(t, \mathbf{x}) = \rho(t, \mathbf{x}) \int_{\mathcal{V}} (\mathbf{v} - \mathbf{U}) \otimes (\mathbf{v} - \mathbf{U}) p(t, \mathbf{x}, \boldsymbol{\nu}) \, d\boldsymbol{\nu}, \qquad (20) $$

being $\mathcal{D}(t, \mathbf{x})$ the variance-covariance matrix of $p$ for each $(t, \mathbf{x})$ fixed.

Moreover, we will use the following definitions of the counterparts of the above quantities for cells in a given phenotypic state $y$:

• the normalised distribution of the cell population in the phenotype space

$$ n(t, \mathbf{x}, y) := \frac{1}{\rho(t, \mathbf{x})} \int_{\mathcal{V}} f(t, \mathbf{x}, \boldsymbol{\nu}, y) \, d\boldsymbol{\nu}, \qquad (21) $$

whereby $\rho(t, \mathbf{x})n(t, \mathbf{x}, y)$ is the number density of cells in the phenotypic state $y$;

• the mean velocity of cells in the phenotypic state $y$

$$ \mathbf{u}(t, \mathbf{x}, y) := \frac{1}{\rho(t, \mathbf{x})n(t, \mathbf{x}, y)} \int_{\mathcal{V}} \mathbf{v} f(t, \mathbf{x}, \boldsymbol{\nu}, y) \, d\boldsymbol{\nu} \qquad (22) $$

and the corresponding momentum

$$ \rho(t, \mathbf{x})n(t, \mathbf{x}, y) \mathbf{u}(t, \mathbf{x}, y) = \int_{\mathcal{V}} \mathbf{v} f(t, \mathbf{x}, \boldsymbol{\nu}, y) \, d\boldsymbol{\nu}; \qquad (23) $$

• the second order tensor of cells in the phenotypic state $y$

$$ \rho(t, \mathbf{x})n(t, \mathbf{x}, y) \mathcal{D}(t, \mathbf{x}, y) = \int_{\mathcal{V}} (\mathbf{v} - \mathbf{u}) \otimes (\mathbf{v} - \mathbf{u}) f(t, \mathbf{x}, \boldsymbol{\nu}, y) \, d\boldsymbol{\nu}, \qquad (24) $$

where $n(t, \mathbf{x}, y)\mathcal{D}(t, \mathbf{x}, y)$ is the variance-covariance matrix of $f$ with respect to $\boldsymbol{\nu}$ for each $(t, \mathbf{x}, y)$ fixed.

Finally, we will use the following definition of the mean phenotypic state of the cell population

$$ \bar{y}(t, \mathbf{x}) := \int_{0}^{1} y n(t, \mathbf{x}, y) \, dy, \qquad (25) $$

which corresponds to the average phenotype of the normalised distribution $n$.

## 3.2 The corresponding mesoscopic and macroscopic models

The mesoscopic model corresponding to the microscopic model presented in section 2 consists of a phenotype-structured kinetic equation for the distribution function $f(t, \mathbf{x}, \boldsymbol{\nu}, y)$ given in (14), the strong form of which reads as

$$ \left\{ \begin{array}{l} \partial_t f + \mathbf{v} \cdot \nabla_{\mathbf{x}} f = \mu (\rho T[\mathcal{S}, \mathcal{S}^\dagger] n - f) \\ \qquad + \lambda \left( \int_0^1 K[\mathcal{S}^\ddagger](y|y') f(t, \mathbf{x}, \boldsymbol{\nu}, y') dy' - f \right), \\ \\ n(t, \mathbf{x}, y) := \frac{1}{\rho(t, \mathbf{x})} \int_{\mathcal{V}} f(t, \mathbf{x}, \boldsymbol{\nu}, y) \, d\boldsymbol{\nu}. \end{array} \right. \qquad (26) $$
This kinetic model is formally derived in Section A in the supplementary material, starting from system (1) complemented with relations (2), using a limiting procedure which generalises the one employed, for instance, in (Conte and Loy, 2024) to the case of phenotypically heterogeneous cell populations.

**Remark 1.** In (26), the term

$$ \int_{0}^{1} K[\mathcal{S}^{\dagger}](y|y') f(t, \mathbf{x}, \mathbf{v}, y') dy' - f $$

could also be rewritten as a differential term by using the quasi-invariant limit approach for transition probabilities, as similarly done in (Loy and Tosin, 2020).

Our aim is now to derive a macroscopic model, that is, a set of PDEs describing the evolution of macroscopic quantities which are conserved by the cellular processes (i.e. reorientations and phenotypic changes) incorporated into the mesoscopic model. In order to obtain a closed macroscopic model, different approaches could be adopted in different regimes (see Section B in the supplementary material). Amongst these, in the light of the experimental results presented in (Goodman et al., 1989) and the considerations made in Section C in the supplementary material, introducing a small parameter $\epsilon \in \mathbb{R}^{+}$, one can consider the following scaling

$$ \lambda = \frac{1}{\epsilon}, \quad \mu = \frac{1}{\epsilon}, \qquad (27) $$

which corresponds to scenarios where phenotypic changes and cell reorientation occur on similar time scales, which are faster than the time scale of spatial dynamics of the cells. The non-dimensionalisation leading to (27) is provided in Section C in the supplementary material.

To this end, denoting the solution to the phenotype-structured kinetic equation (26) under the parameter scaling (27) by $f_{\epsilon}$, we make the following expansion ansatz (i.e. the Chapman-Enskog expansion)

$$ f_{\epsilon} = f_0 + \epsilon f^{\perp}, \quad \rho_0 := \int_{\mathcal{V}} \int_{0}^{1} f_0 dy d\nu = \rho, \quad \rho^{\perp} := \int_{\mathcal{V}} \int_{0}^{1} f^{\perp} dy d\nu = 0. \qquad (28) $$

The ansatz (28) implies that the mass is concentrated in the leading order term, $f_0$, while the correction, $f^{\perp}$, carries no mass, and thus

$$ \rho_{\epsilon} = \rho_{0} = \rho. \qquad (29) $$

From (28), one obtains analogous expansions for the marginals. In fact, integrating $f_{\epsilon}/\rho$ over $[0, 1]$ and using (28) gives the expansion for the marginal $p_{\epsilon}$ along with the mass of the corresponding leading order term and of the correction, i.e.

$$ p_{\epsilon} = p_{0} + \epsilon p^{\perp}, \quad p_{0} = \frac{1}{\rho} \int_{0}^{1} f_{0} dy, \quad p^{\perp} = \frac{1}{\rho} \int_{0}^{1} f^{\perp} dy, \quad \int_{\mathcal{V}} p_{0} d\nu = 1, \quad \int_{\mathcal{V}} p^{\perp} d\nu = 0. \qquad (30) $$

Similarly, integrating $f_{\epsilon}/\rho$ over $\mathcal{V}$ and using (28) gives the expansion for the marginal $n_{\epsilon}$ along with the mass of the corresponding leading order term and of the correction, i.e.

$$ n_{\epsilon} = n_{0} + \epsilon n^{\perp}, \quad n_{0} = \frac{1}{\rho} \int_{\mathcal{V}} f d\nu, \quad n^{\perp} = \frac{1}{\rho} \int_{\mathcal{V}} n^{\perp} d\nu, \quad \int_{0}^{1} n_{0} dy = 1, \quad \int_{0}^{1} n^{\perp} dy = 0. \qquad (31) $$

**Remark 2.** We now introduce the following integral operator defined on the Banach space $L^{1}([0, 1])$ into itself

$$ \mathcal{K}: L^{1}([0, 1]) \to L^{1}([0, 1]) $$

$$ \mathcal{K}[g] = \int_{0}^{1} K[\mathcal{S}^{\dagger}](y|y') g(t, \mathbf{x}, y') dy', $$
which is linear, bounded (cf. assumption (5)), and compact. Because of the positivity assumption in (3), this operator is strictly positive on the cone of positive functions of $L^1([0, 1])$, which is reproducing and also total. Therefore, the Krein-Rutman theorem (see e.g. (Deimling, 1985, Theorem 19.2)) allows one to conclude that this integral operator admits a spectral radius associated with a positive eigenfunction, $g_1(t, \mathbf{x}, y)$, which depends on $(t, \mathbf{x})$ through $S^\dagger(t, \mathbf{x})$ (i.e. $g_1(t, \mathbf{x}, y) \equiv g_1[S^\dagger](y)$) and it is unique up to normalisation, i.e.

$$ \int_{0}^{1} g_{1}[S^{\dagger}](y) \, dy = 1. $$

In particular, note that, due to the normalisation assumption in (3), the eigenfunction $g_1(t, \mathbf{x}, y)$ is associated with the eigenvalue 1. Moreover, if the phenotypic transition kernel $K[S^\dagger]$ also satisfies the following assumption

$$ K[S^\dagger](y|y') \equiv K[S^\dagger](y), \qquad (32) $$

that is, if the new phenotypic states acquired by the cells upon phenotypic changes are independent of the original ones, then

$$ K[g] = \int_{0}^{1} K[S^{\dagger}](y) g(t, \mathbf{x}, y') \, dy' = K[S^{\dagger}](y) \int_{0}^{1} g(t, \mathbf{x}, y') \, dy'. $$

In this case, since $g_1$ is the eigenfunction of $K[g]$ associated with the eigenvalue 1 and it is normalised, then

$$ g_1[S^\dagger](y) = K[S^\dagger](y). \qquad (33) $$

In the remainder of this section, we will be letting assumption (32) hold, so that the result (33) holds as well. Moreover, we will be assuming the distributions of environmental factors $S$, $S^\dagger$, and $S^\ddagger$ to be constant in time, i.e.

$$ S \equiv S(\mathbf{x}), \quad S^\dagger \equiv S^\dagger(\mathbf{x}), \quad S^\ddagger \equiv S^\ddagger(\mathbf{x}), \qquad (34) $$

in order to consider the scenario where the extracellular environment is stationary. Note that, in this scenario, the result (33) implies that $g_1[S^\dagger](y)$ is constant in time.

In the framework of assumptions (32) and (34), exploiting the result (33), we now employ a limiting procedure that generalises the one used, for instance, in (Hillen, 2006), to the case where a continuous phenotype structure and phenotypic changes are incorporated into the kinetic equation. As similarly done in (Hillen, 2006), we define the operator

$$ \mathcal{L}^{\nu}: L^{2}(\mathcal{V} \times [0, 1]) \to L^{2}(\mathcal{V} \times [0, 1]) \\ \qquad \varphi \mapsto T[S, S^{\dagger}](\nu; y) \int_{\mathcal{V}} \varphi(\nu, y) \, d\nu - \varphi, $$

so that $\mathcal{L}^\nu(f) = \rho T[S, S^\dagger](\nu; y)n - f$. Moreover, since, compared to (Hillen, 2006), in our model the microscopic state of the cells comprises both a velocity component and a phenotype component, the dynamics of which are governed by two independent processes, we introduce the additional operator

$$ \mathcal{L}^{y}: L^{2}(\mathcal{V} \times [0, 1]) \to L^{2}(\mathcal{V} \times [0, 1]) \\ \qquad \varphi \mapsto K[S^{\dagger}](y) \int_{0}^{1} \varphi(\nu, y) \, dy - \varphi, $$

so that $\mathcal{L}^y(f) = \rho K[S^\dagger](y)p - f$. Note that, on the same function space $L^2(\mathcal{V} \times [0, 1])$, we can also define the sum of the two operators $\mathcal{L} := \mathcal{L}^\nu + \mathcal{L}^y$, with

$$ \mathcal{L}: L^2(\mathcal{V} \times [0, 1]) \to L^2(\mathcal{V} \times [0, 1]) \\ \qquad \varphi \mapsto T[S, S^\dagger](\nu; y) \int_{\mathcal{V}} \varphi(\nu', y) \, d\nu' + K[S^\dagger](y) \int_{0}^{1} \varphi(\nu, y) \, dy - 2\varphi, \qquad (35) $$
which allows us to rewrite the phenotype-structured kinetic equation (26) under the parameter scaling (27) as

$$ \partial_t f_{\varepsilon} + \mathbf{v} \cdot \nabla_{\mathbf{x}} f_{\varepsilon} = \frac{1}{\varepsilon} \mathcal{L}(f_{\varepsilon}). \quad (36) $$

Substituting (28) into (36), one finds a hierarchy of equations in $\varepsilon$. Specifically, in $\varepsilon^0$

$$ \mathcal{L}(f_0) = 0, \quad (37) $$

while in $\varepsilon^1$

$$ \partial_t f_0 + \mathbf{v} \cdot \nabla_{\mathbf{x}} f_0 = \mathcal{L}(f^{\perp}). \quad (38) $$

The only conserved quantity is the mass (i.e. the number density at $(t, \mathbf{x})$ fixed), as

$\int_V \int_0^1 \mathcal{L}(f) \, dy \, d\nu = 0$, while no higher-order moments of $f$ are conserved. As a consequence, in the regime (27), wherein the operators $\mathcal{L}^\nu$ and $\mathcal{L}^y$ defining the operator $\mathcal{L}$ are of the same order, the macroscopic model will be a closed PDE for the cell density, that is, the number density of the distribution $f_0$ which makes the operator $\mathcal{L}$ vanish.

From (37) one obtains

$$ f_0(t, \mathbf{x}, \nu, y) = \frac{\rho(t, \mathbf{x})}{2} \left( T[\mathcal{S}, \mathcal{S}^{\dagger}](\nu; y) n^0(\mathbf{x}, y) + K[\mathcal{S}^{\ddagger}](y) p^0(\mathbf{x}, \nu) \right), $$

which requires the marginal stationary equilibria $n^0(\mathbf{x}, y)$ and $p^0(\mathbf{x}, \nu)$ that nullify, respectively, the operator corresponding to phenotypic changes and the operator corresponding to cell reorientation.

Integrating (36) with respect to $y$ one recovers the evolution equation for the marginal $p_\varepsilon$, that is,

$$ \partial_t (\rho p_\varepsilon) + \mathbf{v} \cdot \nabla_{\mathbf{x}} (\rho p_\varepsilon) = \frac{1}{\varepsilon} \bar{\mathcal{L}}^\nu(f_\varepsilon), \quad (39) $$

where

$$ \bar{\mathcal{L}}^\nu : L^2(\mathcal{V} \times [0, 1]) \to L^2(\mathcal{V}) $$

$$ \varphi \mapsto \int_0^1 T[\mathcal{S}, \mathcal{S}^{\dagger}](\nu; y) \int_{\mathcal{V}} \varphi(\nu', y) \, d\nu' \, dy - \int_0^1 \varphi(\nu, y) \, dy, \quad (40) $$

while the evolution equation for $n_\varepsilon$ is

$$ \partial_t (\rho n_\varepsilon) + \nabla_{\mathbf{x}} \cdot (\rho n_\varepsilon \mathbf{u}) = \frac{1}{\varepsilon} \bar{\mathcal{L}}^y(f_\varepsilon), \quad (41) $$

where

$$ \bar{\mathcal{L}}^y : L^2(\mathcal{V} \times [0, 1]) \to L^2([0, 1]) $$

$$ \varphi \mapsto K[\mathcal{S}^{\ddagger}](y) \int_{\mathcal{V}} \int_0^1 \varphi(\nu, y) \, dy \, d\nu - \int_{\mathcal{V}} \varphi(\nu, y) \, d\nu. \quad (42) $$

At leading order in $\varepsilon$, from the equations for the marginals (39) and (41) complemented with (30) and (31), one finds

$$ p_0(\mathbf{x}, \nu) = \int_0^1 T[\mathcal{S}, \mathcal{S}^{\dagger}](\nu; y) n_0(\mathbf{x}, y) \, dy, \quad n_0(\mathbf{x}, y) = K[\mathcal{S}^{\ddagger}](y), \quad (43) $$

i.e.

$$ p_0(\mathbf{x}, \nu) = T_K[\mathcal{S}, \mathcal{S}^{\dagger}, \mathcal{S}^{\ddagger}](\nu), \quad T_K[\mathcal{S}, \mathcal{S}^{\dagger}, \mathcal{S}^{\ddagger}](\nu) := \int_0^1 T[\mathcal{S}, \mathcal{S}^{\dagger}](\nu; y) K[\mathcal{S}^{\ddagger}](y) \, dy. \quad (44) $$

Therefore, the distribution $f_0$, i.e. the solution of (37), is

$$ f_0(t, \mathbf{x}, \nu, y) = \rho(t, \mathbf{x}) \mathcal{T}[\mathcal{S}, \mathcal{S}^{\dagger}, \mathcal{S}^{\ddagger}](\nu, y), \quad (45) $$

where

$$ \mathcal{T}[\mathcal{S}, \mathcal{S}^{\dagger}, \mathcal{S}^{\ddagger}](\nu, y) := \frac{K[\mathcal{S}^{\ddagger}](y)}{2} [T_K[\mathcal{S}, \mathcal{S}^{\dagger}, \mathcal{S}^{\ddagger}](\nu) + T[\mathcal{S}, \mathcal{S}^{\dagger}](\nu; y)]. \quad (46) $$
Note that from (46) one sees that the stationarity assumptions (34) are needed so as to ensure that $\mathcal{T}[S, S^\dagger, S^{\ddagger}]$ does not depend on the time variable. Note also that $\mathcal{T}$ is a probability distribution on $V \times [0, 1]$, since

$$ \int_V \int_0^1 \mathcal{T}[S, S^\dagger, S^{\ddagger}](\nu, y) \, dy \, d\nu = 1, \quad \forall \mathbf{x} \in \mathbb{R}^d, $$

due to (3), (7), and to the measurability of $\mathcal{T}$ as a function of $y$, which provides the regularity needed to apply the Fubini-Tonelli theorem and, in so doing, obtain the above normalisation condition. In (45) we have used (29). We also stress the fact that the time dependence of the distribution $f_0$ given by (45) is due exclusively to the time dependence of the number density $\rho$, while $p_0$ and $n_0$, which nullify, respectively, the operator $\bar{\mathcal{L}}^\nu$ and the operator $\bar{\mathcal{L}}^y$, do not depend on the time variable, as a result of the stationarity assumptions (34). Substituting (45) and (46) into (38) and integrating over $V \times [0, 1]$ gives, at leading order in $\epsilon$, the PDE

$$ \partial_t \rho + \nabla_{\mathbf{x}} \cdot (\rho \mathbf{U}_{\mathcal{T}}) = 0, \quad (47) $$

where $\mathbf{U}_{\mathcal{T}}$ is the average velocity of $\mathcal{T}$ integrated over $[0, 1]$, which is defined as

$$ \mathbf{U}_{\mathcal{T}} = \mathbf{U}_{\mathcal{T}}[S, S^{\dagger}, S^{\ddagger}] := \int_{0}^{1} \int_{V} \mathbf{v} \mathcal{T}[S, S^{\dagger}, S^{\ddagger}](\nu, y) \, d\nu \, dy. \quad (48) $$

Exploiting the regularity of the phenotypic transition kernel $K$ and the turning kernel $T$, it is possible to invert the order of integration, and, in so doing, obtain

$$ \mathbf{U}_{\mathcal{T}} = \int_{0}^{1} \mathbf{u}_{\mathcal{T}}[S, S^{\dagger}](y) K[S^{\ddagger}](y) \, dy, \quad (49) $$

with $\mathbf{u}_{\mathcal{T}}$ being defined via (8).

We now seek the first order correction in $\epsilon$ for (47). To this purpose, we need to determine the correction $f^\perp$ to the distribution $f_0$, as $f^\perp$ carries the information about the time evolving mass in the transient. To do so, we will need to invert the operator $\mathcal{L}$ and to define a Fredholm alternative. On $L^2(V \times [0, 1])$ we define the inner product

$$ h, g \in L^2(V \times [0, 1]) \quad \langle h, g \rangle := \int_V \int_0^1 h(\nu, y)g(\nu, y) \mathcal{T}^{-1}(\nu; y) \, dy \, d\nu, \quad (50) $$

so that $L^2(V \times [0, 1])$ endowed with this metric is a Hilbert space. As the kernel of $\mathcal{L}$ is

$$ \ker(\mathcal{L}) = \operatorname{span}\{\mathcal{T}\}, $$

the operator $\mathcal{L}$ will be inverted to the orthogonal to its kernel $\operatorname{span}\{\mathcal{T}\}^\perp$, that is, the subspace of functions which are orthogonal to $\mathcal{T}$ according to the scalar product defined via (50). Such an orthogonality condition implies that the functions in this subspace have zero number density. Then, the pseudo-inverse is defined as $\mathcal{L}^{-1} : \operatorname{span}\{\mathcal{T}\}^\perp \to \operatorname{span}\{\mathcal{T}\}^\perp$. As $\mathcal{L}$ is defined via (35), then $\mathcal{L}(\eta) = -2\eta$ for $\eta \in \operatorname{span}\{\mathcal{T}\}^\perp$, so that $\mathcal{L}^{-1}(\eta) = -\frac{1}{2}\eta$. Similar considerations hold for the operators $\bar{\mathcal{L}}^\nu$ and $\bar{\mathcal{L}}^y$ defined via (40) and (42). Specifically, we define the inner products

$$ h, g \in L^2(V) \quad \langle h, g \rangle := \int_V h(\nu)g(\nu)T_K^{-1}(\nu) \, d\nu $$

and

$$ h, g \in L^2([0, 1]) \quad \langle h, g \rangle := \int_0^1 h(y)g(y)K^{-1}(y) \, dy. $$

Since

$$ \ker(\bar{\mathcal{L}}^\nu) = \left\{ \varphi \in L^2(V \times [0, 1]) : \int_0^1 \varphi(\nu, y) \, dy = \int_{V \times [0, 1]} T(\nu; y) \varphi(\nu', y) \, dy \, d\nu' \right\} $$
and

$$ \ker(\bar{\mathcal{L}}^y) = \left\{ \varphi \in L^2(\mathcal{V} \times [0, 1]) : \int_{\mathcal{V}} \varphi(\boldsymbol{\nu}, y) \mathrm{d}\boldsymbol{\nu} = K(y) \int_{\mathcal{V} \times [0, 1]} \varphi \mathrm{d}y \mathrm{d}\boldsymbol{\nu} \right\}, $$

we have that

$$ \ker(\mathcal{L}) \subset \ker(\bar{\mathcal{L}}^{\nu}) \cap \ker(\bar{\mathcal{L}}^{y}). $$

Therefore, the pseudo-inverse functionals are

$$ \bar{\mathcal{L}}^{\nu^{-1}} : \operatorname{span}\{T_K\}^{\perp} \to \operatorname{span}\{\mathcal{T}\}^{\perp} \\ \eta \mapsto -\varphi \quad \text{with } \varphi \text{ s.t.} \quad \int_{0}^{1} \varphi(\boldsymbol{\nu}, y) \mathrm{d}\boldsymbol{\nu} = \eta(\boldsymbol{\nu}) $$

and

$$ \bar{\mathcal{L}}^{y^{-1}} : \operatorname{span}\{K\}^{\perp} \to \operatorname{span}\{\mathcal{T}\}^{\perp} \\ \eta \mapsto -\varphi \quad \text{with } \varphi \text{ s.t.} \quad \int_{\mathcal{V}} \varphi(\boldsymbol{\nu}, y) \mathrm{d}\boldsymbol{\nu} = \eta(y). $$

Considering (36) alongside (28) and including only the terms up to order 1 in $\varepsilon$ yields

$$ \partial_t f_0 + \mathbf{v} \cdot \nabla_{\mathbf{x}} (f_0 + \varepsilon f^{\perp}) = \frac{1}{\varepsilon} \mathcal{L}(f_0) + \mathcal{L}(f^{\perp}), \quad (51) $$

where the first term on the right-hand side vanishes uniformly in $\varepsilon$ due to (37). Since $\rho^{\perp} = 0$ because of (28), so that $f^{\perp} \in \operatorname{span}\{\mathcal{T}\}^{\perp}$, and thanks to (47), one can apply the pseudo-inverse $\mathcal{L}^{-1}$ to (51). In so doing, from (38) one obtains

$$ f^{\perp} = -\frac{1}{2}(\partial_t f_0 + \mathbf{v} \cdot \nabla_{\mathbf{x}} f_0) + \frac{1}{2}(\rho p^{\perp} K + \rho n^{\perp} T), $$

from which, using (45) along with the fact that $\mathcal{T}$ does not depend on the time variable, which allows one to state that $\partial_t f_0 = \mathcal{T}\partial_t \rho$, and the PDE (47) for expressing $\partial_t \rho = -\nabla_{\mathbf{x}} \cdot (\mathbf{U}_{\mathcal{T}}\rho)$, one finds

$$ f^{\perp} = \frac{1}{2}(\nabla_{\mathbf{x}} \cdot (\mathbf{U}_{\mathcal{T}}\rho)\mathcal{T} - \mathbf{v} \cdot \nabla_{\mathbf{x}}(\rho\mathcal{T})) + \frac{1}{2}(\rho p^{\perp} K + \rho n^{\perp} T). \quad (52) $$

An analogous inversion procedure can be applied to both (39) and (41) so as to determine $p^{\perp}$ and $n^{\perp}$. One has

$$ \partial_t(\rho p_0) + \mathbf{v} \cdot \nabla_{\mathbf{x}} \rho p_0 = \bar{\mathcal{L}}^{\nu}(f^{\perp}), \quad (53) $$

and

$$ \partial_t(\rho n_0) + \nabla_{\mathbf{x}} \cdot \left[ \rho n_0 \frac{1}{2} (\mathbf{u}_T + \mathbf{U}_{\mathcal{T}}) \right] = \bar{\mathcal{L}}^y(f^{\perp}), \quad (54) $$

where $\bar{\mathcal{L}}^{\nu}(f^{\perp}) = -\rho p^{\perp}$, and $\bar{\mathcal{L}}^{y}(f^{\perp}) = -\rho n^{\perp}$ as $\rho^{\perp} = 0$ and $f^{\perp} \in \operatorname{span}\{\mathcal{T}\}^{\perp}$. Using the last one of relations (30), which ensures that $p^{\perp} \in \operatorname{span}\{T_K\}^{\perp}$, the definition (43)-(44) of $p_0$, and the fact that integrating the left-hand side of (53) over $\mathcal{V}$ gives (47), one can apply the pseudo-inverse $\bar{\mathcal{L}}^{y^{-1}}$ and obtain

$$ \rho p^{\perp} = \nabla_{\mathbf{x}} \cdot (\mathbf{U}_{\mathcal{T}} \rho) T_K - \mathbf{v} \cdot \nabla_{\mathbf{x}} (\rho T_K). $$

Using the last one of relations (31), which ensures that $n^{\perp} \in \operatorname{span}\{K\}^{\perp}$, the definition (43) of $n_0$, and the fact that integrating the left-hand side of (54) over $[0, 1]$ gives (47), one can apply the pseudo-inverse $\bar{\mathcal{L}}^{\nu^{-1}}$ and obtain

$$ \rho n^{\perp} = \nabla_{\mathbf{x}} \cdot (\mathbf{U}_{\mathcal{T}} \rho) K - \nabla_{\mathbf{x}} \cdot \left[ \frac{\rho}{2} (\mathbf{u}_T + \mathbf{U}_{\mathcal{T}}) K \right]. $$
In both cases, we have exploited the fact that $p_0$ and $n_0$ do not depend on the time variable, and we have used the PDE (47) for expressing the time derivative $\partial_t \rho$. Substituting the above expressions for $\rho p^\perp$ and $\rho n^\perp$ into (52), using (45) and (46), and rearranging terms one obtains

$$f^{\perp} = \nabla_{\mathbf{x}} \cdot (\mathbf{U}_{\mathcal{T}} \rho) \mathcal{T} - \mathbf{v} \cdot \nabla_{\mathbf{x}} (\rho \mathcal{T}) + \frac{1}{2} \{ \nabla_{\mathbf{x}} \cdot (\mathbf{U}_{\mathcal{T}} \rho) \mathcal{T} - \mathbf{v} \cdot \nabla_{\mathbf{x}} (\rho \mathcal{T}) + \rho(T_K + T) \mathbf{v} \cdot \nabla_{\mathbf{x}} K - \nabla_{\mathbf{x}} \cdot \left[ \frac{\rho}{2} (\mathbf{u}_T + \mathbf{U}_{\mathcal{T}}) K \right] \mathcal{T} \}.$$

Differentiating and using the chain rule one then finds

$$f^{\perp} = -[\mathcal{T}(\mathbf{v} - \mathbf{U}_{\mathcal{T}}) \cdot \nabla_{\mathbf{x}} \rho + (\mathbf{v} \cdot \nabla_{\mathbf{x}} \mathcal{T} - \mathcal{T}\nabla_{\mathbf{x}} \cdot \mathbf{U}_{\mathcal{T}}) \rho] + \frac{1}{2} \{ \nabla_{\mathbf{x}} \cdot (\mathbf{U}_{\mathcal{T}} \rho) \mathcal{T} - \nabla_{\mathbf{x}} \cdot \left[ \frac{\rho}{2} (\mathbf{u}_T + \mathbf{U}_{\mathcal{T}}) K \right] \mathcal{T} + \mathbf{v} \cdot \nabla_{\mathbf{x}} (\rho \mathcal{T}) - \mathbf{v} \cdot \nabla_{\mathbf{x}} (\rho T_K) K \} \quad (55)$$

We now introduce the variance-covariance matrix of $\mathcal{T}$, i.e.

$$\mathbb{D}_{\mathcal{T}} \equiv \mathbb{D}_{\mathcal{T}}[S, S^{\dagger}, S^{\ddagger}] := \int_{\mathcal{V}} \int_{0}^{1} (\mathbf{v} - \mathbf{U}_{\mathcal{T}}) \otimes (\mathbf{v} - \mathbf{U}_{\mathcal{T}}) \mathcal{T}[S, S^{\dagger}](\boldsymbol{\nu}, y) \, dy \, d\boldsymbol{\nu}.$$

Exploiting the regularity of $K$ and $T$ to invert the order of integration over $\mathcal{V} \times [0, 1]$, one sees that

$$\mathbb{D}_{\mathcal{T}} = \int_{\mathcal{V}} \int_{0}^{1} (\mathbf{v} - \mathbf{U}_{\mathcal{T}}) \otimes (\mathbf{v} - \mathbf{U}_{\mathcal{T}}) T[S, S^{\dagger}](\boldsymbol{\nu}; y) K[S^{\ddagger}](y) \, dy \, d\boldsymbol{\nu}. \quad (56)$$

Noting also that $\mathbb{D}_{\mathcal{T}} = \int_{\mathcal{V}} \int_{0}^{1} \mathbf{v} \otimes (\mathbf{v} - \mathbf{U}_{\mathcal{T}}) \mathcal{T}[S, S^{\dagger}](\boldsymbol{\nu}, y) \, dy \, d\boldsymbol{\nu}$ thanks to (48), substituting (55) into (51) and integrating over $\mathcal{V} \times [0, 1]$ yields

$$\partial_t \rho + \nabla_{\mathbf{x}} \cdot [\rho \mathbf{U}_{\mathcal{T}} (1 - \varepsilon \nabla_{\mathbf{x}} \cdot \mathbf{U}_{\mathcal{T}})] = \varepsilon \nabla_{\mathbf{x}} \cdot \nabla_{\mathbf{x}} \cdot (\mathbb{D}_{\mathcal{T}} \rho), \quad (57)$$

where $\mathbf{U}_{\mathcal{T}}$ and $\mathbb{D}_{\mathcal{T}}$ are defined via (49) and (56). The drift and diffusion terms in (57) come from the first line on the right-hand side of (55) through steps analogous to those in (Hillen, 2006), while the second line on the right-hand side of (55) gives

$$\frac{1}{2} [\nabla_{\mathbf{x}} \cdot (\mathbf{U}_{\mathcal{T}} \rho) - \nabla_{\mathbf{x}} \cdot (\mathbf{U}_{\mathcal{T}} \rho)] + \frac{1}{2} \left[ \int_{\mathcal{V}} \int_{0}^{1} K \mathbf{v} \otimes \mathbf{v} (\nabla_{\mathbf{x}} T_K \rho + T_K \nabla_{\mathbf{x}} \rho) \, dy \, d\boldsymbol{\nu} - \int_{\mathcal{V}} \int_{0}^{1} K \mathbf{v} \otimes \mathbf{v} (\nabla_{\mathbf{x}} T_K \rho + T_K \nabla_{\mathbf{x}} \rho) \, dy \, d\boldsymbol{\nu} \right] = 0.$$

Note that, thanks to (29), the PDE (57) is an equation for $\rho = \rho_0$ and not for $\rho_\varepsilon$, and provides a first order in $\varepsilon$ correction to the PDE (47).

**Remark 3.** In the same regime, a more general scaling could be considered by choosing

$$\lambda = \frac{\tilde{\lambda}}{\varepsilon}, \quad \mu = \frac{\tilde{\mu}}{\varepsilon}, \quad \tilde{\mu}, \tilde{\lambda} = \mathcal{O}(1)$$

with $\tilde{\mu} \neq \tilde{\lambda}$, instead of (27). The formal derivation carried out in this section could be adapted to this more general case by noting that, under this scaling,

$$\mathcal{T}[S, S^{\dagger}, S^{\ddagger}](\boldsymbol{\nu}, y) := \frac{K[S^{\ddagger}](y)}{2} \left[ \tilde{\mu} T_K[S, S^{\dagger}, S^{\ddagger}](\boldsymbol{\nu}) + \tilde{\lambda} T[S, S^{\dagger}](\boldsymbol{\nu}; y) \right],$$

while the pseudo-inverse would be $\mathcal{L}^{-1}(\phi) = -(\tilde{\lambda} + \tilde{\mu})\phi$.

# 4 Numerical simulations

In this section, we present the results of numerical simulations of the microscopic model introduced in section 2, which is defined by the system (1)-(2), under the parameter scaling (27), and numerical solutions of the corresponding macroscopic model formally derived in section 3 from the phenotype-structured kinetic equation of the mesoscopic model, that is, the PDE (57) complemented with definitions (49) and (56).
## 4.1 Set-up of numerical simulations of the macroscopic model

### Preliminaries

We pose the PDE (57) on a bounded spatial domain $\Omega$ and complement it with the following zero-flux boundary conditions

$$ \{\rho \mathbf{U}_{\tau} - \varepsilon [\mathbb{D}_{\tau} \nabla_{\mathbf{x}} \rho + \rho (\nabla_{\mathbf{x}} \cdot \mathbb{D}_{\tau} + \mathbf{U}_{\tau} \nabla_{\mathbf{x}} \cdot \mathbf{U}_{\tau})]\} \cdot \mathbf{n} = 0 \quad \text{on } \partial\Omega, \quad (58) $$

where $\mathbf{n}$ is the outer unit normal on $\partial\Omega$. These boundary conditions can be derived from appropriate zero-flux kinetic boundary conditions as similarly done in (Loy and Perthame, 2024). We carry out numerical simulations both in a one-dimensional setting and in a two-dimensional setting. In the one-dimensional setting we take $\mathbf{x} \equiv x \in \Omega := [-L_m, L_M]$, while in the two-dimensional setting we take $\mathbf{x} \equiv (x_1, x_2) \tau \in \Omega := [0, L_M] \times [-L_m, L_M]$, where $L_m, L_M \in \mathbb{R}^+$ with $L_m < L_M$. In particular, coherently with the experimental set-up employed in (Goodman et al., 1989, Figure 1A), we take $L_m = 0.06\text{cm}$ and $L_M = 0.06\text{cm}$. As for the scaling parameter $\varepsilon$ in (57), in the one-dimensional setting we explore the effect of considering different values of this parameter by taking $\varepsilon \in \{10^{-2}, 10^{-3}, 10^{-4}\}$, while in the two-dimensional setting simulations are carried out with $\varepsilon = 10^{-3}$.

Moreover, consistently with subsection 3.2, we focus on the case where the phenotypic transition kernel $K[\mathcal{S}^{\dagger}](y|y')$ satisfies assumption (32) – i.e. $K[\mathcal{S}^{\dagger}](y|y') \equiv K[\mathcal{S}^{\dagger}](y)$ – so that definitions (49) and (56) hold.

Finally, to reproduce the experimental set-up of Figure 1A in (Goodman et al., 1989), an illustration of which is given in Figure 8 in the supplementary material, assuming the concentrations of the adhesive glycoproteins laminin and fibronectin and the ECM density in the system to be fixed, we introduce the non-negative, real functions $C_L(\mathbf{x})$ and $C_F(\mathbf{x})$ to model the concentrations (in non-dimensional form) of laminin and fibronectin, respectively, and the non-negative, real function $M(\mathbf{x})$ to model the ECM density (also in non-dimensional form). We then use these functions to identify $\mathcal{S}$, $\mathcal{S}^{\dagger}$, and $\mathcal{S}^{\ddagger}$ when modelling cell reorientation and phenotypic changes, as detailed below.

### Modelling cell reorientation

We assume, for simplicity, that the reorientation of cell polarity is solely dependent on the sensed ECM density; hence, we set $\mathcal{S}(t, \mathbf{x}) \equiv M(\mathbf{x})$ in (9) and (10). Note that, coherently with assumptions (34), this implies that $\mathcal{S}$ is constant in time. Moreover, we assume that the direction of polarity of the cells is biased towards the direction where they sense a higher ECM density, phenomenon known as haptotaxis (Smith et al., 2004). Under these assumptions, we define the function $b$ in (10) as

$$ b[\mathcal{S}](r, \hat{\mathbf{v}}') \equiv b[M](r, \hat{\mathbf{v}}') := M(\mathbf{x} + r\hat{\mathbf{v}}'). \quad (59) $$

Moreover, we let the cell post-reorientation speed be determined by the amount and strength of adhesion sites sensed by the cells, which we assume to be directly linked to the concentrations of laminin and fibronectin. Hence, we set $\mathcal{S}^{\dagger}(t, \mathbf{x}) \equiv (C_L(\mathbf{x}), C_F(\mathbf{x}))$ in (9) and (11), which implies, coherently with assumptions (34), that $\mathcal{S}^{\dagger}$ is constant in time. In particular, we assume the probability distribution $\psi[\mathcal{S}^{\dagger}](v) \equiv \psi[C_L, C_F](v)$ in (11) to be such that the mean post-reorientation speed acquired by cells at position $\mathbf{x}$, which is defined via (13), is

$$ u_{\psi}[\mathcal{S}^{\dagger}] \equiv u_{\psi}[C_L, C_F] = w_L[C_L, C_F] v^L + w_F[C_L, C_F] v^F. \quad (60) $$

Here $v^F, v^L \in \mathbb{R}_0^+$ with $0 \le v^F \le v^L \le V_{\text{max}}$ denote, respectively, the mean post-reorientation speeds preferentially acquired by cells adhering to laminin and fibronectin. The assumption that $v^F \le v^L$ follows from experimental results in (Goodman et al., 1989). Specifically, in order to match the corresponding cell velocities measured in (Goodman et al., 1989), we take

$$ v^L = 2\text{cm/h} \quad \text{and} \quad v^F = 0.6\text{cm/h}. \quad (61) $$
The weights $w_L[C_L, C_F]$ and $w_F[C_L, C_F]$ in (60) provide a measure of the adhesive strength of cells to laminin and fibronectin, and we define them as the local fractions of these two proteins, i.e.

$$w_L[C_L, C_F] := \frac{C_L}{C_L + C_F}, \qquad w_F[C_L, C_F] := \frac{C_F}{C_L + C_F}. \qquad (62)$$

Furthermore, we define the probability distribution $\psi[S^\dagger](v) \equiv \psi[C_L, C_F](v)$ as a Dirac delta distribution centered in $v = u_\psi[C_L, C_F]$, i.e.

$$\psi[C_L, C_F](v) := \delta(v - u_\psi[C_L, C_F]).$$

Finally, in line with assumptions (6), we define the sensing radius of cells in the phenotypic state $y$ as

$$R(y) := R_{\min} + y(R_{\max} - R_{\min}). \qquad (63)$$

Consistently with (Sen et al., 2009), where the authors observed cells with radii as small as one twentieth and as large as one half of the width of each stripe, we set $R_{\min} = 5 \times 10^{-4} \text{cm}$ and $R_{\max} = 5 \times 10^{-3} \text{cm}$ (see also Figure 8 in the supplementary material for a visual comparison with the physical domain).

## Modelling phenotypic changes

We assume phenotypic changes undergone by the cells to be influenced by the local environmental conditions, which are determined by the concentrations of laminin and fibronectin. Hence, we set $S^\dagger(t, \mathbf{x}) \equiv (C_L(\mathbf{x}), C_F(\mathbf{x}))$ in the phenotypic transition kernel, that is, $K[S^\dagger](y) \equiv K[C_L, C_F](y)$. Again, note that this implies that $S^\dagger$ is constant in time, which is coherent with assumptions (34).

We let $K(y)$ be such that the phenotypic state in which cells enter, on average, at position $\mathbf{x}$ as a result of phenotypic changes, which is defined via (4), is

$$y_K[S^\dagger] \equiv y_K[C_L, C_F] = w_L[C_L, C_F] y_K^L + w_F[C_L, C_F] y_K^F. \qquad (64)$$

Here $y_K^L, y_K^F \in \mathbb{R}_0^+$ with $0 \le y_K^F \le y_K^L \le 1$ are, respectively, the phenotypic states in which, on average, cells are preferentially led to enter by phenotypic changes driven by signalling cascades triggered by laminin and fibronectin binding to cell surface receptors, and the assumption that $y_K^F \le y_K^L$ follows from experimental evidence on cytoskeletal organisation reported in (Goodman et al., 1989). Moreover, the weights $w_L[C_L, C_F]$ and $w_F[C_L, C_F]$, which are defined via (62), provide a measure of the strength of the mechanotransductive signals of laminin and fibronectin.

Specifically, in (64) we choose

$$y_K^F = 0 \qquad (65)$$

to model the scenario where cells on fibronectin are induced to enter phenotypic states corresponding to a stiffer cytoskeleton, and thus a smaller sensing radius (cf. definition (63)), which is consistent with what reported in (Goodman et al., 1989).

Furthermore, we investigate two cases for the phenotypic state in which, on average, cells are preferentially led to enter by phenotypic changes driven by signalling cascade triggered by laminin in (64):

* the case where

$$y_K^L = 0, \qquad (66)$$

which models the scenario where cells on laminin are induced to enter phenotypic states corresponding to a stiffer cytoskeleton, and thus a smaller sensing radius (cf. definition (63));

* the case where

$$y_K^L = 1, \qquad (67)$$

which models the scenario, consistent with experimental evidence of (Goodman et al., 1989), where cells on laminin are induced to enter phenotypic states corresponding to a more flexible cytoskeleton, and thus a larger sensing radius (cf. definition (63)).
Moreover, building on the modelling strategies presented in (Conte and Loy, 2024; Loy and Preziosi, 2020b), we consider two possible definitions of the phenotypic transition kernel $K[C_L, C_F](y)$:

* a Dirac delta distribution centered in $y_K[C_L, C_F]$, i.e.

$$K[C_L, C_F](y) := \delta(y - y_K[C_L, C_F]); \qquad (68)$$

* a unimodal von Mises distribution with concentration parameter $k_y$ (specifically, we take $k_y = 2$) and circular mean $y_K[C_L, C_F]$, i.e.

$$K[C_L, C_F](y) := \frac{1}{I_0(k_y)} \exp \left[ k_y \cos \left( 2\pi(y - y_K[C_L, C_F]) \right) \right], \qquad (69)$$

where $I_0$ is the modified Bessel function of the first kind.

Definition (68) translates into mathematical terms the idea that cells exposed to the same environmental conditions are led by phenotypic changes to enter exactly the same phenotypic state (i.e. there is zero phenotypic variance). On the other hand, definition (69) takes into account the fact that, even if the environmental conditions are the same, due to variability in intra-cellular regulatory dynamics amongst cells, there can be variability in the phenotypic state acquired by the cells undergoing phenotypic changes (i.e. there is non-zero phenotypic variance).

For the sake of brevity, in the remainder of this section:

* in the light of definitions (63) and (65), we will refer to the scenario corresponding to definition (66) as the case of “small sensing radius on fibronectin and laminin”, while the scenario corresponding to definition (67) will be referred to as the case of “small sensing radius on fibronectin and large sensing radius on laminin”;

* the scenario corresponding to definition (68) will be referred to as the case of “zero phenotypic variance”, while we will refer to the scenario corresponding to definition (69) as the case of “non-zero phenotypic variance”.

**Remark 4.** Note that, since numerical simulations are carried out in the framework of assumptions (32) and (34), and thus the normalised distribution of the cell population in the phenotype space is stationary and coincides with $K[S^\dagger](y)$ (cf. (43)), in the case of small sensing radius on fibronectin and laminin the mean phenotypic state of the cell population is expected to be $y_K^F$ over the fibronectin stripe and $y_K^L$ over the laminin stripes. Moreover, in the zero phenotypic variance case all cells are expected to be in the phenotypic state $y_K^F$ over the fibronectin stripe and in the phenotypic state $y_K^L$ over the laminin stripes. On the other hand, in the non-zero phenotypic variance case, the phenotypic states of the cells are expected to be scattered around $y_K^F$ over the fibronectin stripe and around $y_K^L$ over the laminin stripes.

## Definitions of the laminin and fibronectin concentrations and the ECM density

Consistently with the experimental set-up employed in (Goodman et al., 1989, Figure 1A), in the two-dimensional setting, we let laminin and fibronectin be distributed along parallel stripes which run along the $x_2$ direction. In the one-dimensional setting, we let laminin be distributed along the stripe running along the $x$ direction. Full details on the definitions of the laminin and fibronectin concentrations and the ECM density in the one- and two-dimensional settings are provided in Section D, and illustrated in Figure 8, in the supplementary material.

## Initial conditions

We complement the PDE (57) with initial conditions that are in line with the experimental set-up employed in (Goodman et al., 1989, Figure 1A), see also Figure 8 in the supplementary material.
Specifically, in the two-dimensional setting we use the following initial condition

$$ \rho^0(x_1, x_2) := \begin{cases} \bar{\rho}^0 \left( 1 - \left( \frac{-L_m + 2l_m - x_2}{l_m} \right)^2 \right) & -L_m + l_m \le x_2 < -L_m + 2l_m, \\ \bar{\rho}^0 & -L_m + 2l_m \le x_2 \le -l_m, \\ \bar{\rho}^0 \left( 1 - \left( \frac{x_2 + l_m}{l_m} \right)^2 \right) & -l_m < x_2 \le 0, \\ 0 & x_2 < -L_m + l_m \lor x_2 > 0, \end{cases} \quad (70) $$

for all $x_1 \in [0, L_M]$, where $\bar{\rho}^0 \in \mathbb{R}^+$ and $l_m \in (0, L_m)$. Similarly, in the one-dimensional setting we use the following initial condition

$$ \rho^0(x) := \begin{cases} \rho^0 \left( 1 - \left( \frac{-L_m + 2l_m - x}{l_m} \right)^2 \right) & -L_m + l_m \le x < -L_m + 2l_m, \\ \bar{\rho}^0 & -L_m + 2l_m \le x \le -l_m, \\ \bar{\rho}^0 \left( 1 - \left( \frac{x + l_m}{l_m} \right)^2 \right) & -l_m < x \le 0, \\ 0 & x < -L_m + l_m \lor x > 0. \end{cases} \quad (71) $$

In particular, we choose $l_m = 0.01\text{cm}$ and $\bar{\rho}^0 = 2 \times 10^3\text{cells/cm}^2$.

**Numerical methods** Numerical solutions are constructed on a uniform grid by using a mixed finite-difference and finite-volume scheme, as detailed in Section E in the supplementary material. The numerical scheme is implemented in MATLAB® and the source code is available on GitHub<fnref n="4" />.

## 4.2 Set-up of numerical simulations of the microscopic model

For consistency with the macroscopic model (57), which is formally obtained under the parameter scaling (27) and making assumption (32), numerical simulations of the microscopic model defined by the system (1)-(2) are carried out letting assumptions (27) and (32) hold and using a set-up corresponding to the one employed for the macroscopic model. To carry out numerical simulations, we implement a Nanbu-Babovski Monte Carlo scheme, as similarly done in (Conte and Loy, 2024; Loy and Tosin, 2021). We use a sufficiently large number of particles (i.e. $10^6$ particles) and a sufficiently small time-step (i.e. $\Delta t = 10^{-5}$), so as to ensure that the implemented Monte Carlo scheme can provide a good numerical approximation of the solution to the phenotype-structured kinetic equation (26). The Monte Carlo scheme is implemented in MATLAB® and the source code is available on GitHub<fnref n="4" />.

## 4.3 Main results of numerical simulations

### 4.3.1 Comparison between the microscopic and the macroscopic models

The plots in Figure 2 present a comparison between the cell number density obtained from numerical simulations of the microscopic model and the numerical solution of the macroscopic model in the one-dimensional setting corresponding to the set-up described in subsections 4.1 and 4.2, for different values of the scaling parameter $\varepsilon$. These plots demonstrate that there is an excellent quantitative agreement between the two models, for all the values of the scaling parameter $\varepsilon$ here considered. Note that, as it can be expected from the form of the PDE (57), larger values of $\varepsilon$ lead to broader spatial distributions of cells.

<footnote n="4">⁴https://github.com/ChiaraVilla/LorenziEtAl2024Modelling</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 2: Comparison between the microscopic and the macroscopic models in 1D. Plots of the cell number density, normalised with respect to $\bar{\rho}^0$, obtained from 1D numerical simulations of the microscopic model (red dots) and the numerical solution of the corresponding macroscopic model (blue lines), also normalised with respect to $\bar{\rho}^0$, at $t = 0h$, $t = 12h$, and $t = 24h$ (decreasing transparency levels). The microscopic model consists of the system (1)-(2) subject to the parameter scaling (27), while the macroscopic model comprises the PDE (57) complemented with definitions (49) and (56). Numerical simulations are carried out under the set-up detailed in subsections 4.1 and 4.2, in the one-dimensional setting, with definition (68), assumption (67), and the scaling parameter $\varepsilon = 10^{-4}$ (top row), $\varepsilon = 10^{-3}$ (middle row) or $\varepsilon = 10^{-2}$ (bottom row). The variable $x$ in the plots is in units of $10^{-2}\text{cm}$.</figcaption></figure>

An excellent quantitative agreement between the microscopic and the macroscopic models is also observed in the two-dimensional setting corresponding to the set-up described in subsections 4.1 and 4.2, as shown by the plots in Figure 3, which display cross-sections of the cell number density resulting from numerical simulations of the microscopic model alongside cross-sections of the numerical solution of the corresponding macroscopic model, when $\varepsilon = 10^{-3}$. Note that the agreement between the two models deteriorates slightly in the regions where the cell number density undergoes sharper changes, as it can be expected due to the histogram method through which the cell number density for the microscopic model is reconstructed from the results of Monte Carlo simulations.

Taken together, these results validate the formal procedures employed to derive first the mesoscopic model, defined by the phenotype-structured kinetic equation (26), and then the macroscopic model, comprising the PDE (57) complemented with definitions (49) and (56), from the microscopic model consisting of the system (1)-(2) in the regime (27).

### 4.3.2 Comparison between the macroscopic model and experimental results

The plots in Figures 4 and 5 display numerical solutions of the macroscopic model defined by the PDE (57) complemented with definitions (49) and (56), under the set-up detailed in subsection 4.1, in the two-dimensional setting mimicking the experimental set-up of the stripe migration assay of Goodman et al. (1989). In more detail, Figure 4 refers to the case of zero phenotypic variance
<figure><img src="image_3.png" /><figcaption>Figure 3: Comparison between the microscopic and the macroscopic models in 2D. Plots of the cross-sections of the cell number density, normalised with respect to $\bar{\rho}^0$, obtained from 2D numerical simulations of the microscopic model (red dots) and the numerical solution of the corresponding macroscopic model (blue lines), also normalised with respect to $\bar{\rho}^0$, at $t = 0h$, $t = 24h$, and $t = 48h$ (decreasing transparency levels). Longitudinal cross-sections at $x_1 = 0.975 \times 10^{-2}\text{cm}$ (top) and $x_1 = 3.025 \times 10^{-2}\text{cm}$ (bottom) – i.e. over laminin and fibronectin, respectively – are displayed on the left, while transversal cross-sections at $x_2 = -0.225 \times 10^{-2}\text{cm}$ are displayed on the right. The full plot of the numerical solution of the macroscopic model at $t = 48h$ is displayed in Figure 4B. The microscopic model consists of the system (1)-(2) subject to the parameter scaling (27), while the macroscopic model comprises the PDE (57) complemented with definitions (49) and (56). Numerical simulations are carried out under the set-up detailed in subsections 4.1 and 4.2, in the two-dimensional setting, with definition (68), assumption (67), and the scaling parameter $\varepsilon = 10^{-3}$. The variables $x_1$ and $x_2$ in the plots are in units of $10^{-2}\text{cm}$.</figcaption></figure>

(i.e. when the phenotypic transition kernel, $K$, is defined via (68)), while Figure 5 refers to the case of non-zero phenotypic variance (i.e. when $K$ is defined via (69)). Moreover, the results in Figures 4A and 5A are for the case of small sensing radius on fibronectin and laminin (i.e. when the phenotypic states in which, on average, cells are preferentially led to enter by signalling cascades triggered by fibronectin, $y_K^F$, and laminin, $y_K^L$, are defined via (65) and (66)). On the other hand, the results in Figures 4B and 5B are for the case of small sensing radius on fibronectin and large sensing radius on laminin (i.e. when $y_K^F$ is still defined via (65) but $y_K^L$ is defined via (67) instead).

The numerical solutions displayed in Figure 4 and Figure 5 show that there is accumulation of cells on the laminin stripes. This is due to the fact that, since cell reorientation is driven by non-local sensing of the surrounding environment, cells near the boundary between the fibronectin stripe and the laminin stripes will eventually move to the laminin stripes, where their motion can occur at a higher speed (cf. assumption (61)). However, in the case of small sensing radius on fibronectin and large sensing radius on laminin, which is the one consistent with the experimental study in Goodman et al. (1989) (cf. Figure 4B and Figure 5B), accumulation of cells over laminin stripes is more pronounced than in the case of small sensing radius on fibronectin and laminin (cf. Figure 4A and Figure 5A). Furthermore, the part of the cell population adhering to laminin undergoes appreciably faster migration than the part of the cell population adhering to fibronectin, an emergent property that was also experimentally observed in (Goodman et al., 1989).

The numerical solutions displayed in Figure 4 and Figure 5 also indicate that migration of the cell population is slightly faster, both on fibronectin and on laminin, when the phenotypic variance in non-zero (cf. Figure 5) than when there is zero phenotypic variance (cf. Figure 4), as it is especially evident in the case of small sensing radius on fibronectin and laminin (cf. Figure 4A and Figure 5A). This can be explained by noting that in this case cells will have a tendency to
<figure><img src="image_4.png" /><figcaption>Figure 4: **Numerical solutions of the macroscopic model in 2D under zero phenotypic variance.** Plots of the numerical solution of the macroscopic model, normalised with respect to $\bar{\rho}^0$, at $t = 0h$ (left column) and at $t = 48h$ (central and right columns). The macroscopic model comprises the PDE (57) complemented with definitions (49) and (56). Numerical simulations are carried out under the set-up detailed in subsection 4.1, in the two-dimensional setting, with definition (68), assumption (66) (central column) or assumption (67) (right column), and the scaling parameter $\varepsilon = 10^{-3}$. The variables $x_1$ and $x_2$ in the plots are in units of $10^{-2}\text{cm}$. Laminin and fibronectin stripes are highlighted by LN and FN, respectively.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 5: **Numerical solutions of the macroscopic model in 2D under non-zero phenotypic variance.** Plots of the numerical solution of the macroscopic model, normalised with respect to $\bar{\rho}^0$, at $t = 0h$ (left column) and at $t = 48h$ (central and right columns). The macroscopic model comprises the PDE (57) complemented with definitions (49) and (56). Numerical simulations are carried out under the set-up detailed in subsection 4.1, in the two-dimensional setting, with definition (69), assumption (66) (central column) or assumption (67) (right column), and the scaling parameter $\varepsilon = 10^{-3}$. The variables $x_1$ and $x_2$ in the plots are in units of $10^{-2}\text{cm}$. Laminin and fibronectin stripes are highlighted by LN and FN, respectively.</figcaption></figure>
enter the phenotypic state $y = 0$ (which corresponds to a stiff cytoskeleton and the smallest sensing radius, $R_{min}$) both on laminin and on fibronectin. Hence, phenotypic variance will be the sole catalyst for the appearance of cells in phenotypic states $y > 0$ (which correspond to a more flexible cytoskeleton and a larger sensing radius). Cells on the fibronectin stripe which are in such phenotypic states will then reach more easily the laminin stripes, where they will move, on average, at a higher speed, thus promoting the emergence of faster migration of the cell population. This is slightly less apparent in the case of small sensing radius on fibronectin and large sensing radius on laminin. In fact, in this case, if the phenotypic variance is non-zero then the sensing radii of the cells will be more uniformly distributed between $R_{min}$ and $R_{max}$ than when there is zero phenotypic variance.

### 4.3.3 Comparison between the original macroscopic model and a reduced model

<figure><img src="image_6.png" /><figcaption>Figure 6: Numerical solutions of the reduced macroscopic model in 2D under zero phenotypic variance. Plots of the numerical solution of the reduced macroscopic model, normalised with respect to $\bar{\rho}^0$, at $t = 0h$ (left column) and at $t = 48h$ (central and right columns). The reduced macroscopic model comprises the PDE (72) complemented with definition (49) and subject to the boundary conditions (73). Numerical simulations are carried out under the set-up detailed in subsection 4.1, in the two-dimensional setting, with definition (68), assumption (66) (central column) or assumption (67) (right column), and the scaling parameter $\varepsilon = 10^{-1}$. The variables $x_1$ and $x_2$ in the plots are in units of $10^{-2}cm$. Laminin and fibronectin stripes are highlighted by LN and FN, respectively.</figcaption></figure>

In the framework of the macroscopic model defined by the PDE (57), cell dynamics result from the superposition of inhomogeneous, anisotropic diffusion, with diffusion tensor $\varepsilon D_T$, and advective transport, with velocity field $U_T(1 - \varepsilon \nabla_x \cdot U_T)$. In particular, the transport term represents the effect, at the cell-population level, of the interplay between directional migration, which is led by cell reorientation driven by non-local sensing of the surrounding environment, and environment-induced changes in the cytoskeletal structure. In order to assess the robustness of the patterns of cell migration produced by the macroscopic model, we investigate whether features qualitatively similar to those displayed by the numerical solutions presented in Figures 4 and 5 can emerge when the transport term in (57) is left unchanged whilst the inhomogeneous, anisotropic diffusion term is replaced by a linear diffusion term with diffusion coefficient $\varepsilon$. Hence, we consider the reduced macroscopic model given by the following PDE

$$ \partial_t \rho + \nabla_x \cdot [\rho U_T (1 - \varepsilon \nabla_x \cdot U_T)] = \varepsilon \Delta_x \rho, \quad (72) $$
<figure><img src="image_7.png" /><figcaption>Figure 7: Numerical solutions of the reduced macroscopic model in 2D under non-zero phenotypic variance. Plots of the numerical solution of the reduced macroscopic model, normalised with respect to $\bar{\rho}^0$, at $t = 0h$ (left column) and at $t = 48h$ (central and right columns). The reduced macroscopic model comprises the PDE (72) complemented with definition (49) and subject to the boundary conditions (73). Numerical simulations are carried out under the set-up detailed in subsection 4.1, in the two-dimensional setting, with definition (68), assumption (66) (central column) or assumption (67) (right column), and the scaling parameter $\epsilon = 10^{-1}$. The variables $x_1$ and $x_2$ in the plots are in units of $10^{-2}cm$. Laminin and fibronectin stripes are highlighted by LN and FN, respectively.</figcaption></figure>

where $U_{\tau}$ is still defined via (49). This simplified model is obtained from the original model by replacing the variance-covariance matrix of $T$, $D_{\tau}$, in the PDE (57) by the identity matrix. In analogy with what done to solve numerically the PDE (57), we pose the PDE (72) on the spatial domain $\Omega$ and complement it with the following zero-flux boundary conditions

$$[\rho U_{\tau} - \epsilon (\nabla_{x}\rho + \rho U_{\tau}\nabla_{x} \cdot U_{\tau})] \cdot \mathbf{n} = 0 \quad \text{on } \partial\Omega. \qquad (73)$$

Moreover, we carry out numerical simulations under the set-up detailed in subsection 4.1, in the two-dimensional setting mimicking the experimental set-up of the stripe migration assay from (Goodman et al., 1989). Numerical solutions are constructed using methods similar to those employed to solve numerically the PDE (57), which are detailed in Section E in the supplementary material.

Numerical solutions of the reduced macroscopic model (72) in the cases of zero and non-zero phenotypic variance are displayed in Figure 6 and Figure 7, respectively. These numerical solutions display the same key qualitative features as those of the numerical solutions in Figures 4 and 5 (i.e. emergence of higher cell density and faster cell migration on laminin stripes), which testifies to the robustness of the patterns of cell migration produced by the macroscopic model.

# 5 Discussion and research perspectives

We developed a modelling framework for migration of heterogeneous cell populations driven by non-local environmental sensing. In the vein of previous work on phenotype-structured models of cell movement reviewed in Lorenzi et al. (2024b), this framework generalises the modelling approach proposed in Loy and Preziosi (2020a,b) by incorporating a continuous phenotype structure that makes it possible to take into account: inter-cellular variability in structural properties of the cytoskeleton, and thus in the length of the surface projections through which the cells sense
the surrounding environment; the occurrence of environment-induced changes in structural prop-
erties of the cytoskeleton, which result in dynamical changes in the length of the cells' surface
projections.

We started by formulating a microscopic model in which single cell dynamics are described by means of stochastic processes, which represent cell movement and environment-induced changes in structural properties of the cytoskeleton. Through a limiting procedure, we formally derived a phenotype-structured kinetic equation that governs the dynamics of the cell distribution in the phase and phenotype spaces, which constitutes the mesoscopic counterpart of the microscopic model. From the mesoscopic model we formally derived a closed PDE for the cell density, which defines the corresponding macroscopic model. Such a PDE comprises an inhomogeneous, anisotropic diffusion term and an advective transport term, with the latter modelling the interplay between environment-induced changes in the cytoskeletal structure and directional migration led by cell reorientation, which is driven by non-local sensing of the surrounding environment.

We compared numerical solutions of the PDE defining the macroscopic model and the results of Monte Carlo simulations of the microscopic model, and showed that there is an excellent quantitative agreement between them, under a simulation set-up which reproduces the experimental set-up of the cell locomotion assays of Goodman et al. (1989). We also showed that numerical solutions of the macroscopic model recapitulate qualitative features of experimental observations presented in Goodman et al. (1989) by demonstrating that the interplay between cell reorientation driven by non-local sensing of the surrounding environment and environment-induced changes in the cytoskeletal structure lead to faster migration of cells adhering to laminin stripes. To corroborate the robustness of the patterns of cell migration emerging in numerical solutions of the macroscopic model, we showed that qualitatively similar patterns are displayed by numerical solutions of a reduced model, wherein the advective transport term is left unchanged while the inhomogeneous, anisotropic diffusion term is replaced by a linear diffusion term.

We conclude with an outlook on possible research perspectives. Informed by the numerical solutions presented in Figures 4-7, it would be interesting to investigate the existence of travelling wave solutions for the PDE (57) or the PDE (72) in the form of travelling pulses. Moreover, we carried out numerical simulations under assumption (32) on the phenotypic transition kernel, which facilitates the explicit identification of $U_\tau$ and $D_\tau$ thanks to the result (33), but it would also be relevant to explore how the patterns of cell migration presented here may change when this assumption is relaxed. Furthermore, while our theoretical study has eschewed specific mechanisms driving phenotypic changes, consistently with the fact that these are not detailed in the experimental study of Goodman et al. (1989), modelling more precisely the effect of chemical or mechanical external stimuli that drive cytoskeletal changes would be another avenue for future research. Although, guided by the experiments of Goodman et al. (1989), in this work we focused on the microscopic model under the parameter scaling (27) and the corresponding macroscopic model defined by the PDE (57), another track to follow would be to integrate the results of our simulation-based study with Monte Carlo simulations of the microscopic model under different parameter scalings, and complement these with numerical solutions of the corresponding macroscopic models. Finally, while our attention has been focused on the microscopic and macroscopic models, it would certainly be relevant to consider also the mesoscopic model defined by the phenotype-structured kinetic equation (26). In this regard, it would be interesting to investigate whether such a kinetic equation admits travelling wave solutions that exhibit phenotype structuring, whereby cells with different cytoskeleton properties dominate different parts of the wave (Lorenzi et al., 2024b); for this, combining techniques similar to those employed in Bouin and Caillerie (2019); Bouin and Calvez (2014); Bouin et al. (2012, 2015) and Freingruber et al. (2025); Lorenzi et al. (2024a); Lorenzi and Painter (2022); Lorenzi et al. (2022) may prove useful. This would allow for further investigation into how phenotypic heterogeneity affects migration of cell populations driven by non-local environmental sensing.
A Formal derivation of the mesoscopic model

Starting from system (1) complemented with relations (2), we formally derive the mesoscopic counterpart of the microscopic model presented in (2). The derivation method relies on classical limiting procedures of kinetic theory for multi-agent systems, as detailed, in the context of velocity-jump processes for cell migration, in Conte and Loy (2024).

Since the components of the quadruple $(\mathbf{X}_{t+\Delta t}, V_{t+\Delta t}, \hat{\mathbf{V}}_{t+\Delta t}, Y_{t+\Delta t})$ are given by the system (1), for any observable $\phi: \mathbb{R}^d \times \mathcal{V} \times [0, 1] \to \mathbb{R}$, in the asymptotic regime $\Delta t \to 0^+$, the expectation

$$ \langle \phi (\mathbf{X}_t, V_t, \hat{\mathbf{V}}_t, Y_t) \rangle = \int_{\mathbb{R}^d} \int_{\mathcal{V}} \int_0^1 \phi(\mathbf{x}, \mathbf{v}, y) f(t, \mathbf{x}, \mathbf{v}, y) \, dy \, d\mathbf{v} \, dx \quad (74) $$

formally satisfies, see for example Pareschi and Toscani (2013), the following differential equation

$$ \begin{aligned} \frac{d}{dt} \langle \phi (\mathbf{X}_t, V_t, \hat{\mathbf{V}}_t, Y_t) \rangle + \nabla_{\mathbf{x}} \cdot \langle V_t \hat{\mathbf{V}}_t \phi(\mathbf{X}_t, V_t, \hat{\mathbf{V}}_t, Y_t) \rangle = & \\ & \mu \langle \phi (\mathbf{X}_t, V'_t, \hat{\mathbf{V}}'_t, Y_t) - \phi (\mathbf{X}_t, V_t, \hat{\mathbf{V}}_t, Y_t) \rangle \\ & + \lambda \langle \phi (\mathbf{X}_t, V_t, \hat{\mathbf{V}}_t, Y'_t) - \phi (\mathbf{X}_t, V_t, \hat{\mathbf{V}}_t, Y_t) \rangle \end{aligned} \quad (75) $$

with

$$ \begin{aligned} \langle \phi (\mathbf{X}_t, V'_t, \hat{\mathbf{V}}'_t, Y_t) \rangle &= \int_{\mathbb{R}^d} \int_{\mathcal{V}} \int_0^1 \int_{\mathcal{V}} \Psi[\mathcal{S}^\dagger](v'; \hat{\mathbf{v}}', y) \mathcal{B}[\mathcal{S}](\hat{\mathbf{v}}'; y) \phi(\mathbf{x}, \mathbf{v}', y) \, d\mathbf{v}' f(t, \mathbf{x}, \mathbf{v}, y) \, dy \, d\mathbf{v} \, dx \\ &= \int_{\mathbb{R}^d} \int_{\mathcal{V}} \int_0^1 \int_{\mathcal{V}} T[\mathcal{S}, \mathcal{S}^\dagger](\mathbf{v}'; y) \phi(\mathbf{x}, \mathbf{v}', y) \, d\mathbf{v}' f(t, \mathbf{x}, \mathbf{v}, y) \, dy \, d\mathbf{v} \, dx, \end{aligned} \quad (76) $$

where the kernel $T$ is defined via (9), and

$$ \langle \phi (\mathbf{X}_t, V_t, \hat{\mathbf{V}}_t, Y'_t) \rangle = \int_{\mathbb{R}^d} \int_{\mathcal{V}} \int_0^1 \int_0^1 K[\mathcal{S}^\dagger](y'|y) \phi(\mathbf{x}, \mathbf{v}, y') \, dy' f(t, \mathbf{x}, \mathbf{v}, y) \, dy \, d\mathbf{v} \, dx, \quad (77) $$

with the kernel $K$ satisfying assumptions (4). By using (74), (76), and (77), we rewrite (75) as

$$ \begin{aligned} \frac{d}{dt} \int_{\mathbb{R}^d} \int_{\mathcal{V}} \int_0^1 \phi(\mathbf{x}, \mathbf{v}, y) f(t, \mathbf{x}, \mathbf{v}, y) \, dy \, d\mathbf{v} \, dx + \int_{\mathbb{R}^d} \int_{\mathcal{V}} \int_0^1 \phi(\mathbf{x}, \mathbf{v}, y) v \hat{\mathbf{v}} f(t, \mathbf{x}, \mathbf{v}, y) \, dy \, d\mathbf{v} \, dx = \\ \mu \int_{\mathbb{R}^d} \int_{\mathcal{V}} \int_0^1 \left( \int_{\mathcal{V}} T[\mathcal{S}, \mathcal{S}^\dagger](\mathbf{v}'; y) \phi(\mathbf{x}, \mathbf{v}', y) \, d\mathbf{v}' - \phi(\mathbf{x}, \mathbf{v}, y) \right) f(t, \mathbf{x}, \mathbf{v}, y) \, dy \, d\mathbf{v} \, dx \\ + \lambda \int_{\mathbb{R}^d} \int_{\mathcal{V}} \int_0^1 \left( \int_0^1 K[\mathcal{S}^\dagger](y'|y) \phi(\mathbf{x}, \mathbf{v}, y') \, dy' - \phi(\mathbf{x}, \mathbf{v}, y) \right) f(t, \mathbf{x}, \mathbf{v}, y) \, dy \, d\mathbf{v} \, dx. \end{aligned} $$

Then, using the fact that

$$ \int_{\mathcal{V}} \int_0^1 \int_{\mathcal{V}} T[\mathcal{S}, \mathcal{S}^\dagger](\mathbf{v}'; y) \phi(\mathbf{x}, \mathbf{v}', y) \, d\mathbf{v}' f(t, \mathbf{x}, \mathbf{v}, y) \, dy \, d\mathbf{v} = \\ \int_{\mathcal{V}} \int_0^1 \int_{\mathcal{V}} T[\mathcal{S}, \mathcal{S}^\dagger](\mathbf{v}; y) \phi(\mathbf{x}, \mathbf{v}, y) \, d\mathbf{v} f(t, \mathbf{x}, \mathbf{v}', y) \, dy \, d\mathbf{v}' $$

and

$$ \int_{\mathcal{V}} \int_0^1 \int_0^1 K[\mathcal{S}^\dagger](y'|y) \phi(\mathbf{x}, \mathbf{v}, y') \, dy' f(t, \mathbf{x}, \mathbf{v}, y) \, dy \, d\mathbf{v} = \\ \int_{\mathcal{V}} \int_0^1 \int_0^1 \int_0^1 K[\mathcal{S}^\dagger](y|y') \phi(\mathbf{x}, \mathbf{v}, y) \, dy f(t, \mathbf{x}, \mathbf{v}, y') \, dy' \, d\mathbf{v} $$

along with definition (21) of the number density of cells in the phenotypic state $y$, $\rho(t, \mathbf{x})n(t, \mathbf{x}, y)$, and the notation $\mathbf{v} = v\hat{\mathbf{v}}$, choosing $\phi(\mathbf{x}, \mathbf{v}, y) := \xi(\mathbf{x})\varphi(\mathbf{v}, y)$, where $\xi(\mathbf{x})$ and $\varphi(\mathbf{v}, y)$ are test
functions, we obtain a weak formulation (in the physical space) of the following weak form of the phenotype-structured kinetic equation (26) for the distribution function $f(t, \mathbf{x}, \mathbf{v}, y)$

$$ 
\left\{ 
\begin{aligned}
& \int_{\mathcal{V}} \int_{0}^{1} \varphi(\mathbf{v}, y) \partial_t f \, dy \, d\mathbf{v} + \nabla_{\mathbf{x}} \cdot \int_{\mathcal{V}} \int_{0}^{1} \varphi(\mathbf{v}, y) \mathbf{v} \, f \, dy \, d\mathbf{v} = \\
& \qquad \mu \int_{\mathcal{V}} \int_{0}^{1} \varphi(\mathbf{v}, y) (\rho T[\mathcal{S}, \mathcal{S}^{\dagger}](\mathbf{v}; y) n - f) \, dy \, d\mathbf{v} \\
& \qquad + \lambda \int_{\mathcal{V}} \int_{0}^{1} \varphi(\mathbf{v}, y) \left( \int_{0}^{1} K[\mathcal{S}^{\dagger}](y|y') f(t, \mathbf{x}, \mathbf{v}, y') dy' - f \right) dy \, d\mathbf{v}, \\
& n(t, \mathbf{x}, y) := \frac{1}{\rho(t, \mathbf{x})} \int_{\mathcal{V}} f(t, \mathbf{x}, \mathbf{v}, y) \, d\mathbf{v},
\end{aligned}
\right. 
\quad (78) 
$$

which in the strong form is (26).

## B Additional considerations on macroscopic models

From the weak form (78) of the phenotype-structured kinetic equation for $f$, under suitable choices of the test function $\varphi(\mathbf{v}, y)$, one can formally derive the governing equations for the moments of $f$, which provide a macroscopic counterpart of the underlying microscopic model. In summary:

* choosing $\varphi(\mathbf{v}, y) := \mathbf{1}_{\mathcal{V}}(\mathbf{v}) \zeta(y)$, where $\mathbf{1}_{(\cdot)}$ is the indicator function of the set $(\cdot)$ and $\zeta(y)$ is a test function, we obtain a weak formulation of the following governing equation for the number density of cells in the phenotypic state $y$, which is defined via (21),

$$ \partial_t (\rho n) + \nabla_{\mathbf{x}} \cdot (\rho n \mathbf{u}) = \lambda \rho \left( \int_0^1 K[\mathcal{S}^{\dagger}](y|y') n(t, \mathbf{x}, y') dy' - n \right), \quad (79) $$

where $\mathbf{u}$ is the mean velocity of cells in the phenotypic state $y$, which is defined via (22);

* choosing $\varphi(\mathbf{v}, y) := \mathbf{v} \zeta(y)$, with $\mathbf{v} = v \hat{\mathbf{v}}$, we obtain a weak formulation of the following governing equation for the momentum of cells in the phenotypic state $y$, which is defined via (23),

$$ 
\begin{aligned}
\partial_t (\rho n \mathbf{u}) + \nabla_{\mathbf{x}} \cdot (\rho n \mathbf{u} \otimes \mathbf{u} + \rho n \mathbf{d}) ={}& \mu \rho (n \mathbf{u}_T - n \mathbf{u}) \\
& + \lambda \rho \left( \int_0^1 K[\mathcal{S}^{\dagger}](y|y') n(t, \mathbf{x}, y') \mathbf{u}(t, \mathbf{x}, y') dy' - n \mathbf{u} \right),
\end{aligned} 
\quad (80) 
$$

where $\mathbf{u}_T$ is defined via (8) and $\rho n \mathbf{d}$ is the second order tensor of cells in the phenotypic state $y$, which is defined via (24);

* choosing $\varphi(\mathbf{v}, y) := \mathbf{1}_{\mathcal{V}}(\mathbf{v}) \mathbf{1}_{[0,1]}(y)$, we obtain the following governing equation for the cell number density, which is defined via (16),

$$ \partial_t \rho + \nabla_{\mathbf{x}} \cdot (\rho \mathbf{U}) = 0, \quad (81) $$

where $\mathbf{U}$ is the mean velocity of the cell population, which is defined via (18);

* choosing $\varphi(\mathbf{v}, y) := \mathbf{v} \mathbf{1}_{[0,1]}(y)$, we obtain the following governing equation for the momentum corresponding to the mean velocity of the cell population, which is defined via (19),

$$ \partial_t (\rho \mathbf{U}) + \nabla_{\mathbf{x}} \cdot (\rho \mathbf{U} \otimes \mathbf{U} + \rho \mathbf{D}) = \mu (\rho \mathbf{U}_T - \rho \mathbf{U}), \quad (82) $$

where

$$ \mathbf{U}_T \equiv \mathbf{U}_T(t, \mathbf{x}) := \int_0^1 \mathbf{u}_T[\mathcal{S}, \mathcal{S}^{\dagger}](y) n(t, \mathbf{x}, y) \, dy, \quad (83) $$

with $\mathbf{u}_T$ defined via (8), and $\rho \mathbf{D}$ is the second order tensor of the cell population, which is defined via (20).
In order to obtain a closed macroscopic model, we can consider appropriately rescaled versions of the system (79)-(82) corresponding to different biological scenarios. We first focus on biological scenarios where phenotypic changes and cell reorientation occur on different time scales. Under these scenarios, introducing a small parameter $\epsilon \in \mathbb{R}^{+}$, we consider the case where phenotypic changes occur more frequently than cell reorientation by assuming

$$ \lambda = \frac{1}{\epsilon}, \quad \mu = O(1) \quad \text{for } \epsilon \to 0^{+}, \qquad (84) $$

and the opposite case, i.e. we alternatively assume

$$ \mu = \frac{1}{\epsilon}, \quad \lambda = O(1) \quad \text{for } \epsilon \to 0^{+}. \qquad (85) $$

Then we consider biological scenarios where phenotypic changes and cell reorientation occur on similar time scales, which are faster than the time scale of collective spatial dynamics of the cells, i.e. we assume

$$ \lambda = \frac{1}{\epsilon}, \quad \mu = \frac{1}{\epsilon}. \qquad (86) $$

Throughout this section, we denote by $f$ the limit of $f_{\epsilon}$ (i.e. the solution to the phenotype-structured kinetic equation (26) under one of the scalings (84)-(86)) as $\epsilon \to 0^{+}$, with the corresponding macroscopic quantities being defined via (16)-(20) and (21)-(24). Moreover, we refer the reader to (2) for the definition of $g_1[S^{\ddagger}](y)$.

**Closure of the macroscopic model under the scaling (84)** Under the scaling (84), in the asymptotic regime $\epsilon \to 0^{+}$:

* from (79) one formally finds

$$ n(t, \mathbf{x}, y) = g_1[S^{\ddagger}](y); \qquad (87) $$

* from (81)-(83), using the expression for $n$ given by (87), one formally obtains the following system

$$ \left\{ \begin{array}{l} \partial_t \rho + \nabla_{\mathbf{x}} \cdot (\rho \mathbf{U}) = 0, \\ \\ \partial_t (\rho \mathbf{U}) + \nabla_{\mathbf{x}} \cdot (\rho \mathbf{U} \otimes \mathbf{U} + \rho \mathbb{D}) = \mu (\rho \mathbf{U}_T - \rho \mathbf{U}), \end{array} \right. \qquad (88) $$

where

$$ \mathbf{U}_T \equiv \mathbf{U}_T[\mathcal{S}, \mathcal{S}^{\dagger}, \mathcal{S}^{\ddagger}] := \int_0^1 \mathbf{u}_T[\mathcal{S}, \mathcal{S}^{\dagger}](y) g_1[\mathcal{S}^{\ddagger}](y) \, dy \qquad (89) $$

with $\mathbf{u}_T$ being defined via (8).

In order to close the system (88), we need to find a closed-form expression for $\mathbb{D}$. To this end, building on a moment closure method that is commonly used for transport models of cell migration – see, for instance, the review Hillen and Painter (2013) – we make the ansatz

$$ p(t, \mathbf{x}, \boldsymbol{\nu}) = \int_{0}^{1} T[\mathcal{S}, \mathcal{S}^{\dagger}](\boldsymbol{\nu}; y) g_{1}[\mathcal{S}^{\ddagger}](y) \, dy, \qquad (90) $$

which follows from assuming that the normalised distribution of the cell population in the phase space, $p$ is fully determined by the environmental conditions. Substituting this ansatz into (20) gives the following closed-form expression for the variance-covariance matrix

$$ \mathbb{D} \equiv \mathbb{D}[\mathcal{S}, \mathcal{S}^{\dagger}, \mathcal{S}^{\ddagger}] := \int_{\mathcal{V}} (\mathbf{v} - \mathbf{U}) \otimes (\mathbf{v} - \mathbf{U}) \int_{0}^{1} T[\mathcal{S}, \mathcal{S}^{\dagger}](\boldsymbol{\nu}; y) g_{1}[\mathcal{S}^{\ddagger}](y) \, dy \, d\boldsymbol{\nu}, \qquad (91) $$

which makes it possible to close system (88).
**Closure of the macroscopic model under the scaling (85)** Under the scaling (85), in the asymptotic regime $\epsilon \to 0^+$:

• from (80) and (82) one formally finds, respectively,

$$ \mathbf{u}(t, \mathbf{x}, y) = \mathbf{u}_T[\mathcal{S}, \mathcal{S}^\dagger](y) \qquad (92) $$

and

$$ \mathbf{U}(t, \mathbf{x}) = \mathbf{U}_T(t, \mathbf{x}), \qquad (93) $$

where $\mathbf{u}_T$ and $\mathbf{U}_T$ are defined via (8) and (83);

• from (79) and (81), using the expressions for $\mathbf{u}$ and $\mathbf{U}$ given by (92) and (93), one formally obtains the following closed system

$$ \left\{ \begin{aligned} \partial_t (\rho n) + \nabla_{\mathbf{x}} \cdot (\rho n \mathbf{u}_T) &= \lambda \rho \left( \int_0^1 K[\mathcal{S}^\dagger](y|y') n(t, \mathbf{x}, y') dy' - n \right), \\ \partial_t \rho + \nabla_{\mathbf{x}} \cdot (\rho \mathbf{U}_T) &= 0. \end{aligned} \right. \qquad (94) $$

**Closure of the macroscopic model under the scaling (86)** Under the scaling (86), in the asymptotic regime $\epsilon \to 0^+$:

• from (79) and (82) one formally finds, respectively, (87) and (93), where $\mathbf{U}_T$ is defined via (89);

• from (81), using the expression for $\mathbf{U}$ given by (93) and (89), one formally obtains

$$ \partial_t \rho + \nabla_{\mathbf{x}} \cdot (\rho \mathbf{U}_T) = 0, \qquad (95) $$

with $\mathbf{U}_T$ being defined via (89).

**Remark 5.** Note that, under the scaling (86), in the asymptotic regime $\epsilon \to 0^+$ from (78) one formally finds the expression for $p$ given by (90) – i.e. if assumptions (86) hold then (90) is the actual expression for $p$ and not only a closure approximation – which makes it possible to compute explicitly all moments of $p$, including the average velocity and the variance-covariance matrix defined via (18) and (20). Moreover, note that, under the scaling (85), making the ansatz (87) – i.e. assuming that the normalised distribution of the cell population in the phenotype space $n$ is fully determined by the local environmental conditions – one could formally reduce the model (94) to the model (95).

**Hyperbolic limit of the phenotype-structured kinetic equation** Under the scaling (86), it is also possible to derive a more accurate closed macroscopic model through the hyperbolic limit of the phenotype-structured kinetic equation (26). This model consists of a first order in $\epsilon$ correction to the PDE (95) (i.e. the PDE (47)), which is derived by means of an asymptotic procedure whereby one starts from the phenotype-structured kinetic equation (26) under the scaling (86) and makes the Chapman-Enskog expansion (28) for $f_\epsilon$ – see (3.2).

## C Considerations about the scaling (27)

Denoting by $t_0$ and $L$ some characteristic time and length scales of the system and by $V$ and $\bar{\rho}$ some reference values of the cell speed and the cell number density, letting

$$ t \to \frac{t}{t_0}, \quad \mathbf{x} \to \frac{\mathbf{x}}{L}, \quad \mathbf{v} \to \frac{\mathbf{v}}{V}, \quad \rho \to \frac{\rho}{\bar{\rho}}, \quad f \to \frac{f}{\bar{\rho}/V^d} $$

and

$$ T[\mathcal{S}, \mathcal{S}^\dagger](\boldsymbol{\nu}; y) \to \frac{T[\mathcal{S}, \mathcal{S}^\dagger](\boldsymbol{\nu}; y)}{V^d}, $$
we obtain the following non-dimensionalised form of the phenotype-structured kinetic equation (26) for the distribution $f(t, \mathbf{x}, \mathbf{v}, y)$

$$ \text{St} \partial_t f + \mathbf{v} \cdot \nabla_{\mathbf{x}} f = \frac{1}{\text{Kn}_{\nu}} (\rho T[\mathcal{S}, \mathcal{S}^{\dagger}](\mathbf{v}; y) n - f) + \frac{1}{\text{Kn}_y} \left( \int_0^1 K[\mathcal{S}^{\dagger}](y|y') f(t, \mathbf{x}, \mathbf{v}, y') dy' - f \right), \quad (96) $$

where the Strouhal number, St, and the Knudsen numbers, Kn$_{ν}$ and Kn$_{y}$, are defined as

$$ \text{St} := \frac{L}{Vt_0}, \quad \text{Kn}_{\nu} := \frac{V}{L\mu}, \quad \text{Kn}_y := \frac{V}{L\lambda}. $$

Consistently with the experiments corresponding to Figure 1A in Goodman et al. (1989), we choose $t_0 = 48h$, $V = 0.5\text{cm/h}$, $L = 10\text{cm}$, and $\mu = 0.5\text{h}^{-1}$, and we also consider phenotypic changes and spatial movement of cells to occur on similar time scales (i.e. $\lambda \approx \mu$). Hence

$$ \frac{L}{Vt_0} \approx 1, \quad \frac{V}{L\mu} = 10^{-1}, \quad \frac{V}{L\mu} \approx \frac{V}{L\lambda}. $$

It is then natural to introduce a small parameter $\varepsilon \in \mathbb{R}^+$ and set

$$ \text{St} = 1, \quad \text{Kn}_{\nu} = \varepsilon, \quad \text{Kn}_y = \varepsilon. $$

Substituting into (96) yields

$$ \partial_t f + \mathbf{v} \cdot \nabla_{\mathbf{x}} f = \frac{1}{\varepsilon} (\rho T[\mathcal{S}, \mathcal{S}^{\dagger}](\mathbf{v}; y) n - f) + \frac{1}{\varepsilon} \left( \int_0^1 K[\mathcal{S}^{\dagger}](y|y') f(t, \mathbf{x}, \mathbf{v}, y') dy' - f \right), $$

which corresponds to (26) under the scaling (27), thus indicating that, in the light of the experiments of Goodman et al. (1989), it is reasonable to consider the scaling (27).

## D Definitions of the laminin and fibronectin concentrations and the ECM density used in numerical simulations

Consistently with the experimental set-up employed in (Goodman et al., 1989, Figure 1A), illustrated in Figure 8, in the two-dimensional setting, we let laminin and fibronectin be distributed along parallel stripes which run along the $x_2$ direction. Specifically, we use the following definitions, corresponding to the situation where two stripes of laminin are separated by one stripe of fibronectin:

$$ C_L(x_1, x_2) := \begin{cases} \bar{C}_L & x_1 \in [0, \frac{1}{3}L_M] \cup (\frac{2}{3}L_M, L_M], \quad \forall x_2 \in [-L_m, L_M], \\ 0 & \text{otherwise,} \end{cases} \quad (97) $$

$$ C_F(x_1, x_2) := \begin{cases} \bar{C}_F & x_1 \in (\frac{1}{3}L_M, \frac{2}{3}L_M], \quad \forall x_2 \in [-L_m, L_M], \\ 0 & \text{otherwise,} \end{cases} \quad (98) $$

with $\bar{C}_L, \bar{C}_F \in \mathbb{R}^+$. Moreover, we assume the ECM density to be uniformly increasing in the $x_2$ direction, and thus use the following definition

$$ M(x_1, x_2) := M_{\min} + M_{\text{gr}}(x_2 + L_m) \quad \forall x_1 \in [0, L_M], \quad (99) $$

with $M_{\min}, M_{\text{gr}} \in \mathbb{R}^+$.

Analogously, in the one-dimensional setting we use the following definitions

$$ C_L(x) = \bar{C}_L, \quad C_F(x) = 0, \quad M(x) := M_{\min} + M_{\text{gr}}(x + L_m). \quad (100) $$

Since the laminin and fibronectin concentrations and the ECM density are in non-dimensional form, we choose $\bar{C}_L = 1$, $\bar{C}_F = 1$, $M_{\min} = 0.1$, and $M_{\text{gr}} = 1/\text{cm}$.
<figure><img src="image_8.png" /><figcaption>Figure 8: Initial conditions for the numerical simulations in the two-dimensional (top row, right panel) and one-dimensional (bottom row) settings, reproducing the experimental set-up of Goodman *et al.* Goodman *et al.* (1989). In their *in vitro* stripe assay of cell migration, Goodman *et al.* arranged parallel stripes of laminin (cf. top row, left panel) and fibronectin (top row, central panel), and plated cells (myoblasts) at one end of the stripes, evenly distributing them over laminin and fibronectin (top row, right panel). In this experimental set-up, during cell migration, they observed cells with radius as small as $R_{min} = 5 \times 10^{-4}$ cm and as large as $R_{max} = 5 \times 10^{-3}$ cm, marked in red in the top row, right panel for visual comparison. In the one-dimensional setting, the cells are assumed to be arranged over a laminin stripe. The variables $x_1, x_2$, and $x$ in the plots are in units of $10^{-2}$cm.</figcaption></figure>
## E Numerical methods used for solving the PDE (57)

We begin by rewriting the PDE (57) as

$$ \partial_t \rho + \nabla_x \cdot (\rho \mathbf{U}_T^\varepsilon) = \nabla_x \cdot (\nabla_x \cdot (\mathbb{D}_T^\varepsilon \rho)), \quad (101) $$

where

$$ \mathbf{U}_T^\varepsilon = \mathbf{U}_T (1 - \varepsilon \nabla_x \cdot \mathbf{U}_T) \quad \text{and} \quad \mathbb{D}_T^\varepsilon = \varepsilon \mathbb{D}_T, \quad (102) $$

and, under the set-up of numerical simulations described in (4.1), $\mathbf{U}_T$ and $\mathbb{D}_T$ are defined via (49) and (56) with

$$ S(t, \mathbf{x}) \equiv M(\mathbf{x}), \quad S^\dagger(t, \mathbf{x}) \equiv (C_L(\mathbf{x}), C_F(\mathbf{x})), \quad S^\ddagger(t, \mathbf{x}) \equiv (C_L(\mathbf{x}), C_F(\mathbf{x})). $$

Hence, $\mathbf{U}_T(t, \mathbf{x}) \equiv \mathbf{U}_T(\mathbf{x})$ and $\mathbb{D}_T(t, \mathbf{x}) \equiv \mathbb{D}_T(\mathbf{x})$, which also imply that $\mathbf{U}_T^\varepsilon(t, \mathbf{x}) \equiv \mathbf{U}_T^\varepsilon(\mathbf{x})$ and $\mathbb{D}_T^\varepsilon(t, \mathbf{x}) \equiv \mathbb{D}_T^\varepsilon(\mathbf{x})$.

We discretise the spatial domain with a uniform grid of step $\Delta x = 0.05$ and solve numerically the PDE (101) using an explicit first-order in time mixed finite-difference and finite-volume scheme. We choose the time-step $\Delta t$ such that the following CFL condition is satisfied

$$ \Delta t \le \min \left( \frac{\Delta x}{\max_{\mathbf{x} \in \Omega} \mathbf{U}_T^\varepsilon(\mathbf{x})}, \frac{(\Delta x)^2}{2 \max_{\mathbf{x} \in \Omega} \mathbb{D}_T^\varepsilon(\mathbf{x})} \right) $$

so as to ensure stability of the scheme.

The numerical approximation of the advection term in (101) relies on a finite-volume scheme with a first-order upwind approximation for the advective flux $\rho \mathbf{U}_T^\varepsilon$. This is implemented within a MUSCL scheme Van Leer (1979) by setting the flux limiter function to zero. Since here $\mathbf{U}_T^\varepsilon$ does not depend on $t$, the scheme includes a correction term obtained from employing the average flux between the current time-step and the next one – see Dullemond (2008) for details. The numerical approximation of $\mathbf{U}_T^\varepsilon$ relies on a first-order central finite-difference approximation of the first-order spatial derivatives of $\mathbf{U}_T^\varepsilon$.

The right-hand side of the PDE (101) is treated with second-order central finite-difference approximations of the second-order derivatives of the product $\mathbb{D}_T^\varepsilon \rho$. The resulting scheme requires a three-point stencil for the approximation at each grid point for simulations in 1D, while a nine-point stencil for the approximation at each grid point is required for simulations in 2D – this is due to the presence of second-order mixed derivatives in space, the approximation of which requires the value of $\mathbb{D}_T^\varepsilon \rho$ at the centres of the diagonal neighbouring cells.

The zero-flux boundary conditions (58) are implemented with the use of ghost points.

This scheme requires the approximation of $\mathbf{U}_T$, defined via (49), at the grid cell interfaces and the approximation of $\mathbb{D}_T$, defined via (56), at the grid cell centres. For both terms, we employ a uniform discretisation of step $\Delta y = 0.1$ for the phenotypic domain and a uniform discretisation for the sensing region, which consists of $N\theta = 60$ points for $\theta \in [0, \pi]$, and then approximate integrals using a midpoint double Riemann sum. For consistency with the zero-flux boundary conditions (58), we replace the points within a sensing region that fall outside the spatial domain $\Omega$ with the nearest boundary points. For further details and numerical optimisation of the computation of non-local terms like those in (49) and (56), we refer the interested reader to Gerisch (2010).

Simulations are performed in MATLAB® and the numerical tests carried out include: mass conservation check; employing different flux limiters in the MUSCL scheme (superbee, Lax-Wendroff, minimod, ospre, Koren, MC, and van Leer); investigating numerical convergence by varying grid step for the spatial variable as well as discretisation step for the phenotypic domain and the sensing region; verifying stability by comparison with the numerical solution obtained using the MATLAB function ode45, which is based on an explicit Runge-Kutta method of higher order with time-step dynamically adjusted to control accuracy.
References

E. Bouin and N. Caillerie. Spreading in kinetic reaction-transport equations in higher velocity dimensions. *European Journal of Applied Mathematics*, 30(2):219-247, 2019.

E. Bouin and V. Calvez. Travelling waves for the cane toads equation with bounded traits. *Nonlinearity*, 27(9):2233, 2014.

E. Bouin, V. Calvez, N. Meunier, S. Mirrahimi, B. Perthame, G. Raoul, and R. Voituriez. Invasion fronts with variable motility: phenotype selection, spatial sorting and wave acceleration. *Comptes Rendus Mathématique*, 350(15-16):761-766, 2012.

E. Bouin, V. Calvez, and G. Nadin. Propagation in a kinetic reaction-transport equation: travelling waves and accelerating fronts. *Archive for Rational Mechanics and Analysis*, 217:571-617, 2015.

A. Chauviere, T. Hillen, and L. Preziosi. Modeling the motion of a cell population in the extra-cellular matrix. *Discrete and Contiguous Dynamical Systems*, 2007:250-259, 2007.

L. Chen, K. Painter, C. Surulescu, and A. Zhigun. Mathematical models for cell migration: a non-local perspective. *Philosophical Transactions of the Royal Society B*, 375(1807):20190379, 2020.

M. Conte and N. Loy. A non-local kinetic model for cell migration: a study of the interplay between contact guidance and steric hindrance. *SIAM Journal on Applied Mathematics*, 84(3): S429-S451, 2024.

K. Deimling. *Nonlinear Functional Analysis*. Springer, 1985.

C. Dullemond. Advection algorithms ii: Flux conservation, subgrid models and flux limiters, 2008.

V. Freingruber, T. Lorenzi, K. J. Painter, and M. Ptashnyk. Trait-structured chemotaxis: Exploring ligand-receptor dynamics and travelling wave properties in a keller-segel model. *arXiv:2502.18947*, 2025.

P. Friedl and D. Gilmour. Collective cell migration in morphogenesis, regeneration and cancer. *Nature Reviews Molecular Cell Biology*, 10(7):445-457, 2009.

P. Friedl, K. S. Zänker, and E.-B. Bröcker. Cell migration strategies in 3-d extracellular matrix: Differences in morphology, cell matrix interactions, and integrin function. *Microscopy Research and Technique*, 43(5):369-378, 1998.

A. Gerisch. On the approximation and efficient evaluation of integral terms in pde models of cell adhesion. *IMA Journal of Numerical Analysis*, 30(1):173-194, 2010.

V. Givant-Horwitz, B. Davidson, and R. Reich. Laminin-induced signaling in tumor cells. *Cancer Letters*, 223(1):1-10, 2005.

V. Gkretsi and T. Stylianopoulos. Cell adhesion and matrix stiffness: coordinating cancer cell invasion and metastasis. *Frontiers in Oncology*, 8:145, 2018.

S. L. Goodman, G. Risse, and K. von der Mark. The e8 subfragment of laminin promotes locomotion of myoblasts over extracellular matrix. *The Journal of Cell Biology*, 109(2):799-809, 1989.

T. Hillen. M 5 mesoscopic and macroscopic models for mesenchymal motion. *Journal of Mathematical Biology*, 53(4):585-616, 2006.

T. Hillen and K. J. Painter. *Transport and Anisotropic Diffusion Models for Movement in Oriented Habitats*, pages 177-222. Springer Berlin Heidelberg, Berlin, Heidelberg, 2013.
M. K. Jolly, S. A. Mani, and H. Levine. Hybrid epithelial/mesenchymal phenotype (s): The 'fittest' for metastasis? *Biochimica et Biophysica Acta (BBA)-Reviews on Cancer*, 1870(2): 151-157, 2018.

T. Lorenzi and K. J. Painter. Trade-offs between chemotaxis and proliferation shape the phenotypic structuring of invading waves. *International Journal of Non-Linear Mechanics*, 139:103885, 2022.

T. Lorenzi, B. Perthame, and X. Ruan. Invasion fronts and adaptive dynamics in a model for the growth of cell populations with heterogeneous mobility. *European Journal of Applied Mathematics*, 33(4):766-783, 2022.

T. Lorenzi, F. R. Macfarlane, and K. J. Painter. Derivation and travelling wave analysis of phenotype-structured haptotaxis models of cancer invasion. *European Journal of Applied Mathematics*, pages 1-33, 2024a.

T. Lorenzi, K. J. Painter, and C. Villa. Phenotype structuring in collective cell migration: a tutorial of mathematical models and methods. *arXiv preprint arXiv:2410.13629*, 2024b.

N. Loy and B. Perthame. A Hamilton-Jacobi approach to nonlocal kinetic equations. *Nonlinearity*, 37(10):105019, 2024.

N. Loy and L. Preziosi. Kinetic models with non-local sensing determining cell polarization and speed according to independent cues. *Journal of Mathematical Biology*, 80(1):373-421, 2020a.

N. Loy and L. Preziosi. Modelling physical limits of migration by a kinetic model with non-local sensing. *Journal of Mathematical Biology*, 80(6):1759-1801, 2020b.

N. Loy and A. Tosin. Markov jump processes and collision-like models in the kinetic description of multi-agent systems. *Communications in Mathematical Sciences*, 18(6):1539-1568, 2020.

N. Loy and A. Tosin. Boltzmann-type equations for multi-agent systems with label switching. *Kinetic and Related Models*, 14(5):867-894, 2021.

T. Mitchison and L. Cramer. Actin-based cell motility and cell locomotion. *Cell*, 84(3):371-379, 1996.

S. Miyamoto, B.-Z. KATHZ, R. M. Lafrenie, and K. M. Yamada. Fibronectin and integrins in cell adhesion, signaling, and morphogenesis. *Annals of the New York Academy of Sciences*, 857(1): 119-129, 1998.

H. G. Othmer, S. R. Dunbar, and W. Alt. Models of dispersal in biological systems. *Journal of Mathematical Biology*, 26(3):263-298, 1988.

L. Pareschi and G. Toscani. *Interacting multiagent systems: kinetic equations and Monte Carlo methods*. OUP Oxford, 2013.

P. Rousselle and J. Y. Scoazec. Laminin 332 in cancer: When the extracellular matrix turns signals from cell anchorage to cell movement. In *Seminars in Cancer Biology*, volume 62, pages 149-165. Elsevier, 2020.

S. Sen, A. J. Engler, and D. E. Discher. Matrix strains induced by cells: computing how far cells can feel. *Cellular and Molecular Bioengineering*, 2:39-48, 2009.

J. T. Smith, J. K. Tomfohr, M. C. Wells, T. P. Beebe, T. B. Kepler, and W. M. Reichert. Measurement of cell migration on surface-bound fibronectin gradients. *Langmuir*, 20(19):8279-8286, 2004.

D. W. Stroock. Some stochastic processes which arise from a model of the motion of a bacterium. *Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete*, 28(4):305-315, 1974.
X. Trepat, Z. Chen, and K. Jacobson. Cell migration. *Comprehensive Physiology*, 2(4):2369, 2012.

B. Van Leer. Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method. *Journal of Computational Physics*, 32(1):101–136, 1979.