# A 2D fluid-based model for individual cell crawling trajectories

Christèle Etchegaray, Nicolas Meunier, Raphaël Voituriez

► To cite this version:

Christèle Etchegaray, Nicolas Meunier, Raphaël Voituriez. A 2D fluid-based model for individual cell crawling trajectories. 2025. hal-05102874

HAL Id: hal-05102874

https://inria.hal.science/hal-05102874v1

Preprint submitted on 7 Jun 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
A 2D fluid-based model for individual cell crawling
trajectories

Christèle Etchegaray<fnref n="1" />, Nicolas Meunier<fnref n="2" />, Raphaël Voituriez<fnref n="3" />

Contributing authors: christele.etchegaray@inria.fr;
nicolas.meunier@univ-evry.fr; raphael.voituriez@upmc.fr;

Abstract

This paper is devoted to the analysis and numerical simulations of a 2D model of single-cell crawling migration focusing on the intracellular dynamics. The fluid-based model relies on the interaction between the actin cytoskeleton and a generic molecular regulator of its activity. It writes as a Poisson problem with Darcy law coupled with a nonlocal advection-diffusion problem on the molecular concentration, where the nonlocal term involves boundary values. Existence and uniqueness of resting and traveling stationary states are established in some parameter domains, together with the linear stability in the resting case. Numerical simulations allow to illustrate both cell regimes. Additive stochastic perturbations of the advection field by a colored Gaussian noise allow to generate numerical trajectories ranging from diffusive to ballistic behaviours, in coherence with experimental data and previous modelling from [1]. A positive correlation between velocity and persistence is also recovered. This highlights the model's ability to capture key features of cell migration.

**Keywords:** Cell migration, Darcy law, Nonlinear advection-diffusion equation, Cell trajectories

<footnote n="1">¹Inria, Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251,
Talence, F-33400, France.</footnote>

<footnote n="2">²LaMME, CNRS UMR 8071, Université Évry Val d'Essonne, 23
boulevard de France, Évry, 91 037, France.</footnote>

<footnote n="3">³Institut de Biologie Paris-Seine (IBPS), Laboratoire Jean Perrin
(LJP), Sorbonne Université, CNRS, Paris, F-75005, France.</footnote>
# 1 Introduction

## 1.1 Motivations: a general rule of cell migration

Cell migration is a fundamental process involved in many physiological and pathological functions (immune response, morphogenesis, cancer metastasis, etc. [2, 3]). It is based on a complex intracellular machinery [4]. Therefore, understanding its key features in spite of the variety of cellular behaviours is a challenging task. Recently, a biophysical approach showed that even if different migration modes coexist, cell migration obeys to a very general principle.

In [1], following the First World Cell Race [5], the authors performed and analyzed extensive experimental assays (both in geometrical settings and in migration modes). They found that cell persistence, that is the directional correlation of trajectories, is exponentially correlated to cell velocity: faster cells have a more directional migration. This coupling was denoted by "Universal Coupling between cell Speed and Persistence" (UCSP)). Further experiments suggested that it relies on actin flows, shared by all cell types. Actin filaments are components of the cell cytoskeleton, the dynamic set of biopolymers responsible for the integrity and force generation of the cell. These filaments are polar: they grow by polymerizing at one end, and shrink by depolymerizing at the other end. In migrating cells, polymerizing ends are located near the cell membrane that resists the filaments growth. As a result, in the frame of reference of the cell, actin filaments drift away from the membrane, forming the so-called actin retrograde flows. Large and fast actin flows are known to improve cell polarity. Experiments suggest that this could follow from the advection of polarity cues, molecules involved in the regulation of the cytoskeleton activity. Such molecules could be involved in the regulation of various phenomena like contraction or regulation of the cytoskeleton activity, depending on the cell type.

## 1.2 Physical model for a minimal feedback loop

A physical model based on the feedback loop between actin flows and a generic molecular species was introduced in [1]. The dynamics is described along the 1D cell polarity axis $[0, L]$, $L > 0$, and in the frame of reference of the moving cell. In this model, the molecules are advected by actin flows, and can be activated at the cell membrane. The molecular concentration $c$ follows a conservative reaction-advection-diffusion dynamics, where the advection velocity is the global actin flows velocity norm $u$. This writes

$$ \partial_t c - \partial_x [D \partial_x c - u_t c + \zeta_c] = 0, \quad (1) $$

where $\zeta_c$, a Gaussian white noise in time carrying space correlations, describes the fluctuations in the current of molecules. The concentration in activated molecules at the boundaries follows a Hill response of index $n \in \mathbb{N}^*$ and saturation threshold $C_s$.

The inward actin flows are described by a global actin velocity defined by the polar variable $(U_t)_t = (u_t, \phi_t)_t$, where $\phi$ is the orientation of the polarity axis. It follows a stochastic diffusion, where dependence of the fluctuations on the orientation
is neglected. One has

$$ \left\{ \begin{aligned} du_t &= -\lambda(u_t - u_t^*) dt + \frac{K}{2u_t} dt + \sqrt{K} dB_t^u, \\ d\phi_t &= \frac{\sqrt{K}}{u_t} dB_t^\phi, \end{aligned} \right. $$

with $dB_t^u$ and $dB_t^\phi$ independent gaussian white noises of unit variance, $K$ the amplitude of fluctuations, $\lambda^{-1}$ their typical timescale. The terms involving $u_t^{-1}$ come from the polar formulation of the SDE. The velocity module is assumed to fluctuate around a mean value $u_t^*$, which phenomenologically depends on the difference in concentration of activated molecules between both cell ends:

$$ u_t^* = -\beta \left( \frac{c(t, 0)^n}{C_s^n + c(t, 0)^n} - \frac{c(t, L)^n}{C_s^n + c(t, L)^n} \right), \quad (2) $$

with $\beta$ the tightness of coupling between actin flows and molecular asymmetry. Hence, the steeper the activated concentration gradient across the cell, the faster the mean actin flow. The negative correlation depicts the molecular cue as a rear marker.

Considering the molecular concentration at steady state with respect to the actin velocity, a series of approximations give rise to a closed SDE on $U$, or equivalently on the cell velocity $V := -U$. This problem carries a complex noise term combining molecular and actin contributions. Its analysis led to typical quantities (mean polarization time, mean persistence time) obtained either analytically or numerically. The authors were able to reproduce the exponential correlation of the UCSP law and provided a bifurcation diagram distinguishing Brownian, Persistent, and Intermittent behaviours. It was based on both the coupling parameter $\beta$ and the activation saturation parameter $C_s$.

This work puts to light the powerful idea of a generic mechanism for cell migration, and allows further analysis of cell trajectories from microscopic features. However, the model includes some limitations. The first point is that the molecular dynamics is one dimensional, hence not affected by the angular variations of the polarity axis. The second point is that the actin velocity is a global quantity, whose dynamics is defined as a stochastic diffusion drifting towards its mean value. Its regulation by the molecular specie is defined phenomenologically from the boundary molecular concentration. Lastly, the derivation of the closed SDE on cell velocity relies on formal approximations involving two sources of noise at different time scales. Providing a justification in 2D and elucidating the importance of combining these two sources of stochasticity remain to decipher. With that in mind, we address the two first points by developing a two-dimensional problem based on a fluid description of the cytoskeleton. This yields a fully spatial description of cell polarisation at both molecular and cytoskeletal scales. Their coupling is then justified by the force balance principle applied on the fluid problem. This results in a minimal model of the intracellular physicochemical interactions leading to cell motion. This context is however less suited to address the third point, which should be primarily investigated with a global actin velocity. This will be the object of a future work.
## 1.3 Context and plan of the paper

Existing physical models based on a fluid description of the cytoskeleton mainly consist in free boundary problems, and aim at investigating the cell shape stability. In [6], the authors use a Darcy law for the fluid depicting the actin cytoskeleton and its polar order for a crawling cell. They show that the cytoskeleton’s reorganisation is sufficient to induce cell motion, and investigate symmetry-breaking in the system. In [7], the authors model a cytoskeletal fragment and show the nonlinear instability of the center of mass of the system. They also relate the fragment’s shape to its velocity.

Numerous mathematical models were also developed (see the references in [8] for an overview). In [8], the authors build and study a free boundary model for cell deformation in response to the interaction between the cell and its environment. In [9], the authors write a free boundary Keller-Segel type model for cell crawling without external cues. The model, related to a previous work [10] takes into account the retrograde actin flow as well as cell contraction arising in some cell types. They put to light both steady and non radial stationary solutions, as well as travelling waves stationary solutions.

In our study that follows the work in [11], we focus on the minimal coupling that is shared among cell types. Namely, we consider a polymerization-based feedback loop between the actin retrograde flow and one of its molecular regulator. Our goal is to focus on this interaction that ought to occur for very different cells having different shapes. For this reason, we neglect the variation in cell shape and study the system in a solid discoidal domain.

In the context of polarization and motility of eukaryotic cells on substrates, some convection-diffusion 2d models were designed in [9, 12–17], but either they were mathematically studied in the 1d case in [18–21], or they were only partially mathematically studied in the 2d case [12, 16, 22] (traveling wave solutions) and [23] for existence results.

This paper is organized as follows: in Section 2, we present the model under study, and give some energy estimates. In Section 3, we put to light both resting and traveling stationary regimes under given conditions. In Section 4, we perform numerical simulations of the system using the finite volumes method. The simulations illustrate the two regimes and provide phase diagrams depicting transition between both. In Section 5, we perform numerical simulations of the system perturbed by a colored Gaussian noise. The simulations yield a variety of numerical trajectories that are characterized and discussed in Section 6. Finally, Section 7 summarizes our findings.

# 2 The 2D model and some energy estimates

In this Section, we introduce the model under study and provide related mathematical estimates. Subsection 2.1 is devoted to the model presentation and introduction of the mathematical framework. Then, some estimates are given in Subsection 2.2.
## 2.1 Modelling approach

Denote $\Omega = B(0, 1)$ the solid disk defining the motionless cell in its frame of reference, and $\mathbf{v} : \mathbb{R}_+ \to \mathbb{R}^2$ the cell velocity. The actin cytoskeleton is described by a fluid of velocity $\mathbf{u} : \mathbb{R}_+ \times \Omega \to \mathbb{R}^2$ and pressure $p : \mathbb{R}_+ \times \Omega \to \mathbb{R}$. We consider a molecular specie of cytosolic concentration $c : \mathbb{R}_+ \times \Omega \to \mathbb{R}_+$, diffusing in the cell body, and carried by the actin flows. The problem for the molecular dynamics writes as follows:

$$ \left\{ \begin{array}{ll} \partial_t c + \operatorname{div} (c\mathbf{u} - D\nabla c) = 0 & \text{in } \Omega, \\ (D\nabla c - c\mathbf{u}) \cdot \mathbf{n} = 0 & \text{on } \partial\Omega. \end{array} \right. \qquad (3) $$

In these equations, $D$ is the diffusion coefficient, and the boundary condition ensures that the molecular mass $M$ is preserved. In the following, it is fixed to $M = 1$. The fluid problem consists in the coupling between a Darcy law and a Poisson-Dirichlet problem. It writes

$$ \mathbf{u}(t, \mathbf{x}) = -\nabla p(t, \mathbf{x}) - \mathbf{v}(t) \text{ in } \Omega, \qquad (4) $$

$$ \left\{ \begin{array}{ll} -\Delta p = -1 & \text{in } \Omega, \\ p = -\chi k_p(c) & \text{on } \partial\Omega, \end{array} \right. \qquad (5) $$

where $-\chi k_p$ is a Dirichlet boundary term and $-1$ a sink term quantifying respectively actin polymerization and depolymerization. Together, they are responsible for the generation of an inward flow from the boundary, seen in the pressure gradient in Equation (4). The dependence of $k_p$ on the concentration $c$ describes its regulating effect on the actin polymerization. When $k_p$ is an increasing function of $c$, the accumulation of activated molecules on a part of the boundary locally decreases the pressure, characterizing molecular inhibitors of actin polymerization. Finally, the cell velocity is obtained from the force balance principle and writes

$$ \mathbf{v}(t) = \gamma \int_{\partial\Omega} p \mathbf{n} \, d\sigma \mathbf{x} = -\gamma\chi \int_{\partial\Omega} k_p(c) \mathbf{n} \, d\sigma. \qquad (6) $$

In the case of an increasing function $k_p$, for example for $k_p(c) = c$, the cell velocity rewrites $\mathbf{v}(t) = -\gamma\chi \int_{\partial\Omega} c \mathbf{n} \, d\sigma$ and the molecular specie characterizes the cell rear. Hence, from the modelling viewpoint, the necessary assumptions on $k_p$ are its regularity, being increasing, and such that $k_p(0) = 0$. Let us now introduce the following assumptions on $k_p$:

**Hypotheses 1**

$k_p \in C^3(\mathbb{R}_+)$, and $\lim_{x \to +\infty} k_p(x) = 1$,

$k_p(0) = 0$ and $k_p$ is monotone increasing,

There exists a positive constant $A$ such that for all $x \ge 0$:

$$ xk'_p(x) \le A. $$
A prototype example of a function satisfying Assumptions 1-1-1 is

$$k_p(x) = \frac{x}{c_t + x}, \qquad (7)$$

for $c_t > 0$. In such a case, (A3) holds with $A = 1/4$.

## Discussion

This model represents an attempt to reach a minimal description of the chemical and mechanical machinery in a single crawling cell. Its complete derivation is detailed in Appendix A. The mechanical aspects of the model relate the cytoskeleton dynamics and the cell displacement as follows. The Poisson problem (5) on the pressure presents a nonpositive Dirichlet boundary condition, leading to a negative pressure with larger values at the boundary due to polymerization, and lower values in the cell center, hence driving a pressure gradient towards the cell boundary. The actin velocity given by the Darcy problem (4) involves two components. The local term corresponds to the opposite of the pressure gradient, and depicts the inward actin flow resulting from the polymerization dynamics. The nonlocal term corresponds to the opposite of the cell velocity, and arise from the actin flows friction forces expressed in the frame of reference of the cell. It acts as an additional directional bias on the actin flows resulting from cell motion itself. The cell velocity (6) follows from the force balance principle applied to the whole cell, involving friction forces arising from both the inner actin flows and the cell displacement itself. As a consequence, the cell displacement follows from asymmetry in the boundary pressure. All in all, inhomogeneous polymerization at the boundary drives polarisation of the actin flows from both local and nonlocal effects, mediating a link between cell velocity and actin flows, in coherence with the observation that *faster cells have larger actin flows*. Depolymerization has only a local depolarizing effect.

The molecular species is transported by the actin flows, and is responsible for polymerization inhibition in the pressure boundary condition. As in [1], it is written using a nonlinear Hill function to quantify an instantaneous activation with saturation. The Hill coefficient is taken equal to 1 for simplicity, as it was noticed to have little effect within the range [0.7, 1.5]. This results in a cell velocity directly expressed from the integration of the boundary molecular concentration, so that faster cell velocities follow from steeper concentration gradient across the cell, for example by accumulation at the cell rear.

This model recalls a yeast cells polarisation model developed and studied in [19, 20, 24–26], where some polarity cues are actively transported towards the cell membrane by the cytoskeleton. In [27], we performed the theoretical study of a close model in dimension 1, with two main differences. First, the boundary pressure writes as an affine function of the concentration (e.g $p = 1 - c$). Moreover, the actin velocity depends only on the pressure gradient, so that the cell displacement does not induce direct feedback on the dynamics. For this 1D model, resting and moving stationary states coexist. Moreover, the system shows a dichotomy phenomenon: for a surcritical
molecular mass, a sufficiently polarised initial condition leads to the blowup of solutions in finite time. In the subcritical case, weak solutions exist globally, and converge towards the resting stationary state with explicit rate.

## 2.2 Some estimates for (3) – (6)

In this Subsection, we provide energy estimates for (3) – (6) for sufficiently smooth solutions. We first recall some results concerning the Dirichlet-to-Neumann operator, and then prove energy estimates using these properties.

### 2.2.1 The Dirichlet-to-Neumann operator

**Definition 1** Let $\Omega$ be a smooth enough open set of $\mathbb{R}^2$, the Dirichlet-to-Neumann operator $\mathcal{I}_{\Omega}$ is defined by

$$ \mathcal{I}_{\Omega}: H^1(\partial\Omega) \to L^2(\partial\Omega) \\ \varphi \to \nabla q \cdot \mathbf{n}, $$

where the function $q$ denotes the harmonic extension of $\varphi$ to the domain $\Omega$, that is

$$ \begin{aligned} \Delta q &= 0 \quad \text{in } \Omega, \\ q &= \varphi \quad \text{on } \partial\Omega. \end{aligned} $$

The name of this operator comes from the fact that it makes the transition between a problem verifying a Dirichlet boundary condition towards a problem with a Neumann boundary condition.

**Proposition 1** Let $\Omega$ and $\varphi$ be smooth enough and let $q$ be the harmonic extension of $\varphi$, then it holds that

$$ \int_{\partial\Omega} \varphi \mathcal{I}_{\Omega} \varphi \, d\sigma = \int_{\Omega} |\nabla q|^2 \, dx. $$

Hence

$$ \int_{\partial\Omega} \varphi \mathcal{I}_{\Omega} \varphi \, d\sigma = \inf \left\{ \int_{\Omega} |\nabla \psi|^2 \, dx, \forall \psi \in H^1(\Omega), \text{ s.t. } \psi = \varphi \text{ on } \partial\Omega \right\} $$

In the case where $\Omega$ is the disk of radius $R_0 > 0$, by expanding the function $\varphi$ in Fourier series, it is possible to define $\mathcal{I}_{R_0} := \mathcal{I}_{B_{R_0}}$ as follows:

**Definition 2** The Dirichlet-to-Neumann operator $\mathcal{I}_{R_0} : H^1(\mathbb{R}/2\pi\mathbb{Z}) \to L^2(\mathbb{R}/2\pi\mathbb{Z})$ is defined by

$$ \mathcal{I}_{R_0}(\varphi)(\theta) = \sum_{n=1}^{\infty} n R_0^{n-1} (a_n \cos(n\theta) + b_n \sin(n\theta)), $$

for $\varphi(\theta) = a_0 + \sum_{n=1}^{\infty} R_0^n (a_n \cos(n\theta) + b_n \sin(n\theta))$.
This is simply because harmonic functions $q$ on $B_{R_0}$ write in polar coordinates for $(r, \theta) \in [0, R] \times \mathbb{R}/2\pi\mathbb{Z}$ as:

$$q(r, \theta) = \sum_{n \ge 0} a_n \cos(n\theta)r^n + b_n \sin(n\theta)r^n,$$

when we discarded solutions that diverge at $r = 0$, and then

$$\nabla q(R_0, \theta) \cdot \mathbf{n} = \partial_r q(R_0, \theta) = \sum_{n \ge 1} n R_0^{n-1} (a_n \cos(n\theta) + b_n \sin(n\theta)).$$

We deduce the following Lemma.

**Lemma 2** The operator $\mathcal{I}_{R_0}$ has $R_0^{-1}\mathbb{N}^*$ for set of eigenvalues. Furthermore, Ker $(-\mathcal{I} + nI_d) = \text{span}\{\cos(n\theta)\}$ has dimension 1.

For brevity, in the sequel $\mathcal{I}_{R=1}$ will be simply denoted by $\mathcal{I}$.

### 2.2.2 Energy estimates

We now introduce an energy for which we provide a priori estimates.

**Proposition 3** Assume $c$ is a solution to (3) – (6) and $f$ is a smooth function on $\mathbb{R}_+$. Then, the following equality holds

$$\frac{d}{dt} \int_{B_1} f(c) \, dx = - \int_{\partial B_1} (f(c) - cf'(c)) \mathcal{I}(\chi k_p(c)) \, d\sigma - D \int_{B_1} f''(c) |\nabla c|^2 \, dx \\ - \gamma \chi \left( \int_{\partial B_1} k_p(c) \mathbf{n} \, d\sigma \right) \cdot \left( \int_{\partial B_1} (f(c) - cf'(c)) \mathbf{n} \, d\sigma \right) \\ + \frac{1}{2} \int_{B_1} \begin{pmatrix} x \\ y \end{pmatrix} \cdot \nabla (f(c) - cf'(c)) \, dx.$$

*Proof* We first see that

$$\frac{d}{dt} \int_{B_1} f(c) \, dx = \int_{B_1} f'(c) \operatorname{div} (D \nabla c - c \mathbf{u}) \, dx.$$

Then, we have that

$$\int_{B_1} f'(c) \operatorname{div}(c\mathbf{u}) \, dx = - \int_{B_1} cf'(c) \, dx - \int_{B_1} \nabla f(c) \nabla p \, dx \\ + \gamma \chi \left( \int_{\partial B_1} k_p(c) \mathbf{n} \, d\sigma \right) \left( \int_{B_1} \nabla f(c) \, dx \right) \\ = \int_{B_1} (f(c) - cf'(c)) \, dx - \int_{\partial B_1} f(c) \nabla p \cdot \mathbf{n} \, d\sigma \\ + \gamma \chi \left( \int_{\partial B_1} k_p(c) \mathbf{n} \, d\sigma \right) \cdot \left( \int_{\partial B_1} f(c) \mathbf{n} \, d\sigma \right).$$
On the other hand, using the boundary condition in (3), we have that

$$
\begin{align*}
D \int_{B_1} f'(c) \Delta c \, dx &= -D \int_{B_1} f''(c) |\nabla c|^2 \, dx + D \int_{\partial B_1} f'(c) \nabla c \cdot \mathbf{n} \, d\sigma \\
&= -D \int_{B_1} f''(c) |\nabla c|^2 \, dx - \int_{\partial B_1} c f'(c) \nabla p \cdot \mathbf{n} \, d\sigma \\
&\quad + \gamma \chi \left( \int_{\partial B_1} k_p(c) \mathbf{n} \, d\sigma \right) \cdot \left( \int_{\partial B_1} c f'(c) \mathbf{n} \, d\sigma \right),
\end{align*}
$$

which leads to

$$
\begin{align*}
\frac{d}{dt} \int_{B_1} f(c) \, dx = & -\gamma\chi \left( \int_{\partial B_1} k_p(c) \mathbf{n} \, d\sigma \right) \cdot \left( \int_{\partial B_1} (f(c) - cf'(c)) \mathbf{n} \, d\sigma \right) \\
& -D \int_{B_1} f''(c) |\nabla c|^2 \, dx + \int_{\partial B_1} (f(c) - cf'(c)) \nabla p \cdot \mathbf{n} \, d\sigma \\
& - \int_{B_1} (f(c) - cf'(c)) \, dx.
\end{align*}
$$

Moreover, defining $\bar{p}$ in $B_1$ such that $\bar{p}(t, \mathbf{x}) := p(t, \mathbf{x}) - \frac{1}{4} \| \mathbf{x} \|_2^2$, we obtain that for any smooth function $g$ on $R_+$,

$$
\begin{align*}
\int_{\partial B_1} g(c) \nabla p \cdot \mathbf{n} \, d\sigma &= \int_{\partial B_1} g(c) \nabla \bar{p} \cdot \mathbf{n} \, d\sigma + \frac{1}{2} \int_{\partial B_1} g(c) \begin{pmatrix} x \\ y \end{pmatrix} \cdot \mathbf{n} \, d\sigma \\
&= - \int_{\partial B_1} g(c) \mathcal{I} \left( \chi k_p(c) + \frac{1}{4} \right) d\sigma + \frac{1}{2} \int_{\partial B_1} g(c) \begin{pmatrix} x \\ y \end{pmatrix} \cdot \mathbf{n} \, d\sigma.
\end{align*}
$$

Consequently,

$$
\begin{align*}
\frac{d}{dt} \int_{B_1} f(c) \, dx &= -\gamma\chi \left( \int_{\partial B_1} k_p(c) \mathbf{n} \, d\sigma \right) \cdot \left( \int_{\partial B_1} (f(c) - cf'(c)) \mathbf{n} \, d\sigma \right) \\
&\quad -D \int_{B_1} f''(c) |\nabla c|^2 \, dx - \int_{\partial B_1} (f(c) - cf'(c)) \mathcal{I} \left( \chi k_p(c) + \frac{1}{4} \right) \, d\sigma \\
&\quad +\frac{1}{2} \int_{\partial B_1} (f(c) - cf'(c)) \begin{pmatrix} x \\ y \end{pmatrix} \cdot \mathbf{n} \, d\sigma - \int_{B_1} (f(c) - cf'(c)) \, dx, \\
&= -\gamma\chi \left( \int_{\partial B_1} k_p(c) \mathbf{n} \, d\sigma \right) \cdot \left( \int_{\partial B_1} (f(c) - cf'(c)) \mathbf{n} \, d\sigma \right) \\
&\quad -D \int_{B_1} f''(c) |\nabla c|^2 \, dx - \int_{\partial B_1} (f(c) - cf'(c)) \mathcal{I} \left( \chi k_p(c) \right) \, d\sigma \\
&\quad +\frac{1}{2} \int_{B_1} \nabla (f(c) - cf'(c)) \cdot \begin{pmatrix} x \\ y \end{pmatrix} \, dx,
\end{align*}
$$

by integration by part, and because $\mathcal{I}(\frac{1}{4}) = 0$ from Definition 2. $\square$

Next, considering a particular function $f$ yields the following estimation.

**Proposition 4** Assume that

$$
f(c) - cf'(c) = -\chi k_p(c). \qquad (8)
$$

Then,
1. Any $H^1$ solution to (3) – (6) satisfies

$$ \frac{d}{dt} \int_{B_1} f(c) \, dx \le - \int_{B_1} \left( D f''(c) - \chi^2 \left(1 + \frac{1}{4} + \gamma\pi\right) |k'_p(c)|^2 \right) |\nabla c|^2 \, dx \\ + \frac{\pi}{4}. \qquad (9) $$

2. If it holds additionally that

$$ D f''(c) \geq \chi^2 \left(1 + \frac{1}{4} + \gamma\pi\right) |k'_p(c)|^2, \qquad (10) $$

then

$$ \frac{d}{dt} \int_{B_1} f(c) \, dx \le \frac{\pi}{4}. $$

*Proof* Consider a function $f$ such that Equation (8) holds. Then, Proposition 3 rewrites

$$ \frac{d}{dt} \int_{B_1} f(c) \, dx = \gamma \chi^2 \left( \int_{\partial B_1} k_p(c) \mathbf{n} \, d\sigma \right)^2 - D \int_{B_1} f''(c) |\nabla c|^2 \, dx \\ + \chi^2 \int_{\partial B_1} k_p(c) \mathcal{I}(k_p(c)) \, d\sigma - \frac{\chi}{2} \int_{B_1} \nabla k_p(c) \cdot \begin{pmatrix} x \\ y \end{pmatrix} \, dx, $$

From Proposition 1, we obtain that

$$ \int_{\partial B_1} k_p(c) \mathcal{I}(k_p(c)) \, d\sigma \le \int_{B_1} |\nabla k_p(c)|^2 \, dx \\ \le \int_{B_1} |k'_p(c)|^2 |\nabla c|^2 \, dx. $$

The Cauchy-Schwarz inequality yields

$$ \left( \int_{\partial B_1} k_p(c) \mathbf{n} \, d\sigma \right)^2 = \left( \int_{B_1} \nabla k_p(c) \, dx \right)^2 \\ \le \pi \int_{B_1} |k'_p(c)|^2 |\nabla c|^2 \, dx. $$

Finally, using Young's inequality we see that

$$ \left| \frac{\chi}{2} \int_{B_1} \nabla k_p(c) \cdot \begin{pmatrix} x \\ y \end{pmatrix} dx \right| \le \frac{\pi}{8} + \frac{\chi^2}{4} \int_{B_1} |k'_p(c)|^2 |\nabla c|^2 dx. $$

Hence, the Inequality (9) follows. The point 2 easily follows provided that (10) is true. $\square$

Finally, we consider the particular case where $k_p$ is given by Equation (7) and we define the critical value $\chi_c$ by

$$ \chi_c := \frac{4D}{c_t \left(1 + \frac{1}{4} + \gamma\pi\right)}. \qquad (11) $$
**Proposition 5** Assume $k_p$ is given by (7) and $\chi \le \chi_c$. Let $c$ belongs to $H^1 \cap L \log L$ be a solution to (3) – (6), then

$$ \frac{d}{dt} \int_{B_1} c \log \left( \frac{c}{1+c} \right) \mathrm{d}\mathbf{x} \le \frac{\pi}{8}. $$

*Proof* In such a case, let $f$ be such that Equation (8) is verified. Then, for $x > 0$, one has

$$ f''(x) = \chi \frac{k'_p(x)}{x} = \chi \frac{c_t}{x(c_t + x)^2}, $$

and it follows that

$$ f(x) = A_1 + A_2 x + \chi c_t^{-1} x \log \left( \frac{x}{c_t + x} \right), $$

for $A_1, A_2 \in \mathbb{R}$. Then, Condition (10) rewrites

$$ \chi \left(1 + \frac{1}{4} + \gamma\pi\right) c_t \frac{x}{(1+x)^2} \le D, \quad \forall x > 0. \quad (12) $$

Since $\frac{x}{(1+x)^2} \le \frac{1}{4}$, for all $x \ge 0$, a sufficient condition for (12) to be true is $\chi \le \chi_c$, hence the result. $\square$

### 3 Analysis of the stationary problem associated with (3) – (6)

In this Section, we study the stationary problem associated with (3) – (6) in $\Omega = B(0, 1)$, which reads

$$ \left\{ \begin{array}{ll} \operatorname{div} (c\mathbf{u} - D\nabla c) = 0 & \text{in } \Omega, \\ (D\nabla c - c\mathbf{u}) \cdot \mathbf{n} = 0 & \text{on } \partial\Omega, \end{array} \right. \quad (13) $$

with

$$ \mathbf{u}(\mathbf{x}) = -\nabla p(\mathbf{x}) - \mathbf{v} \quad \text{on } \Omega, \quad (14) $$

and

$$ \left\{ \begin{array}{ll} \Delta p = 1 & \text{in } \Omega, \\ p = -\chi k_p(c) & \text{on } \partial\Omega. \end{array} \right. \quad (15) $$

The stationary cell velocity is

$$ \mathbf{v} = \gamma \int_{\partial\Omega} p \mathbf{n} \mathrm{d}\sigma = -\gamma\chi \int_{\partial\Omega} k_p(c) \mathbf{n} \mathrm{d}\sigma. \quad (16) $$

In Subsection 3.1, we prove two preliminary lemmas. In Subsection 3.2, we investigate the existence and uniqueness of a resting stationary state associated with a null cell velocity ($\mathbf{v} = 0$). In Subsection 3.3, we provide conditions under which this stationary solution is linearly stable. Finally, in Subsection 3.4, we prove the existence of a one-parameter family of traveling stationary solutions.
## 3.1 Preliminary lemmas

The first Lemma gives the solution of an auxiliary Laplace-Dirichlet problem involved in the pressure. For a given function $c$, the function defined by $p_0(x, y) := p(x, y) - \frac{1}{4}(x^2 + y^2)$ with $p$ solution of Problem (15), is such that

$$ \begin{cases} \Delta p_0(x, y) = 0, & \text{in } \Omega \\ p_0 = -\frac{1}{4} - \chi k_p(c), & \text{on } \partial\Omega, \end{cases} \qquad (17) $$

which rewrites in polar coordinates, keeping the same notations $c(r, \theta)$ and $p_0(r, \theta)$, as

$$ \begin{cases} -\partial_r (r \partial_r p_0) - \frac{1}{r} \partial_\theta^2 p_0 = 0, & \text{in } (0, 1) \times (0, 2\pi), \\ p_0(1, \theta) = -\frac{1}{4} - \chi k_p(c(1, \theta)) & \text{on } (0, 2\pi). \end{cases} \qquad (18) $$

We have the following classical lemma:

**Lemma 6** The unique solution of (18) in $B(0, R)$, $R > 0$, writes in integral form as

$$ p_0(r, \theta) = -\frac{1}{2\pi R} \int_0^{2\pi R} \left( \frac{1}{4} + \chi k_p(c(R, s)) \right) \frac{R^2 - r^2}{R^2 - 2rR \cos\left(\frac{s}{R} - \theta\right) + r^2} \mathrm{d}s, \qquad (19) $$

and its Fourier expansion writes

$$ \begin{align} p_0(r, \theta) &= \frac{a_0}{2} + \sum_{n \ge 1} r^n (a_n \cos(n\theta) + b_n \sin(n\theta)), \quad \text{with} \nonumber \\ a_0 &= -\frac{1}{\pi} \int_{-\pi}^{\pi} \left( \frac{1}{4} + \chi k_p(c(R, \theta')) \right) \mathrm{d}\theta', \nonumber \\ a_n &= -\frac{1}{\pi} \int_{-\pi}^{\pi} \left( \frac{1}{4} + \chi k_p(c(R, \theta')) \right) \cos(n\theta') \mathrm{d}\theta', \quad n \ge 1, \tag{20} \\ b_n &= -\frac{1}{\pi} \int_{-\pi}^{\pi} \left( \frac{1}{4} + \chi k_p(c(R, \theta')) \right) \sin(n\theta') \mathrm{d}\theta', \quad n \ge 1. \nonumber \end{align} $$

The second lemma gives the general expression of the solutions $c$ of (13) – (14) for given pressure $p$ and velocity $v$.

**Lemma 7** For a given function $p$, and a given $v \in \mathbb{R}$, any nonnegative solution of (13) – (14) is given by

$$ c(x, y) = c_0 e^{-\frac{p(x, y)+v \cdot (x, y)}{D}}, \quad (x, y) \in \Omega, \qquad (21) $$

where $c_0 > 0$.

*Proof* Indeed, it is straightforward that $c(x, y) = e^{-\frac{p(x, y)+v \cdot (x, y)}{D}}$ is a solution of (13) – (14). Let us prove its uniqueness up to a multiplicative constant. Indeed, for an unknown function $\tilde{c}$, write

$$ c(x, y) = \tilde{c}(x, y) e^{-\frac{p(x, y)+v \cdot (x, y)}{D}}. $$
Then, $\tilde{c}$ satisfies a Laplace equation with Neumann boundary conditions given by

$$ \begin{cases} \Delta \tilde{c} = 0 & \text{in } \Omega, \\ \nabla \tilde{c} \cdot \mathbf{n} = 0 & \text{on } \partial \Omega. \end{cases} \qquad (22) $$

Since $\Omega$ is smooth enough (e.g. $\partial\Omega$ is of class $C^2$), the maximum principle yields that $\tilde{c} \equiv c_0$ constant in $\Omega$. $\square$

## 3.2 Existence and uniqueness of the rest state associated with (3)-(6)

Let us define

$$ c^* := \frac{1}{4\pi D \left(e^{\frac{1}{4D}} - 1\right)}, \qquad (23) $$

and

$$ \chi_0 := \frac{D}{A}, \qquad (24) $$

where $A$ is given by (A3). We now prove that there exists a unique resting stationary state for (3) – (6).

**Proposition 8** (Radially symmetric stationary solution) Assume that $k_p$ satisfies Assumption 1. Then, for all $\gamma, D, \chi > 0$, there exists a stationary solution to (3) – (6) such that $\mathbf{v} = 0$. It is given in polar coordinates for $r \in [0, 1]$ by

$$ \bar{c}(r) = c^* e^{\frac{1}{4D}(1-r^2)}. \qquad (25) $$

If additionally $k_p$ satisfies Assumption 1 and $\chi \in (0, \chi_0)$, it is the unique resting stationary solution.

*Proof* The proof is divided in two steps. First, we show that for a given boundary value for $c$, there exists a unique axisymmetric solution $c$. Then, we consider the fully coupled system and prove uniqueness of the solution.

Let us assume that $c$ is axisymmetric, so that $c|_{\partial\Omega} = c^*$ for some $c^* > 0$. Using Lemma 6, we deduce from Equation (20) that only $a_0$ is nonzero, so that

$$ p(r, \theta) = \frac{1}{4}(r^2 - 1) - \chi k_p(c^*). \qquad (26) $$

By the axisymmetric assumption on $c$, Equation (16) yields that $\mathbf{v} = 0$, hence

$$ u(r, \theta) = -\partial_r p e_r - \frac{1}{r} \partial_\theta p e_\theta = -\frac{1}{2} r e_r. \qquad (27) $$

The polar problem associated with (13) writes

$$ \begin{cases} \partial_r (r(cu_r - D\partial_r c)) + \partial_\theta (cu_\theta - D\partial_\theta c) = 0, & \text{for } r < 1, \\ (cu_r - D\partial_r c) = 0, & \text{for } r = 1. \end{cases} \qquad (28) $$

Since $u_\theta = 0$ and the boundary value $c^*$ is constant, it is clear that the solution is axisymmetric ($\partial_\theta c = 0$). We write now $c(r)$ for the unknown, and the corresponding problem rewrites

$$ \begin{cases} \partial_r (\frac{1}{2} r^2 c + D r \partial_r c) = 0, & \text{for } r < 1, \\ \frac{1}{2} r c + D \partial_r c = 0, & \text{for } r = 1, \\ c = c^*, & \text{for } r = 1. \end{cases} \qquad (29) $$
The first equation can be written

$$Dr \partial_r^2 c + \left( D + \frac{1}{2} r^2 \right) \partial_r c + rc = 0. \qquad (30)$$

For $0 < r \le 1$, the solutions of this equation form a bidimensional space. One solution is given by

$$c_1(r) = c^* e^{\frac{1}{4D}(1-r^2)}, \text{ with } c^* = \frac{1}{4\pi D \left(e^{\frac{1}{4D}} - 1\right)},$$

since

$$\int_{0}^{1} c^{*} e^{\frac{1}{4D}(1-r^{2})} r \, dr = \frac{1}{2\pi}.$$

The Wronskien associated with the second order ODE is computed to find another solution $c_2$ such that $\{c_1, c_2\}$ forms a basis of the space of solutions. The Wronskien $W$ is solution of

$$W'(r) = -\frac{D + \frac{1}{2}r^2}{Dr} W(r) = -\left(\frac{1}{r} + \frac{1}{2D}r\right) W(r). \qquad (31)$$

The solution writes

$$W(r) = W(1) \frac{1}{r} e^{\frac{1}{4D}(1-r^2)}. \qquad (32)$$

Moreover, we have $W(r) = c_1(r)c'_2(r) - c_2(r)c'_1(r)$, leading to the following ODE on $c_2$:

$$c'_{2}(r) + \frac{1}{2D}rc_{2}(r) = \frac{W(1)}{c^{*}_{r}}, \qquad (33)$$

A simple computation shows that $c_2$ satisfies the boundary condition in Problem (29) if and only if $W(1) = 0$, leading to $c_2 = c_1$. Therefore, for any given boundary value $c^* > 0$, the unique axisymmetric stationary solution is given by Equation (25).

Now consider the fully coupled problem. Lemma 7 provides the general expression of $c$, so that for $v = 0$, we have the following Poisson-Dirichlet problem:

$$\begin{cases} \Delta p = 1 & \text{in } \Omega, \\ p = -\chi k_p \left(c_0 e^{-\frac{p}{D}}\right) & \text{on } \partial\Omega, \end{cases} \qquad (34)$$

where $c_0$ is a positive constant. Uniqueness of a solution follows from uniqueness of a pressure value satisfying the boundary condition. Let us define the function $g$ in $\mathbb{R}$ by

$$g(x) := x + \chi k_p \left(c_0 e^{-\frac{x}{D}}\right),$$

so that

$$g'(x) = 1 - \frac{\chi}{D} c_0 e^{-\frac{x}{D}} k'_p \left(c_0 e^{-\frac{x}{D}}\right).$$

Therefore, from Assumption 1 together with $\chi < \chi_0 = \frac{D}{A}$, we obtain that $g$ is increasing. Hence, there exists a unique negative root of $g$. This proves the uniqueness of the solution to (34).

□

In the case where $k_p$ is given by Equation (7), the following corollary follows.

**Corollary 9** Assume that the function $k_p$ is given by (7) and that $\chi < \chi_0 = 4D$. Then, there exists a unique stationary solution to (3) – (6) such that $v = 0$ and it is given by Equation (25).

*Proof* The proof follows from the inequality $\frac{x}{(c_t+x)^2} \le \frac{1}{4}$ for all $x \ge 0$. □
### 3.3 Linear stability analysis of the radially symmetric stationary solution for (3) – (6)

Let us consider a radially symmetric reference stationary solution $(\bar{c}, \bar{p})$ where $\bar{c}$ is given by Equation (25) and $\bar{p}$ is solution of

$$ \begin{cases} -\Delta \bar{p} = -1 & \text{in } \Omega, \\ \bar{p} = -\chi k_p(\bar{c}) & \text{on } \partial \Omega. \end{cases} $$

Then, linearizing Problem (3) – (6) around $(\bar{c}, \bar{p})$, and writing

$$ c = \bar{c} + \varepsilon \tilde{c} + \text{h.o.t} , \quad p = \bar{p} + \varepsilon \tilde{p} + \text{h.o.t} , $$

one can get the following problem for $\tilde{c}$ and $\tilde{p}$:

$$ \left\{ \begin{array}{ll} \partial_t \tilde{c} = \text{div} \left( D \nabla \tilde{c} + \bar{c} \nabla \bar{p} + \bar{c} \nabla \tilde{p} + \gamma \bar{c} \int_{\partial \Omega} \tilde{p} \mathbf{n} \, d\sigma \right) & \text{in } \Omega, \\ (D \nabla \tilde{c} + \bar{c} \nabla \bar{p} + \bar{c} \nabla \tilde{p} + \gamma \bar{c} \int_{\partial \Omega} \tilde{p} \mathbf{n} \, d\sigma) \cdot \mathbf{n} = 0 & \text{on } \partial \Omega, \end{array} \right. \quad (35) $$

and

$$ \begin{cases} -\Delta \tilde{p} = 0 & \text{in } \Omega, \\ \tilde{p} = -\chi k'_p(\bar{c}) \tilde{c} & \text{on } \partial \Omega. \end{cases} $$

Some estimates on $\tilde{c}$ and $\tilde{p}$ can be obtained. First, the mass constraint on $c$ and $\bar{c}$ yields $\int_{\Omega} \tilde{c} \, d\mathbf{x} = 0$, and the Poincaré inequality gives

$$ \int_{\Omega} \tilde{c}^2 \, d\mathbf{x} \le \left(\frac{2}{\pi}\right)^2 \int_{\Omega} |\nabla \tilde{c}|^2 \, d\mathbf{x}. $$

Furthermore, since the harmonic extension minimizes the Dirichlet energy for a fixed boundary data, we have

$$ \int_{\Omega} |\nabla \tilde{p}|^2 \, d\mathbf{x} \le \chi^2 (k'_p(c^*))^2 \int_{\Omega} |\nabla \tilde{c}|^2 \, d\mathbf{x}. \quad (36) $$

Finally, let us define

$$ \chi_1 := \frac{4\pi D(1 - e^{-\frac{1}{4D}})(D - \frac{1}{\pi})}{(1 + \gamma\pi) |k'_p(c^*)|} . \quad (37) $$

We can prove the following Proposition.

**Proposition 10 (Linear stability)** Assume that $k_p$ satisfies assumptions 1-1 and that $\chi < \chi_1$. Then, the radially symmetric reference stationary solution $\bar{c}$ given by (23)–(25) is linearly stable.
Proof We first see that

$$
\begin{align*}
\frac{d}{dt} \int_{\Omega} |\tilde{c}|^2 \, \mathrm{d}\mathbf{x} &= \int_{\Omega} \tilde{c} \operatorname{div} \left( D \nabla \tilde{c} + \tilde{c} \nabla \tilde{p} + \tilde{c} \nabla \tilde{p} + \gamma \tilde{c}(\mathbf{x}) \int_{\partial \Omega} \tilde{p}(t, \mathbf{x}) \, \mathrm{d}\sigma \right) \, \mathrm{d}\mathbf{x} \\
&= -D \int_{\Omega} |\nabla \tilde{c}|^2 \, \mathrm{d}\mathbf{x} - \int_{\Omega} \tilde{c} \nabla \tilde{c} \cdot \nabla \tilde{p} \, \mathrm{d}\mathbf{x} - \int_{\Omega} \tilde{c} \nabla \tilde{c} \cdot \nabla \tilde{p} \, \mathrm{d}\mathbf{x} \\
&\quad - \gamma \int_{\Omega} \tilde{c}(\mathbf{x}) \nabla \tilde{c}(t, \mathbf{x}) \, \mathrm{d}\mathbf{x} \int_{\partial \Omega} \tilde{p}(t, \mathbf{x}) \, \mathrm{d}\sigma.
\end{align*}
$$

Using the Poincaré inequality, we obtain

$$
\begin{aligned}
\left| \int_{\Omega} \tilde{c} \nabla \tilde{c} \cdot \nabla \tilde{p} \, \mathrm{d}\mathbf{x} \right| &\le \|\nabla \tilde{p}\|_{L^{\infty}(\Omega)} \left( \int_{\Omega} |\nabla \tilde{c}|^2 \, \mathrm{d}\mathbf{x} \right)^{1/2} \left( \int_{\Omega} |\tilde{c}|^2 \, \mathrm{d}\mathbf{x} \right)^{1/2} \\
&\le \frac{2}{\pi} \|\nabla \tilde{p}\|_{L^{\infty}(\Omega)} \int_{\Omega} |\nabla \tilde{c}|^2 \, \mathrm{d}\mathbf{x}.
\end{aligned}
$$

Furthermore, using (36), we have

$$
\begin{aligned}
\left| \int_{\Omega} \bar{c} \nabla \bar{c} \cdot \nabla \bar{p} \, \mathrm{d}\mathbf{x} \right| &\le \|\bar{c}\|_{L^{\infty}(\Omega)} \int_{\Omega} |\nabla \bar{c} \cdot \nabla \bar{p}| \, \mathrm{d}\mathbf{x} \\
&\le \|\bar{c}\|_{L^{\infty}(\Omega)} \left( \int_{\Omega} |\nabla \bar{c}|^2 \, \mathrm{d}\mathbf{x} \right)^{1/2} \left( \int_{\Omega} |\nabla \bar{p}|^2 \, \mathrm{d}\mathbf{x} \right)^{1/2} \\
&\le \chi |k'_{p}(c^*)| \|\bar{c}\|_{L^{\infty}(\Omega)} \int_{\Omega} |\nabla \bar{c}|^2 \, \mathrm{d}\mathbf{x}.
\end{aligned}
$$

Finally, since

$$
\begin{align*}
-\gamma \int_{\Omega} \bar{c} \nabla \bar{c} \, \mathrm{d}\mathbf{x} \int_{\partial \Omega} \bar{p} \, \mathrm{d}\sigma &= \gamma \chi k'_{p}(c^*) \int_{\Omega} \bar{c} \nabla \bar{c} \, \mathrm{d}\mathbf{x} \int_{\partial \Omega} \bar{c} \, \mathrm{d}\sigma \\
&= \gamma \chi k'_{p}(c^*) \int_{\Omega} \bar{c} \nabla \bar{c} \, \mathrm{d}\mathbf{x} \int_{\Omega} \nabla \bar{c} \, \mathrm{d}\mathbf{x},
\end{align*}
$$

and

$$
\gamma \chi k'_{p}(c^{*}) \int_{\Omega} \bar{c} \nabla \bar{c} \, \mathrm{d}\mathbf{x} \leq \gamma \chi |k'_{p}(c^{*})| \|\bar{c}\|_{L^{\infty}(\Omega)} \int_{\Omega} |\nabla \bar{c}| \, \mathrm{d}\mathbf{x},
$$

it holds that

$$
-\gamma \int_{\Omega} \bar{c} \nabla \bar{c} \, \mathrm{d}\mathbf{x} \int_{\partial \Omega} \bar{p} \, \mathrm{d}\sigma \leq \gamma \chi \pi |k'_{p}(c^*)| \|\bar{c}\|_{L^{\infty}(\Omega)} \int_{\Omega} |\nabla \bar{c}|^2 \, \mathrm{d}\mathbf{x}.
$$

Consequently, we have

$$
\begin{align*}
& \frac{d}{dt} \int_{\Omega} |\tilde{c}|^2 \, \mathrm{d}\mathbf{x} \\
&\le - \left( D - \chi (1 + \gamma\pi) |k'_{p}(c^*)| \|\tilde{c}\|_{L^{\infty}(\Omega)} - \frac{2}{\pi} \|\nabla \tilde{p}\|_{L^{\infty}(\Omega)} \right) \int_{\Omega} |\nabla \tilde{c}|^2 \, \mathrm{d}\mathbf{x},
\end{align*}
$$

and using Equations (25) and (26), we deduce that $\|\bar{c}\|_{L^{\infty}(\Omega)} = c^* e^{\frac{1}{4D}}$, together with $\|\nabla \bar{p}\|_{L^{\infty}(\Omega)} = \frac{1}{2}$. Finally, Assumption 1 ensures that $|k'_{p}(c^*)| > 0$, from which the sufficient condition follows. $\square$

In the particular case of Equation (7), we obtain a sufficient condition for the linear stability of the rest state. Let us define

$$
\bar{\chi}_1 := \frac{4(\pi D - 1)e^{\frac{1}{4D}}}{\pi (1 + \gamma\pi)}. \qquad (38)
$$
**Corollary 11** In the case where $k_p$ is given by Equation (7), a sufficient condition for the linear stability of the rest state is that $\chi < \bar{\chi}_1$.

*Proof* The result follows from Assumption 1 that yields $|k'_p(c^*)| \le \frac{1}{4c^*}$. □

## 3.4 Existence of travelling wave stationary solutions, $v \ne 0$

Without loss of generality, we assume forward motion in the $x$-direction and write $v = (V, 0)^T$ with $V \in \mathbb{R}_+$. Lemma 7 allows us to eliminate $c$ from (13)–(14)–(15)–(16) and to rewrite the stationary problem with travelling wave in the following concise form:

$$ \begin{cases} -\Delta p = -1 & \text{in } \Omega_0, \\ p(1, \theta) = -\chi k_p \left( c_0 e^{-\frac{p(1, \theta)+V \cos(\theta)}{D}} \right) & \text{on } \partial \Omega, \end{cases} \quad (39) $$

together with

$$ V = -\gamma\chi \int_{0}^{2\pi} k_p \left( c_0 e^{-\frac{p(1,\theta)+V\cos(\theta)}{D}} \right) \cos(\theta) \, d\theta, \quad (40) $$

where $\int_{\Omega} c(x, y) \, dx \, dy = 1$, leading to

$$ c_0 = \frac{1}{\int_{\Omega} e^{-\frac{p(x,y)+xV}{D}} \, dx \, dy} \quad (41) $$

Finally, let us define

$$ \gamma_c(\chi, D) := \pi^{-1} \left[ \left( \frac{\chi}{D} c^* e^{\frac{\chi}{D} k_p(c^*)} k'_p \left( c^* e^{\frac{\chi}{D} k_p(c^*)} \right) \right)^{-1} - 1 \right]. $$

We prove the following proposition.

**Proposition 12** (Traveling wave solutions) Assume that $k_p$ satisfies Assumptions 1-1-1, and that $\chi \in (0, \chi_0)$. Then, for all $D > 0$, for $c^*(D)$ given by Equation (23), and for $\gamma \ne \gamma_c(\chi, D)$, there exist traveling wave solutions of (3) – (6) such that $v = V e_x$ with $V > 0$.

*Proof* First, consider Problem (39) for given $V \ge 0$ and $(\chi, D) \in \mathbb{R}_+^2$. For all $\theta$, consider the function

$$ g_{\theta}(x) := x + \chi k_p \left( c_0 e^{-\frac{x+V \cos(\theta)}{D}} \right). $$

We see that

$$ g'_{\theta}(x) = 1 - \frac{\chi c_0}{D} e^{-\frac{x+V \cos(\theta)}{D}} k'_p \left( c_0 e^{-\frac{x+V \cos(\theta)}{D}} \right). $$

Then, the existence and uniqueness of the solution to (39) follows from Hypothesis 1 together with the assumption $\chi \in (0, \chi_0)$, with arguments similar to those given in the proof of Proposition 8.
The unique solution to Problem (39) depends on $(V, \chi, D)$ and the normalization constant $c_0$ given by (41), which depends on $(V, \chi, D)$. Therefore, we denote it by $p(r, \theta; V, \chi, D)$. Define the function $F$ by

$$F(V, \chi, D, \gamma) := V + \gamma\chi \int_{0}^{2\pi} k_{p} \left( c_{0}(V, \chi, D) e^{-\frac{p(1, \theta; V, \chi, D) + V \cos(\theta)}{D}} \right) \cos(\theta) d\theta.$$

For given parameters $(\chi, D, \gamma)$, the existence of a travelling wave solution, that is Equation (40), rewrites in the form

$$F(V, \chi, D, \gamma) = 0.$$

We apply the implicit function theorem. We first know that $F(0, \chi, D, \gamma) = 0$ for all $D, \gamma > 0$ and $\chi$ in $(0, \chi_0)$. Next, we compute

$$\partial_V F(V, \chi, D, \gamma) = 1 + \gamma\chi \int_0^{2\pi} \left( \partial_V c_0 - \frac{\partial_V p(1, \theta; V, \chi, D) + \cos(\theta)}{D} c_0 \right) \times e^{\frac{-p(1, \theta; V, \chi, D) - V \cos(\theta)}{D}} k'_p \left( c_0 e^{\frac{-p(1, \theta; V, \chi, D) - V \cos(\theta)}{D}} \right) \cos(\theta) d\theta,$$

where $q(1, \theta; V, \chi, D) := \partial_V p(1, \theta; V, \chi, D)$ satisfies the following Laplace problem

$$\left\{ \begin{array}{ll} -\Delta q = 0 & \text{in } \Omega, \\ q(1, \theta; V, \chi, D) = -\chi \left( \partial_V c_0(V, \chi, D) - \frac{q(1, \theta; V, \chi, D) + \cos(\theta)}{D} c_0(V, \chi, D) \right) \\ \times e^{\frac{-p(1, \theta; V, \chi, D) - V \cos(\theta)}{D}} k'_p \left( c_0(V, \chi, D) e^{\frac{-p(1, \theta; V, \chi, D) - V \cos(\theta)}{D}} \right) & \text{on } \partial\Omega, \end{array} \right.$$

hence $\tilde{q}(1, \theta; \chi, D) := q(1, \theta; V = 0, \chi, D)$ satisfies

$$\left\{ \begin{array}{ll} -\Delta \tilde{q} = 0 & \text{in } \Omega, \\ \tilde{q}(1, \theta; \chi, D) = -\chi \left( \partial_V c_0(0, \chi, D) - \frac{\tilde{q}(1, \theta; \chi, D) + \cos(\theta)}{D} c^* \right) \\ \times e^{\frac{-\bar{p}(1; \chi, D)}{D}} k'_p \left( c^* e^{\frac{-\bar{p}(1; \chi, D)}{D}} \right) & \text{on } \partial\Omega. \end{array} \right.$$

since $p(1, \theta; 0, \chi, D) = \bar{p}(1; \chi, D) = -\chi k_p(c^*)$ and $c_0(0, \chi, D) = c^*(D) > 0$. Consequently, we obtain

$$\tilde{q}(1, \theta; \chi, D) = -\partial_V c_0(V = 0, \chi, D) \frac{(1 - \frac{D}{c^*} B)}{B} + \frac{(1 - B)}{B} \cos(\theta),$$

where

$$B := 1 - \frac{\chi}{D} c^* e^{\frac{\chi}{D} k_p(c^*)} k'_p \left( c^* e^{\frac{\chi}{D} k_p(c^*)} \right) \in (0, 1),$$

for $\chi \in (0, \chi_0)$. Hence,

$$\partial_V F(V = 0, \chi, D, \gamma) = 1 - \gamma \pi \frac{1 - B}{B},$$

and $\partial_V F(V = 0, \chi, D, \gamma) \neq 0$ if and only if

$$\gamma \pi \neq \frac{B}{1 - B},$$

hence the result.

□

Assume that $k_p$ is given by Equation (7). Then,

$$\gamma_c(\chi, D) = \pi^{-1} \left[ \left( \frac{\chi}{D} \frac{c_t c^* e^{\frac{\chi}{D} k_p(c^*)}}{(c_t + c^* e^{\frac{\chi}{D} k_p(c^*)})^2} \right)^{-1} - 1 \right],$$

leading to the following Corollary.
**Corollary 13** If $k_p$ is given by Equation (7), then for all $\chi \in (0, \chi_0)$, $D > 0$, for $c^*(D)$ given by Equation (23), and for $\gamma \neq \gamma_c(\chi, D)$, there exist traveling wave solutions of (3) – (6) such that $\mathbf{v} = V\mathbf{e}_x$, with $V > 0$.

*Remark 1* It is of interest to identify the regime when $\gamma_c > 0$. Defining the quantity $E := \frac{\chi}{D}c^*e^{\frac{\chi}{D}k_p(c^*)}k'_p(c^*e^{\frac{\chi}{D}k_p(c^*)})$, it is clear that $\gamma_c > 0$ if and only if $E < 1$. When $k_p$ is given by Equation (7), using the variable $G := c^*e^{\frac{\chi}{D}k_p(c^*)}$, a discriminant analysis yields the following result: $\gamma_c > 0$ is equivalent to

* either $\frac{\chi}{D} < 4$,

* either $\frac{\chi}{D} = 4$ and $c^*(D) \neq c_0^*(c_t) > 0$ unique solution of $xe^{4k_p(x)} = c_t$,

* either $\frac{\chi}{D} > 4$ and $c^*(D) > c_1^*(\frac{\chi}{D}, c_t)$ with $c_1^*(\frac{\chi}{D}, c_t) > 0$ unique solution of $xe^{\frac{\chi}{D}k_p(x)} = \frac{c_t}{2}(\frac{\chi}{D} - 2 + \sqrt{\frac{\chi}{D}(\frac{\chi}{D} - 4)})$,

* or $\frac{\chi}{D} > 4$ and $c^*(D) < c_2^*(\frac{\chi}{D}, c_t)$ with $c_2^*(\frac{\chi}{D}, c_t) > 0$ unique solution of $xe^{\frac{\chi}{D}k_p(x)} = \frac{c_t}{2}(\frac{\chi}{D} - 2 - \sqrt{\frac{\chi}{D}(\frac{\chi}{D} - 4)})$.

# 4 Numerical simulations of the intracellular dynamics

This Section is devoted to numerical simulations of solutions to (3)–(4)–(5)–(6), modelling single-cell intracellular dynamics when it moves on a surface without any external stimulation. The activation function $k_p$ is defined by Equation (7), so that the results from Corollary 9, Corollary 11 and Corollary 13 hold.

In Section 3, a gaussian stationary state (23)–(25) was introduced. Its uniqueness among resting states was shown for $\chi < \chi_0$, as well as its linear stability for $\chi < \chi_1$. Moreover, for $\chi < \chi_0$ and $\gamma \neq \gamma_c$, the existence of traveling wave stationary solutions was established. Here, we both illustrate the cases treated by the theoretical analysis and investigate the other settings.

This Section is organized as follows. Subsection 4.1 gives a brief account of the numerical setting. In Subsection 4.2, the case of a radially symmetric stationary state is illustrated in relation with Corollary 9 and Corollary 11. In Subsection 4.3, Corollary 13 is treated by showing convergence towards a traveling stationary state. Finally, Subsection 4.4 is devoted to the parameter spaces.

## 4.1 Numerical setting

The numerical resolution is performed in polar coordinates. Given the molecular concentration at the boundary, Lemma 6 provides an explicit expression for the pressure. Therefore, the discretization focuses on the advection-diffusion equation. The advection field being nonlinear, the advection is treated explicitly, while the diffusion term is treated implicitly. We use a finite volumes scheme in space. Unless indicated otherwise, parameters of the numerical setting are $\Delta t = 0.01$ and $\Delta r = \Delta \theta = \frac{2\pi}{200} \approx 0.03$. The numerical discretization is further detailed in Appendix C.1.

For $k_p$ given by Equation (7), the parameters under study are the molecular diffusion coefficient $D$, the magnitude $\chi$ of the molecular regulation of the boundary pressure, the parameter $c_t > 0$ such that $k_p(c_t) = \frac{1}{2}$, and $\gamma$ the cell velocity parameter.
Note that the latter arises only as a product with $\chi$ (see Equation (6)). The molecular mass is fixed with $M = 1$.

The simulations are pursued up until the system reaches a numerical stationary state, defined as the time $t^N$ above which the $\ell^1$ distance between the molecular concentration at two consecutive time steps remains lower than $10^{-10}$ during one numerical time unit (or equivalently $\Delta t^{-1}$ time steps). Note that this criterion depends on the spatial discretization. For the parameter spaces, to reduce computational costs, the threshold is set at $10^{-6}$ and the space step is increased such that $\Delta\theta = \frac{2\pi}{120} \approx 0.05$.

the threshold is taken equal to $10^{-6}$ to lighten computational costs.

Two typical quantities can be used as outputs to characterize cell polarisation. The Relative Standard Deviation (RSD) is given by the perimeter-normalized ratio between the standard deviation of boundary values and their mean value. The Resultant Vector Length (RVL) is given by the norm of the normalized discretization of $\int_{\partial\Omega} c(1, \theta)\mathbf{n} d\theta$. Their formulas and a short discussion on their respective meanings can be found in Appendix C.2.

In the following, parameters are chosen as follows. For simulations illustrating stationary states (Figures 1–2–3), parameters $D$ and $\chi$ vary, while $(\gamma, c_t) = (0.5, 1)$. For parameter spaces, $\gamma = 0.1$ is fixed, while the other parameters, when they do not vary, are such that $(\chi, c_t, D) = (1, 0.4, 0.3)$.

## 4.2 Resting stationary state

This Subsection aims at illustrating both the relaxation towards the radially symmetric resting stationary state given by Equation (25), and its linear stability. First, let us consider the case of an initially polarized cell, with

$$c(0, r, \theta) = (1 - \sin(\theta)\mathbf{1}_{r>0})\mathcal{N}, \qquad (42)$$

where $\mathcal{N}$ is the normalization constant. Figure 1 illustrates the corresponding system for $D = 0.3$ and $\chi = 0.3 < \chi_0 = 1.2$. Subfigure (a) shows the initially polarized molecular concentration given by (42). Subfigure (b) shows the corresponding pressure field. The black arrows represent the associated radial and angular advection fields ($\mathbf{u}_r, \mathbf{u}_\theta$), while the red arrow shows the global cell velocity. The molecular asymmetry hence drives the pressure asymmetry from the boundary. Subfigure (c) shows the radially symmetric stationary molecular concentration, with a diffuse accumulation of molecules at the cell center. Subfigure (d) shows the time evolution of both polarity indices that rapidly decrease to zero. We could not find other resting stationary profiles for $\chi > \chi_0$.

Next, let us illustrate the linear stability of such stationary state. For that purpose, the initial condition is taken as a perturbed radially symmetric profile given by

$$c(0, r, \theta) = \bar{c}(r) + 0.1 \sin(\theta)\mathbf{1}_{r>0}, \qquad (43)$$

with $\bar{c}$ given by Equations (23)–(25). The Figure 2 illustrates the dynamics for $D = 0.8$ and several values for $\chi$. In this case, the respective threshold values involved in Proposition 10 and Corollary 11 are $\chi_1 \approx 0.82$ and $\bar{\chi}_1 \approx 1.02$. Subfigure (a) shows the initial
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Illustration of the relaxation of the system to a radially-symmetric stationary state, with $D = \chi = 0.3$ ($\chi < \chi_0$), $c_t = 1$ and $\gamma = 0.5$. (a): Polarized initial concentration given by Equation (42). (b): Initial pressure field given by Lemma 6, together with the corresponding radial (left) and angular (right) advection fields and the global cell velocity (red). (c): Stationary radially-symmetric molecular concentration. (d): Time evolution of the polarity indices (RSD in green, RVL in yellow), both curves overlap</figcaption></figure>

molecular concentration given by Equation (43). Subfigure (b) shows the stationary radial profile obtained for $\chi = 0.3 < \chi_1 < \bar{\chi}_1$ together with its analytical counterpart. Subfigure (c) shows the radially-symmetric stationary concentration obtained for $\chi = 1.85 > \bar{\chi}_1 > \chi_1$. Subfigure (d) shows the polarized stationary concentration profile obtained for $\chi = 1.9 > \bar{\chi}_1 > \chi_1$. These simulations show that the obtained theoretical threshold are not optimal

## 4.3 Traveling stationary states

This Subsection aims at illustrating the numerical convergence of the system towards a traveling stationary state. The initial concentration is given by Equation (42). Figure 3 illustrates the corresponding system for $D = 0.3$ and $\chi = 1 < \chi_0 = 1.2$, such that $\gamma \neq \gamma_c$. Subfigure (a) shows the initially polarized molecular concentration given by (42). Subfigure (b) shows the corresponding pressure field. Subfigure (c) shows the polarized stationary molecular concentration, with accumulation of molecules at a part of the boundary, thus defining the cell rear. Subfigure (d) shows the time evolution of
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Illustration of the relaxation to a radially-symmetric or polarized stationary state from a perturbed initial condition, with $D = 0.8$, $c_t = 1$ and $\gamma = 0.5$. (a): Perturbed initial concentration given by Equation (43). (b): Radial molecular profile at stationary state for $\chi = 0.3 < \chi_1 < \bar{\chi}_1$, together with the corresponding explicit solution. (c): Radially-symmetric stationary concentration obtained for $\chi = 1.85 > \bar{\chi}_1 > \chi_1$. (d): Polarized stationary concentration obtained for $\chi = 1.85 > \bar{\chi}_1 > \chi_1$.</figcaption></figure>

both polarity indices that grow towards a stationary value. We did not find significantly different profiles for either $\gamma = \gamma_c$ or $\chi > \chi_0$.

## 4.4 Parameter spaces

This Subsection is devoted to showing parameter spaces for both the Relative Standard Deviation and the Resultant Vector Length at stationary state, depending on the parameters $\chi$ and either $D$ or $ct$. The initial concentration is given by Equation (42), and the threshold for numerical stationary state is taken equal to $10^{-6}$. The results are given in Figure 4. Subfigures (a)–(b) treat the pair $(c_t, \chi)$ and show both unpolarized and polarized stationary states. One can see that below a threshold value on $\chi$ the concentration converges towards the radially-symmetric profile, and that this threshold increases slowly with $c_t$. Subfigures (c)–(d) treat the pair $(D, \chi)$, and also presents both unpolarized and polarized stationary states. It is observed that both polarity indices decrease as the diffusion coefficient increases, as expected. Moreover, an affine relation between both parameters seems to rule the transition between both types of stationary states. Corollary 13 stated the existence of traveling stationary solutions for $\chi < \chi_0 = 4D$. One can see that it is also the case for $\chi \ge 4D$. All in all, these
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Illustration of the convergence of the system towards a traveling stationary state, with $D = 0.3$, $\chi = 1 < \chi_0 = 1.2$, $c_t = 1$ and $\gamma = 0.5$. (a): Polarized initial concentration given by Equation (42). (b): Initial pressure field given by Lemma 6, together with the corresponding radial (left) and angular (right) advection fields and the global cell velocity (red). (c): Stationary non radially-symmetric molecular concentration. (d): Time evolution of the polarity indices (RSD in green, RVL in yellow)</figcaption></figure>

parameter spaces illustrate numerically the transition from unpolarized to polarized
stationary states depending on the parameters.

**5 Numerical trajectories**

In this Section, we add stochastic fluctuations to the system given by (3)–(4)–(5)–(6)–(7) to generate numerical trajectories. This Section is organized as follows. In Subsection 5.1, we introduce the stochastic model based on the additive perturbation of the advection field by a colored gaussian noise. In Subsection 5.2, single-cell simulations illustrate different regimes. Finally, Subsection 5.3 focuses on their characterization at the population level.

**5.1 Stochastic model**

This Subsection aims at introducing the stochastic version of the model. Its underlying
hypothesis is that the main source of stochasticity apparent at the cell scale lies in
the advection of polarity cues by the actin cytoskeleton. Indeed, this phenomenon is
mediated by numerous other molecules involved in physical and chemical reactions. We
assume that all these processes may be taken into account by an additive perturbation
<figure><img src="image_5.png" /><figcaption>Fig. 4 Parameter spaces of polarity indices at stationary state as a function of the boundary pressure amplitude $\chi$ and either the half-saturation threshold $c_t$ (top) or the molecular diffusion coefficient $D$ (bottom). The initial condition is given by Equation (42). (a-c): Resultant Vector Length. (b-d): Relative Standard Deviation. The relation $\chi = 4D$ is shown in red in the $(D, \chi)$ parameter spaces. Parameters: $(\gamma, D, c_t) = (0.1, 0.3, 0.4)$; threshold for numerical stationary state: $10^{-6}$; $(\Delta t, N_\theta) = (0.01, 120)$, $\Delta r = \Delta \theta = \frac{2\pi}{N_\theta} \approx 0.05$</figcaption></figure>

of the advection velocity using a white noise in time that carries Gaussian space correlations. More precisely, the advection-diffusion equation (3) writes

$$ \left\{ \begin{array}{ll} \partial_t c + \operatorname{div} (c \mathbf{u}_{\text{sto}} - D \nabla c) = 0 & \text{in } \Omega, \\ (D \nabla c - c \mathbf{u}_{\text{sto}}) \cdot \mathbf{n} = 0 & \text{on } \partial \Omega. \end{array} \right. \qquad (44) $$

The advection velocity field $\mathbf{u}_{\text{sto}}$ is given by

$$ \mathbf{u}_{\text{sto}} = \mathbf{u} + \alpha \mathrm{d}\mathbf{w}_t^{Q_{\sigma}}, \qquad (45) $$

where $(\mathbf{w}_t^{Q_{\sigma}})_{t \ge 0} = (w_{1t}^{Q_{\sigma}}, w_{2t}^{Q_{\sigma}})_{t \ge 0}$ is the realization of a $Q_{\sigma}$-Wiener process in $L^2(\Omega)^2$, and $\alpha \ge 0$ is the noise amplitude. Note that the molecular mass remains conserved thanks to the boundary conditions. The operator $Q_{\sigma}$ is assumed to admit a Gaussian kernel $q_{\sigma}$ with correlation length $\sigma > 0$ given for all $(\mathbf{x}, \mathbf{y}) \in \Omega^2$ by

$$ q_{\sigma}(\mathbf{x}, \mathbf{y}) = \frac{1}{4\sigma^2} e^{-\frac{\pi}{4\sigma^2} \|\mathbf{x}-\mathbf{y}\|^2}. \qquad (46) $$
The covariance between two points $(t, \mathbf{x})$ and $(s, \mathbf{y})$ in $\mathbb{R}_+ \times \Omega$ then writes

$$ \mathbb{E}[W_t^{Q_{\sigma}}(\mathbf{x})W_s^{Q_{\sigma}}(\mathbf{y})] = t \wedge s q_{\sigma}(\mathbf{x}, \mathbf{y}). $$

We refer to [28, 29] for the general definition and properties of the Q-Wiener process. The numerical simulations of the noise are realized as follows: at each time step, the advection velocity is computed respectively at the radial and angular mesh interfaces for the radial and angular fluxes. The noise term is computed in cartesian coordinates over all these nodes to account for spatial correlations. The required radial and angular components of the noise are then derived only on the necessary points with respect to the flux, see Appendix C.3 for more details.

In the following simulations, the parameter $\chi$ vary, while the others are fixed with $(\gamma, c_t, D) = (0.1, 0.4, 0.3)$ for the dynamics, and $(\alpha, \sigma) = (0.5, 1)$ for the noise. The final time is $T_f = 100$. The initial condition is taken as the radially-symmetric stationary state of the deterministic model given by (23)–(25). This ensures that all cells are initially unpolarized, so that their dynamics result from stochastic fluctuations either reinforced or attenuated by the deterministic dynamics.

## 5.2 Single-cell dynamics

This Subsection aims at illustrating single-cell dynamics at different regimes, depending on the chosen set of parameters. More precisely, we consider $\chi \in \{0.001, 0.94, 1.1\}$. Figure 5 shows the resulting simulations in the three cases. Subfigure 5-(a) depicts the molecular concentration with stochastic perturbations for increasing value of $\chi$. If no clear distinction can be made in fixed time for the two first cases, the third case shows a clearly polarised state. Subfigure 5-(b) illustrates the time evolution of the cell velocity modules. The case where $\chi = 0.001$ is quantified on the left axis, and shows velocity module values close to zero all the time. The two other cases are quantified on the right axis. When $\chi = 0.94$, the velocity module shows excursions far from zero. When $\chi = 1.1$, there is a net velocity module increase up to a stabilized regime. Subfigure 5-(c) shows the corresponding histograms. The resulting trajectories are given in Subfigure 5-(d). In each case, the trajectory is colored in red (resp. blue) for velocity modules below 0.004 (resp. above). Moreover, trajectories are split in persistent subpaths characterized by a maximal angular drift of $\pi/2$, and delimited by the black horizontal segments. One can see that as $\chi$ increases the trajectories get straighter and roam a widening space. The case $\chi = 0.94$ also illustrates how loss of directional persistence is concomitant with lower velocity modules. This reveals how the dynamics reinforces more durably local polarity spikes with increasing coupling between the fluid and molecular levels.

## 5.3 Population-based characterization

The observed regimes can be characterized at the population level using several quantities.
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Illustration of three types of trajectories, with $(T_f, \Delta t) = (200, 0.01)$, $(\gamma, D, c_t) = (0.1, 0.3, 0.4)$, $(\alpha, \sigma) = (0.5, 1)$ and $\chi \in \{0.001, 0.94, 1.1\}$ (respectively left, middle and right). The initial condition is given by (23)–(25). (a): molecular concentration at final time in each condition. (b): Time evolution of the cell velocity modules, with left axis for $\chi = 0.001$ and right axis for $\chi \in \{0.94, 1.1\}$. (c): Velocity modules histograms. (d): Cell trajectories colored in red for velocity modules below 0.004 (respectively in blue above). The black horizontal lines delimitate persistent subpaths with maximal absolute angular drifts of $\pi/2$</figcaption></figure>

* The Ensemble-Averaged Instantaneous Speed (EAIS) allows both to quantify the mean velocity in the population, and to characterize the ergodicity of the system, that is to say whether the dynamics intrinsically changes during the time period.

* The typical orientation correlation time is the timescale below which the directional orientation is correlated. It is obtained from the tangent-tangent correlation function that maps any time shift $\tau \ge 0$ to the population-averaged and time-averaged cosines of directional angles separated by the time shift, denoted by $\cos((v(t), v(t+\tau)))$. Its tangent at origin meets the time axis at the typical orientation correlation time.
* The *Time-Averaged Squared Displacements* (TASD) is the average of squared displacements over all possible time shifts for a given cell. For cell $\ell$ with position $\mathbf{z}_\ell$ and $1 \le n \le N_f$, it writes:

$$TASD_{\ell}(n\Delta t) := \frac{1}{N_f - n + 1} \sum_{k=0}^{N_f - n} |\mathbf{z}_{\ell}((k+n)\Delta t) - \mathbf{z}_{\ell}(k\Delta t)|^2.$$

* The *Mean Squared Displacements* (MSD) is the average of TASD over the cell population: for $1 \le n \le N_f$, one has

$$MSD(n\Delta t) := \frac{1}{N_{cell}} \sum_{\ell=1}^{N_{cell}} TASD_{\ell}(n\Delta t).$$

Fitting a power law $t \mapsto \delta t^\beta$ to these curves provides the index $\beta$ characterizing motion: it equals 1 for a purely diffusive motion, and 2 for a purely ballistic motion.

* The *velocity histograms* inform on both the range and orientation of cell velocities within the population.

* Finally, the typical directionality indices used in [1] are computed. They are obtained by splitting each trajectory in directional paths: a new subpath is defined whenever the orientation of motion drifts from the initial direction by at least $\pi/2$. Then, each path is associated with its duration $T_{pe}$ and mean velocity $V_{pe}$, respectively called *persistence time* and *persistence velocity*. The UCSP law consists in an exponential relationship between both quantities. Therefore, such a fit is performed.

In the following, these quantities are computed for populations consisting of 100 cells in each cohort, sharing the radially symmetric initial condition given by (23)–(25). The final time is $T_f = 100$.

The case $\chi = 0.001$ is illustrated in Figure 6. In Subfigure 6-(a), the EAIS quickly reaches a plateau indicating ergodicity. The 2D velocity histogram in Subfigure 6-(b) presents a Gaussian distribution centered in the null velocity. In Subfigure 6-(c), the directional correlation of the displacement rapidly drops to zero, while more individual variations appear for larger time shifts, reflecting the smaller number of data points. The associated typical correlation time is $t_c = 0.05$. In the TASD Subfigure 6-(d), some individual curves differentiate from the typical population behaviour and explore a larger distance over time. This is due to the noise driving the dynamics away from the deterministic behaviour. Fitting a power law to the MSD curve yielded $\beta = 1.17$, characterizing a nearly purely diffusive motion on average. This case is therefore refered to as the *diffusive case*. Subfigure 6-(e) shows three trajectories with their individual TASD index, to illustrate population variability. Finally, Subfigure 6-(f) shows the relation between persistence times and velocities. Note that many data points are very close to zero for small persistence velocity values. This suggests to fit a shifted exponential function that is null at $v = 0$, as shown in the red line. One can see that the trend of the data is not well-captured by the exponential function.
The case $\chi = 1.1$ is illustrated in Figure 7. In Subfigure 7-(a), the EAIS shows that ergodicity is reached at $t \approx 60$ with a positive mean velocity module. The 2D velocity histogram given in Subfigure 7-(b) indicates that cells mainly have positive velocities, while lower values correspond to the initial accelerating phase. In Subfigure 7-(c), the mean orientation correlation remains close to 1 for almost all time shifts. The typical correlation time is equal to 2.99, although the individual curves show longer-lasting correlation. The TASD and MSD curves shown in Subfigure 7-(d) are steep with an MSD exponent $\beta = 1.97$, characterizing a ballistic motion. This case is therefore called the *persistent case*. Subfigure 7-(e) illustrates three trajectories together with their individual TASD index. Finally, Subfigure 7-(f) shows a robust distribution among cells, with small values of persistence time and velocity corresponding to the initial phase, and larger values associated with persistent motion. The smaller number of points reflects the presence of long persistent subpaths. Given the observed trend, fitting an exponential function has no meaning in this case.

Lastly, the case where $\chi = 0.94$ ought to represent an in-between dynamics, and is illustrated in Figure 8. In Subfigure 8-(a), one can see small variations in the EAIS, although the system shows a globally ergodic behaviour. The 2D velocity histogram in Subfigure 8-(b) is qualitatively similar as the one in the diffusive case, although the velocity range is wider. The averaged tangent-tangent correlation curve in Subfigure 8-(c) drops to zero slower than the diffusive case, with a typical correlation time $t_c = 0.3$. However, note that the individual curves show a large variability, indicating the frequent occurrence of more persistent displacements. The TASD and MSD curves in Subfigure 8-(d) are associated with an index $\beta = 1.48$, and this case is now refered to as the *intermittent case*. Subfigure 8-(e) illustrates three trajectories together with their individual TASD index. Finally, Subfigure 8-(f) shows the link between persistence velocities and times, with a satisfactory fitting of an exponential curve between both values.

# 6 Discussion

## 6.1 Deterministic model

### Modelling

The construction of the deterministic model relies on several assumptions. Firstly, the choice of a solid disk for the cell domain is a limitation. It was made to focus on the interplay between the actin cytoskeleton and a molecular specie in the intracellular dynamics. Secondly, the main novelty in our approach lies in describing the actin cytoskeleton using an active fluid. The constitutive equations rely on the following ingredients: actin polymerization at the boundary, depolymerization in the domain, and friction forces arising from the actin flowing from the boundary towards the center. Reduction of the problem towards the Darcy-Poisson problem (4)-(5) resulted from neglecting fluid viscosity in comparison with friction forces, which is relevant for a crawling cell, and assuming that the pressure instantaneously forced a constant actin density, following existing works [30]. The limit in constant actin density is a simplifying assumption in line with considering a solid domain: the actin activity translates
<figure><img src="image_7.png" /><figcaption>Fig. 6 Population-based characterization of cell trajectories for $\chi = 0.001$, with initial condition given by (23)–(25). (a) Ensemble-Averaged Instantaneous Velocity. (b) Bidimensional velocity histogram. (c) Tangent-Tangent correlation function with typical orientation correlation time $t_c = 0.05$. (d) Individual Time-Averaged Squared Displacements (TASD) and Mean Squared Displacement (MSD) with index $\beta = 1.17$, indicating a diffusive regime. (e) Three representative trajectories together with their individual TASD index. (f) Persistence time as a function of persistence velocity obtained when sequentially splitting trajectories in maximal subpaths having angular deviation from the initial direction below $\frac{\pi}{2}$. Exponential fitting yields $v \mapsto 0.87e^{3.87 \times 10^5 v} - 1$. Parameters: $(T_f, \Delta t) = (200, 0.01)$, $(\gamma, D, c_t) = (0.1, 0.3, 0.4)$, and $(\delta, \sigma) = (0.5, 1)$</figcaption></figure>

directly into a pressure, as no deformation is possible. The full derivation is detailed in Appendix A. The model can be further compared with [1]. Indeed, in the cited work the mean space-averaged actin velocity $u^*$ in 1D is given by Equation (2). In our model, computing the space-averaged actin velocity given by Equations (4)-(6) for $k_p$ given by (7) yields

$$ \int_{\Omega} u(t, \mathbf{x}) \, d\mathbf{x} = \chi(\gamma\pi - 1) \int_{\partial\Omega} \frac{c(t, \mathbf{x})}{1+c(t, \mathbf{x})} \mathbf{n} \, d\sigma(\mathbf{x}), $$

recovering a similar expression, where, up to spatial normalization, the coupling parameter $\beta$ from Equation (2) corresponds to $\chi(\gamma\pi + 1)$. Moreover, it is possible to go beyond the simplifying assumption $v = -u^*$ from [1] by writing

$$ \int_{\Omega} u(t, \mathbf{x}) \, d\mathbf{x} = -\pi v(t) + \chi \int_{\partial\Omega} \frac{c(t, \mathbf{x})}{1+c(t, \mathbf{x})} \mathbf{n} \, d\sigma(\mathbf{x}). $$
<figure><img src="image_8.png" /><figcaption>**Fig. 7** Population-based characterization of cell trajectories for $\chi = 1.1$, with initial condition given by (23)–(25). (a) Ensemble-Averaged Instantaneous Velocity. (b) Bidimensional velocity histogram. (c) Tangent-Tangent correlation function with typical orientation correlation time $t_c = 2.99$. (d) Individual Time-Averaged Squared Displacements (TASD) and Mean Squared Displacement (MSD) with index $\beta = 1.97$, indicating a diffusive regime. (e) Three representative trajectories together with their individual TASD index. (f) Persistence time as a function of persistence velocity obtained when sequentially splitting trajectories in maximal subpaths having angular deviation from the initial direction below $\pi/2$. Parameters: $(T_f, \Delta t) = (200, 0.01)$, $(\gamma, D, c_t) = (0.1, 0.3, 0.4)$, and $(\delta, \sigma) = (0.5, 1)$</figcaption></figure>

This highlights how the model now captures the dissociation between the internal activity and the external motion. Indeed, on a nonadhesive substrate, a polarised cell could generate sustained actin flows while not moving, since no friction forces convert the internal activity into driving forces. This amounts to taking $\gamma = 0$: the cell velocity is null, while a net global actin velocity may still exist. Furthermore, when $\gamma \neq 0$, then

$$|v(t)| = \frac{\gamma}{1 + \pi\gamma} \left| \int_{\Omega} u(t, \mathbf{x}) \, d\mathbf{x} \right| < \left| \int_{\Omega} u(t, \mathbf{x}) \, d\mathbf{x} \right|,$$

in coherence with the expected cell behaviour moving slower than its driving internal flows. To summarize, our 2D modelling approach allows to justify the 1D model from [1]. Furthermore, our model allows to decipher local effects related to actin flows from global effects linked with cell displacement. The model relies on 4 parameters: the molecular diffusion coefficient $D$, the molecular activation saturation coefficient $c_t$, the minimal boundary pressure value $\chi$, the velocity coefficient $\gamma$.
<figure><img src="image_9.png" /><figcaption>**Fig. 8** Population-based characterization of cell trajectories for $\chi = 0.94$, with initial condition given by (23)–(25). (a) Ensemble-Averaged Instantaneous Velocity. (b) Bidimensional velocity histogram. (c) Tangent-Tangent correlation function with typical orientation correlation time $t_c = 0.3$. (d) Individual Time-Averaged Squared Displacements (TASD) and Mean Squared Displacement (MSD) with index $\beta = 1.48$, indicating a diffusive regime. (e) Three representative trajectories together with their individual TASD index. (f) Persistence time as a function of persistence velocity obtained when sequentially splitting trajectories in maximal subpaths having angular deviation from the initial direction below $\pi/2$. Exponential fitting yields $v \mapsto 0.93e^{600v} - 1$. Parameters: $(T_f, \Delta t) = (200, 0.01)$, $(\gamma, D, c_t) = (0.1, 0.3, 0.4)$, and $(\delta, \sigma) = (0.5, 1)$</figcaption></figure>

## Mathematical analysis

The present work can be compared with the results obtained for the close 1D model studied in [27] for a linear boundary pressure term and an actin velocity depending only on the pressure gradient. In this case, a resting stationary state was put to light for a subcritical mass. In the supercritical case and for some initial condition, it was proved that blowup occured in finite time. From the modelling viewpoint, it was considered as a degenerate polarised state. In the present work, under some parameter conditions, we proved existence of both resting and traveling stationary states and linear stability of the resting stationary state. We did not reproduce numerically a blowup. This may be due to the Hill activation term ensuring boundedness of the pressure boundary terms.
## Numerical simulations

Numerical simulations of the deterministic model yielded both polarized and unpolarized cell regimes, depending on the balance between the strength of molecular inhibition of the inward flow and both molecular diffusion and molecular activation.

The numerical simulations also allowed to investigate the mathematical conditions on parameters obtained in Propositions 8, 10 and 12. They suggested that existence of both stationary states go beyond the theoretical conditions. This suggests that they could reflect technical limitations in the proofs. Indeed, the model presents a complex competition between different trends: the molecular diffusion classically competes with their advection, but the advection field itself results from a trade-off between the self-reinforcing molecular asymmetry at the boundary, and depolymerization drifting the molecules towards the cell center. In this context, identifying strong conditions for the existence of a given regime is a challenging task. Numerical simulations also indicate that the threshold on $\chi$ for linear stability of the resting stationary state is not optimal.

From the numerical viewpoint, using a linear upwind scheme of order 1 prevented having spurious oscillations of the solution, but also induced some numerical diffusion. We found that taking $N_\theta = 200$, hence $\Delta r = \Delta \theta \simeq 0.03$, was sufficient to have a tractable computational cost while reducing the error, as seen in Figure 2.

## 6.2 Stochastic modelling

### Modelling

The key point in the stochastic modelling was choosing the source of stochasticity. The intracellular motility machinery is known to be very complex, multiscale and stochastic, so that several choices could be relevant. In [1], two stochastic sources were considered. Firstly, the global actin velocity was assumed to follow a stochastic diffusion with a velocity module fluctuating around its mean value independently of its orientation, which diffuses with an intensity proportional to the inverse of the velocity module. Secondly, an additive spatially-colored Gaussian noise perturbed the molecular flux. A series of formal steps then led to combining both sources of stochasticity in an Active Brownian Particle model.

If the rigorous justification of the derivation of the SDE is of interest, it goes beyond the scope of the present work based on a fully 2D framework. Therefore, we considered an additive perturbation of the advection velocity in the molecular problem using a space-colored noise derived from a $Q$-Wiener process. The noise amplitude was taken small enough with respect to the deterministic dynamics. Spatial correlations were taken large enough to generate cell-scale variations in the actin-based molecular transport, reflecting variability in the cytoskeletal dynamics. This yields several differences with the preceding work. First, only one source of noise is considered and applies to the advection velocity. However, in this case it jointly acts on the velocity magnitude and orientation. Also, stochasticity is no longer directly involved in the cell velocity apart from terms arising from the molecular scale.

Other sources of spatial stochasticity can be considered. The other natural option was to consider a similar spatial noise directly in the molecular flux. Numerical simulations in this setting yielded qualitatively similar behaviours. One can also assume that
major fluctuations are directly occuring in the actin polymerization at the cell membrane. This assumption could be related to other Active Brownian Particle models relying on polymerization that could justify a boundary noise term structured by the orientation [11, 31]. Finally, a more straightforward assumption relies on the stochastic perturbation of the cell velocity term directly using a white noise. This would result in a global perturbation of the bias induced by cell motion. However, derivation of a SDE would be easier and should constitute a first step in this direction. These ideas will be explored in future works.

## Single-cell dynamics

In the stochastic framework, the intracellular molecular concentration exhibited cell-scale perturbations associated with transient reduction of polarity or change in the direction of motion, in competition with the fluid-molecular coupling sustaining polarity (Figure 5). Typical trajectories illustrate how their morphology change with increasing coupling, with cell polarization associated with more persistent trajectories. They can be qualitatively compared with the experimental trajectories presented in [1]. Figure 9 depicts the diffusive, persistent and intermittent regimes in trajectories and velocity module histograms. Coloration in red (resp. blue) corresponds to a velocity module below 4 µm/min (resp. above). We observe similarities in our simulations in terms of both straightness of trajectories and explored space. Moreover, splitting trajectories with respect to either velocity modules or angular drift showed that changes in direction are associated with lower velocities. In the individual velocity module histograms in Subfigure 5(c), we did not find the bimodal histogram corresponding to an intermittent regime, as observed experimentally.

<figure><img src="image_10.png" /><figcaption>**Fig. 9** Experimental trajectories from [1], together with their velocity module histograms. The red circles account for confidence interval of 3 µm. Coloration in red (resp. blue) accounts for $||V|| > 4$ µm/min (resp. $||V|| < 4$ µm/min). The histograms are shown together with the theoretical distribution as found in [1]</figcaption></figure>
Population-based characterization

At the population scale, it was possible to derive more robust quantities in the three
considered cases (Figures 6–8–7). Exponents obtained from Mean Squared Displace-
ments characterized behaviours from diffusive to ballistic. The Ensemble-Average
Instantaneous Velocities showed the overall ergodic dynamics – a bit less in the inter-
mediate regime due to the competition between both trends. Depending on both the
deterministic coupling and the stochastic pertubation, the directional persistence of
motion can be either rapidly lost or strongly sustained. This translates in velocity his-
tograms ranging from a Gaussian-like distribution centered in the null velocity to an
annulus indicating consistent net displacement. The typical correlation times derived
from the tangent-tangent correlation curves were ordered as expected. However, it was
undervalued in the persistent case, where the curve itself showed sustained directional
correlation for most part of the simulation. The interplay between velocity magnitude
and directional persistence were further quantified jointly by splitting trajectories with
respect to their angular drifts, as in [1]. In the diffusive case, there was no point in
testing an exponential correlation as many values were very close to zero. In the persis-
tent case, the net separation between initial diffusive subpaths and persistent subpaths
made the fitting of an exponential function irrelevant, although no incompatibility
is observed either. In the intermittent case, we were able to recover the exponential
coupling between typical persistence time and their associated persistence velocity.

7 Conclusions

This work aims at formulating a 2D model for single-cell migration based on the mechano-chemical coupling between the actin cytoskeleton and a regulating molecular specie as presented in [1]. The fulfilled objective was to provide a fully spatial framework and to justify the phenomenological formulation of the actin velocity considered therein. The actin velocity field is described as an active fluid with constitutive equations based on both actin polymerization and friction forces arising from actin flows and cell displacements. Consequently, cell polarization occurs thanks to the coupling between the fluid and the molecular specie, and it is associated with persistence of motion. The mathematical analysis of the resulting nonlinear nonlocal model characterizes both resting and traveling stationary regimes, as illustrated by numerical simulations. Numerical trajectories obtained from stochastic perturbations of the model allows to recover a variety of behaviours ranging from diffusive to ballistic, characterized in terms of velocity magnitude and directionality. In an intermediate setting, the model is also able to recover a generic exponential coupling between cell speed and persistence as observed experimentally.

Further works ought to take advantage of both mechanical and chemical descrip-
tions of the intracellular dynamics to investigate the cell interaction with its environ-
ment. Indeed, it allows to take into account each type of interaction at the right level.
The whole intracellular activity then adapts to the external cue. As a consequence,
investigating phenomena such as chemotaxis or guidance by mechanical obstacles will
be of interest. More precise comparison to trajectories will be investigated in a future
work. Finally, perspectives include investigating the link between such kind of models
with given stochastic perturbation with active particle model as in [1].

**Acknowledgements.** The authors are very grateful to C. Poignard and C. Alamichel for very helpful discussions and remarks, and to L. Facq for his help in script parallelization.

**Declarations**

**Funding**

The authors did not receive support from any organization for the submitted work

**Competing interests**

The authors have no competing interests to declare that are relevant to the content
of this article

**Ethics approval and consent to participate**

Not applicable

**Consent for publication**

All authors gave their consent for publication

**Data availability**

Not applicable

**Materials availability**

Not applicable

**Code availability**

The Python script used to perform the numerical simulations is available upon request to the corresponding author.

**Author contribution**

Conceptualization: Nicolas Meunier and Raphaël Voituriez ; Methodology: Nicolas Meunier ; Formal analysis and investigation: Christèle Etchegaray and Nicolas Meunier ; Writing - original draft preparation: Christèle Etchegaray and Nicolas Meunier ; Writing - review and editing: Christèle Etchegaray and Nicolas Meunier ; Supervision: Nicolas Meunier and Raphaël Voituriez
# Appendix A Modelling approach

This Section aims at introducing the model based on physical laws from which we formally derive the reduced model (3)-(4)-(5)-(6). First, the modelling principles underlying the approach are discussed.

## The cytoskeleton

The actin cytoskeleton is a highly complex and dynamical medium composed of actin filaments organized in bundles and networks by their transient molecular interaction. The chosen approach bypasses the individual filaments dynamics and adopts a continuum description based on the active gel theory. Initially applied to the study of liquid crystals [32], it consists in assuming large length and long time scales compared to the cross-linking dynamics, and in preserving conservation laws and symmetry properties [33, 34]. Several models of migration have used this approach [6, 7, 35], which have helped explaining the actin retrograde flow, building a qualitative description of an advancing lamellipodium, and estimating active stresses on the keratocyte cell [34]. We refer to [34, 36–38] for a general description of the theory. Under the assumption of long time scales compared with the individual filaments dynamics, the cytoskeleton is described by a fluid [7]. Considering a crawling cell, has two additional consequences. Firstly, any thickness effect is neglected, yielding a bidimensional description. Secondly, the cytosol is not described, as its interaction with the cytoskeleton is negligible compared with the cytoskeletal forces applied on the substrate [38, 39]. This leads to considering a monophasic fluid description of the cytoskeleton.

## Interaction with a molecular regulator of the actin dynamics

It is known that a large number of molecular species interact and form signalling pathways that regulate cellular functions related to motility. As a consequence, spatial heterogeneity in molecular concentration profiles results in a polarised activity that may induce a displacement. In the present approach, a single generic molecule is considered as a regulator of the actin cytoskeleton activity. The main assumption is that it locally inhibits actin polymerization at the membrane: its aggregation in an activated form in a part of the membrane characterizes the cell rear. This specie could be an antagonist to polymerization-inducing molecules (Rac1, Cdc42), such as RhoA, Arpin, or even myosin II (see [40, 41]). Moreover, it has been observed that some molecules can bind to actin filaments, and be transported by the inward actin flows. This mechanism was proposed in [1] as a way of initiating and maintaining a spatial molecular gradient (see e.g the Graphical Abstract in the cited work. This is the chosen interaction here.

*Remark 2* Let us emphasize that the aim is not to identify a single molecular specie responsible for universally downregulating actin polymerization. The key idea is that the nature of this regulator can vary a lot among cell types, while the general mechanism remains the same.
## A.1 Fluid equations

Denote $t \mapsto \Omega(t)$ the moving solid disk that defines the cell, and $V : \mathbb{R}_+ \to \mathbb{R}^2$ its global velocity. For $t > 0$ and $X \in \Omega(t)$, let us define the fluid pressure $\tilde{P}(t, X)$, and the actin velocity field $\tilde{U}(t, X)$. The viscous fluid that models the actin cytoskeleton is characterized by two physical conservation laws: the force balance principle and the continuity equation on the actin density.

### A.1.1 Force balance principle

In a microscopic setting with low Reynolds number, and in the frame of reference of the lab, the force balance principle writes

$$ -\mathrm{div} \, \boldsymbol{\sigma} = f(t, \mathbf{X}) \quad \text{in } \Omega(t), $$

where the Stokes stress tensor $\boldsymbol{\sigma}$ is given by

$$ \begin{aligned} \boldsymbol{\sigma} &= \eta \left( \nabla(\tilde{\mathbf{U}}(t, \mathbf{X}) + \mathbf{V}(t)) + {}^t\nabla(\tilde{\mathbf{U}}(t, \mathbf{X}) + \mathbf{V}(t)) \right) - \tilde{P}(t, \mathbf{X}) \mathrm{Id}, \\ &= \eta \left( \nabla\tilde{\mathbf{U}}(t, \mathbf{X}) + {}^t\nabla\tilde{\mathbf{U}}(t, \mathbf{X}) \right) - \tilde{P}(t, \mathbf{X}) \mathrm{Id}, \end{aligned} $$

with $\eta \in \mathbb{R}_+$ is the fluid viscosity. The force $f$ arising from the actin filaments rubbing on the substrate when the actin flow is given by

$$ f(t, \mathbf{X}) = -\xi(\tilde{\mathbf{U}}(t, \mathbf{X}) + \mathbf{V}(t)), $$

where $\xi \in \mathbb{R}_+$ is an effective friction coefficient. This leads to the Stokes equation

$$ -\eta \Delta \tilde{\mathbf{U}}(t, \mathbf{X}) + \nabla \tilde{P}(t, \mathbf{X}) = -\xi(\tilde{\mathbf{U}}(t, \mathbf{X}) + \mathbf{V}(t)) \quad \text{in } \Omega(t). $$

The boundary condition is defined using the normal component of the stress tensor:

$$ \boldsymbol{\sigma} \cdot \mathbf{n} = \eta \left( \nabla \tilde{\mathbf{U}}(t, \mathbf{X}) + {}^t\nabla \tilde{\mathbf{U}}(t, \mathbf{X}) \right) \cdot \mathbf{n} - \tilde{P}(t, \mathbf{X}) \mathbf{n} = -K_p \mathbf{n}, $$

where $K_p$ corresponds to the force related to actin polymerization at the boundary.

### A.1.2 Continuity equation

For $\varepsilon > 0$ a small parameter homogeneous to the inverse of a pressure, we introduce the actin density $\tilde{\rho}_\varepsilon(t, X)$, defined for $t > 0$ and $X \in \Omega(t)$. The continuity equation inside the domain writes

$$ \frac{\partial}{\partial t} \tilde{\rho}_{\varepsilon}(t, \mathbf{X}) + \mathrm{div} (\tilde{\rho}_{\varepsilon}(t, \mathbf{X})(\tilde{\mathbf{U}}(t, \mathbf{X}) + \mathbf{V}(t))) = -K_d \tilde{\rho}_{\varepsilon}(t, \mathbf{X}) \quad \text{in } \Omega(t), $$

where $K_d$ is a constant depolymerization rate acting on the whole domain and the actin advection field $\tilde{\mathbf{U}}(t, X) + \mathbf{V}(t)$ is expressed in the frame of reference of the lab.
To account for polymerization taking place at the edge of the cell, we impose the following Dirichlet boundary condition:

$$ \tilde{\rho}_{\varepsilon}(t, \mathbf{X}) = \rho_0 + \varepsilon K_p \quad \text{on } \partial\Omega(t), $$

where $\rho_0 \in \mathbb{R}_+$ is a constant.

*Remark 3* A more realistic approach would consist in assuming that the polymerization occurs in an annulus initially defined as $\bar{\Omega}(0) \setminus B(0, R - \lambda)$, of width $\lambda > 0$, see e.g [39].

Following [30], we make the assumption that the pressure is linear in the difference between the actin density and the resting value $\rho_0$:

$$ \tilde{P}_{\varepsilon} = \frac{1}{\varepsilon}(\tilde{\rho}_{\varepsilon} - \rho_0) \quad \text{in } \Omega(t). $$

This means that the pressure tends to maintain a fluid density close to $\rho_0$.

### A.1.3 Cell velocity

The global cell velocity derives from the force balance principle applied to the cell. Two types of friction force arise, generated respectively by the actin flows inside the cell and by the cell displacement itself. The latter is related to both the adhesion on the substrate, and the drag of inert organites as a passive load. It writes

$$ \xi_1 \mathbf{V}(t) + \xi \int_{\Omega(t)} \tilde{\mathbf{U}}(t, \mathbf{x}) + \mathbf{V}(t) \, d\mathbf{x} = 0, $$

where $\xi_1$ is a global cell friction coefficient. This is equivalent to

$$ \mathbf{V}(t) = - \frac{\xi}{\xi_1 + \xi |\Omega(0)|} \int_{\Omega(t)} \tilde{\mathbf{U}}(t, \mathbf{X}) \, d\mathbf{X}. $$

### A.1.4 Formulation in a motionless domain

The fluid problem is now written on the fixed domain $\Omega := \Omega(0) = B(0, R)$, with $R > 0$. As $(\tilde{P}, \tilde{\mathbf{U}}, \tilde{\rho}_{\varepsilon})$ are defined on $\Omega(t)$, we denote $(P, \mathbf{U}, \rho_{\varepsilon})$ their analogous on $\Omega$. Since the domain velocity is a global quantity, only the advection velocity in the continuity equation is modified by the change of reference domain, and the expression for the cell velocity remains the same in the new reference domain. It writes

$$ \begin{cases} -\eta \Delta \mathbf{U}(t, \mathbf{x}) + \nabla P(t, \mathbf{x}) = -\xi(\mathbf{U}(t, \mathbf{x}) + \mathbf{V}(t)) & \text{in } \Omega, \\ \eta (\nabla \mathbf{U}(t, \mathbf{x}) + {}^t\nabla \mathbf{U}(t, \mathbf{x})) \cdot \mathbf{n} - P(t, \mathbf{x})\mathbf{n} = -K_p \mathbf{n} & \text{on } \partial\Omega, \end{cases} $$
together with

$$
\left\{
\begin{array}{l}
\displaystyle \frac{\partial}{\partial t} \rho_{\varepsilon}(t, \mathbf{x}) + \operatorname{div} (\rho_{\varepsilon}(t, \mathbf{x}) \mathbf{U}(t, \mathbf{x})) = -K_d \rho_{\varepsilon}(t, \mathbf{x}) \quad \text{in } \Omega, \\
\displaystyle \rho_{\varepsilon}(t, \mathbf{x}) = \rho_0 + \varepsilon K_p \quad \text{on } \partial\Omega, \\
\displaystyle P_{\varepsilon} = \frac{1}{\varepsilon}(\rho_{\varepsilon} - \rho_0) \quad \text{in } \Omega.
\end{array}
\right.
$$

A.2 Molecular dynamics

For $t \in \mathbb{R}_+$ and $\mathbf{x} \in \Omega$, we denote by $C(t, \mathbf{x})$ the concentration of a molecular specie inside $\Omega$. It evolves following an advection-diffusion dynamics with no-flux boundary condition, ensuring conservation of the total molecular mass $M$. The active transport by the cytoskeleton is described by the advection field equal to the fluid velocity field. Note that a justification of this equation for molecular transport can be found in [42]. For $t > 0$, the model writes

$$
\left\{
\begin{array}{ll}
\partial_t C + \operatorname{div} (C\mathbf{U} - D'\nabla C) = 0 & \text{for } \mathbf{x} \in \Omega, \\
(C\mathbf{U} - D'\nabla C) \cdot \mathbf{n} = 0 & \text{for } \mathbf{x} \in \partial\Omega,
\end{array}
\right.
$$

where $D'$ is the diffusion coefficient, and $\mathbf{n}$ the outward unit normal to the boundary. Finally, we assume that activated molecules at the cell boundary are able to inhibit actin polymerization. More precisely, $K_p$ is taken as a Hill function given by $K_p(C(t, \mathbf{x})) = \bar{k}C(t, \mathbf{x})/(c_t + C(t, \mathbf{x}))$, where $\bar{k}, c_t > 0$.

A.3 Formal derivation of the reduced model

A.3.1 Non-dimensionalization

The non-dimensionalization of the problem is based on the following typical quantities:

$$
\begin{align*}
\bar{L} &= R, & \bar{T} &= K_d^{-1}, & \bar{V} &= \frac{\bar{L}}{T} = RK_d, \\
\bar{F} &= \xi \bar{V} &= \xi R K_d, & \bar{P} &= \bar{F} \bar{L} = \xi R^2 K_d, & \bar{C} &= M.
\end{align*}
$$

The dimensionless equations are now written using lower case letters for the dimen-
sionless variables and keeping the notation $\epsilon$ for simplicity. Let $\Omega := B(0, 1)$ be the
cell domain. Then, the Stokes Problem rewrites

$$
\left\{
\begin{array}{ll}
-\tilde{\eta}\Delta\mathbf{u} + \nabla p = -(\mathbf{u} + \mathbf{v}) & \text{in } \Omega, \\
\tilde{\eta}(\nabla\mathbf{u} + {}^t\nabla\mathbf{u}) \cdot \mathbf{n} - p\mathbf{n} = -\chi k_p(c)\mathbf{n} & \text{on } \partial\Omega,
\end{array}
\right.
$$

with $\tilde{\eta} := \frac{\eta}{\xi R^2}$, $\chi := \frac{\bar{k}}{\bar{P}}$, and $k_p(c) = \frac{c}{c_t+c}$. The continuity equation rewrites

$$
\left\{
\begin{array}{l}
\displaystyle \frac{\partial}{\partial t} \rho_{\varepsilon} + \operatorname{div} (\rho_{\varepsilon} \mathbf{u}) = -\rho_{\varepsilon} \quad \text{in } \Omega, \\
\rho_{\varepsilon} = \rho_0 + \varepsilon \bar{k} k_p(c) \quad \text{on } \partial \Omega, \\
\bar{P} p_{\varepsilon} = \frac{1}{\varepsilon} (\rho_{\varepsilon} - \rho_0) \quad \text{in } \Omega.
\end{array}
\right.
$$
The molecular problem now writes

$$ \begin{cases} \partial_t c + \operatorname{div} (c\mathbf{u} - D\nabla c) = 0 & \text{for } \mathbf{x} \in \Omega, \\ (c\mathbf{u} - D\nabla c) \cdot \mathbf{n} = 0 & \text{for } \mathbf{x} \in \partial\Omega, \end{cases} $$

with $D := \frac{D'}{R^2 K_d}$. The cell velocity equation remains unchanged.

### A.3.2 Reduction to a Darcy-Poisson problem

In the context of crawling cell migration, it is known that the friction forces associated with crawling are very large compared to the cytoskeleton viscosity [38, 39]. With that in mind, formally taking $\tilde{\eta} \to 0$ in the Stokes equatino reduces to a Darcy law

$$ \begin{cases} \mathbf{u}(t, \mathbf{x}) = -\nabla p(t, \mathbf{x}) - \mathbf{v}(t) & \text{in } \Omega, \\ p(t, \mathbf{x})\mathbf{n} = -\chi k_p(c)\mathbf{n} & \text{on } \partial\Omega. \end{cases} $$

*Remark 4* The Darcy law can also be obtained from a 3D Stokes problem with Dirichlet boundary conditions under an infinitesimal width assumption [35]. Such a formal justification is displayed in Appendix B.

Now, using the continuity equation on $\rho_\varepsilon$ combined with the Darcy law allows to write

$$ \begin{aligned} \varepsilon \partial_t p_\varepsilon &= \operatorname{div} (\rho_0(\nabla p_\varepsilon + \mathbf{v}) + \varepsilon p_\varepsilon(\nabla p_\varepsilon + \mathbf{v})) - (\rho_0 + \varepsilon p_\varepsilon) & \text{in } \Omega, \\ &= \operatorname{div} (\rho_0 \nabla p_\varepsilon + \varepsilon p_\varepsilon(\nabla p_\varepsilon + \mathbf{v})) - (\rho_0 + \varepsilon p_\varepsilon), \end{aligned} $$

since $\mathbf{v}$ does not depend on the space variable. Therefore, in the formal limit $\varepsilon \to 0$, the fluid density is constant, and the fluid pressure satisfies the following Poisson-Dirichlet problem:

$$ \begin{cases} -\Delta p(t, \mathbf{x}) = -1 & \text{in } \Omega, \\ p(t, \mathbf{x}) = -\chi k_p(c(t, \mathbf{x})) & \text{on } \partial\Omega. \end{cases} $$

The boundary condition accounts for polymerization at the edge of the cell, while the sink term $-1$ describes depolymerization in the cell body. Finally, the cell velocity now writes

$$ \mathbf{v}(t) = \gamma \int_{\partial\Omega} p(t, \mathbf{x})\mathbf{n} \, d\sigma(\mathbf{x}) = -\gamma\chi \int_{\partial\Omega} k_p(c(t, \mathbf{x}))\mathbf{n} \, d\sigma(\mathbf{x}), $$

with $\gamma := \frac{\xi}{\xi_1}$ the balance of friction. This expression shows that the cell displacement follows from the boundary coupling to the molecular dynamics. The new parameters are listed below.
<table>
  <caption>Non-dimensionalized parameters</caption>
  <tbody>
    <tr>
      <td>Viscosity</td>
      <td>η̃ := η / ξR<sup>2</sup></td>
    </tr>
    <tr>
      <td>Polymerization magnitude</td>
      <td>χ := k / ξR<sup>2</sup>K<sub>d</sub></td>
    </tr>
    <tr>
      <td>Diffusion coefficient</td>
      <td>D := D' / (L<sup>2</sup>V) = D' / (R<sup>2</sup>K<sub>d</sub>)</td>
    </tr>
    <tr>
      <td>Balance of frictions</td>
      <td>γ := ξ / ξ<sub>1</sub></td>
    </tr>
  </tbody>
</table>

Appendix B The 2D Darcy law originating in the
3D Stokes problem in small width

It is known that either for a confined cell [43] or in the context of cell crawling [35],
that the force balance equation reduces to a Darcy law with an effective coefficient.
Let us formally show the derivation in the spirit of [44] with the addition of the friction
force. Consider the Stokes problem on a 3D domain Ω × (0, ε) given by

$$
\begin{gather*}
-\eta \Delta \mathbf{u} + \nabla p = -\xi(\mathbf{u} + \mathbf{v}), \\
\mathbf{u}(t, x_1, x_2, 0) = \mathbf{u}(t, x_1, x_2, \varepsilon) = 0.
\end{gather*}
$$

The boundary condition for **u** on $x_3 \in \{0, \varepsilon\}$ means that there is no actin flow directly at the upper and lower parts of the membrane, As **ε** is very small, friction forces arise from actin flow rubbing against various molecular complexes all along the cell width. Using the change of variable $(x_1, x_2, x_3) \mapsto (x_1, x_2, y_3 = \frac{x_3}{\varepsilon})$ with $y_3 \in (0, 1)$ leads to

$$
\frac{\partial^2}{\partial x_3^2} = \frac{1}{\varepsilon^2} \frac{\partial^2}{\partial y_3^2} .
$$

Now, let us assume that we can write

$$
\begin{align*}
\mathbf{u}(t, x_1, x_2, y_3) &= \varepsilon^2 \mathbf{u}^0(t, x_1, x_2, y_3) + \varepsilon^3 \mathbf{u}^1(t, x_1, x_2, y_3) + \dots \\
p(t, x_1, x_2, y_3) &= p^0(t, x_1, x_2) + \varepsilon p^1(t, x_1, x_2, y_3) + \dots
\end{align*}
$$

Then, rewriting the Stokes problem for the *i*-th coordinate of the vectorial problem
(*i* ∈ {1, 2}), and taking the 0-th order term in *ε*, we obtain that

$$
-\eta \partial_{y_3}^2 u_i^0 + \partial_{x_i} p^0 = -\xi v_i(t),
$$
so that the friction related to **u** disappears. Now, since $p^0$ is independent of $y_3$, integrating this equation on $y_3$ on $(0, y)$ and $(y, 1)$ leads to

$$
\begin{align*}
\partial_{y_3} u_i^0 (y_3 = y) &= \partial_{y_3} u_i^0 (y_3 = 0) + \left( \frac{1}{\eta} \partial_{x_i} p^0 + \xi v_i \right) y, \\
&= \partial_{y_3} u_i^0 (y_3 = 1) - \left( \frac{1}{\eta} \partial_{x_i} p^0 + \xi v_i \right) (1-y),
\end{align*}
$$

so that

$$
\partial_{y_3} u_i^0 (y_3 = 0) = \partial_{y_3} u_i^0 (y_3 = 1) - \frac{1}{\eta} \partial_{x_i} p^0 - \xi v_i .
$$

Integrating again and using the Dirichlet boundary conditions on **u** yields

$$
u_i^0(y_3 = y) = \partial_{y_3} u_i^0(y_3 = 0)y + \left( \frac{1}{\eta} \partial_{x_i} p^0 + \xi v_i \right) \frac{y^2}{2}
$$

or equivalently

$$
u_i^0(y_3 = y) = -\partial_{y_3} u_i^0(y_3 = 0) + \partial_{y_3} u_i^0(y_3 = 0)y - \frac{1}{2} \left( \frac{1}{\eta} \partial_{x_i} p^0 + \xi v_i \right) (1 - y^2).
$$

By equality,

$$
\partial_{y_3} u_i^0 (y_3 = 0) = -\frac{1}{2\eta} \partial_{x_i} p^0 - \frac{\xi v_i}{2},
$$

$$
\partial_{y_3} u_i^0 (y_3 = 1) = \frac{1}{2\eta} \partial_{x_i} p^0 + \frac{\xi v_i}{2}.
$$

Finally,

$$
u_i^0(y_3 = y) = -\frac{1}{2} \left( \frac{1}{\eta} \partial_{x_i} p^0 + \xi v_i \right) y (1-y),
$$

so that integrating in $y_3$ leads to

$$
\int_{0}^{1} u_{i}^{0}(\cdot, y_{3}) \, dy_{3}=-\frac{1}{12} \left(\frac{1}{\eta} \partial_{x_{i}} p^{0}+\xi v_{i}\right),
$$

which concludes the derivation.

**Appendix C Numerical framework and cell polarity**
**indices**

In this part, we detail the finite volumes numerical framework used to perform the
numerical simulations. Then, we detail and illustrate the cell polarity indices used to
characterize the internal cell dynamics.
## C.1 Numerical framework

The radial discretization is given by {$r_j = j \Delta r, j = 1...N_r$}, where $r_{N_r+\frac{1}{2}} = 1$. Similarly, the angular discretization is given by {$\theta_k = k \Delta \theta, k \in 1...N_\theta$}, with $\theta_{N_\theta} = 2\pi$. For $j \in \{1...N_r\}$ and $k \in \{1...N_\theta\}$, the control volumes $W_{(j,k)} \subset \mathbb{R}^2$ are given by

$$W_{(j,k)} = \left(r_{j-\frac{1}{2}}, r_{j+\frac{1}{2}}\right) \times \left(\theta_{k-\frac{1}{2}}, \theta_{k+\frac{1}{2}}\right),$$

such that $|W_{(j,k)}| = r_j \Delta r \Delta \theta$. The center volume is given by $W_0 = B(0, r_{\frac{1}{2}})$ such that $|W_0| = \pi r_{\frac{1}{2}}^2$. We use a finite volumes scheme in space. More precisely, keeping the same notation for simplicity, the discretized values of the molecular concentration at time $t > 0$ are given by

$$\begin{cases} c_0(t) := \frac{1}{\pi r_{\frac{1}{2}}^2} \iint_{B(0,r_{\frac{1}{2}})} c(t,r,\theta)r \, dr \, d\theta, \\ c_{j,k}(t) := \frac{1}{|W_{(j,k)}|} \iint_{W_{(j,k)}} c(t,r,\theta)r \, dr \, d\theta, \quad j \in \{1...N_r\}, k \in \{1...N_\theta\}, \end{cases}$$

and their time dynamics are obtained using Equation (3). The advection term based on Equations (4)–(6) is approximated using a first-order upwind scheme. Finally, let {$t^n = n \Delta t, n = 0...N_f$} be the time discretization, with $T_f := t^{N_f}$ the final time of the simulation. Then, ($c_{j,k}^n := c_{j,k}(t^n)$) $1 \le i \le N_r, 1 \le j \le N_\theta$ and $c_0^n := c_0(t^n)$ are the fully discretized values of the solution at time $t^n$ obtained using an Euler scheme. Unless indicated otherwise, we take $\Delta t = 0.01$, and $N_\theta = 200$, so that $\Delta r = \Delta \theta = \frac{2\pi}{N_\theta} \approx 0.03$.

## C.2 Cell polarity indices

Two cell polarity indices are considered in this paper.

* the Relative Standard Deviation (RSD) or Coefficient of Variation (CV) is given by the perimeter-normalized ratio between the standard deviation of {$c_{N_r,k}^n$} $1 \le k \le N_y$ and its mean value $\bar{c}^n := \frac{1}{N_\theta} \sum_{k=1}^{N_\theta} c_{N_r,k}^n$. More precisely,

$$RSD^n := \frac{1}{2\pi\bar{c}^n} \sqrt{\frac{1}{N_y} \sum_{k=1}^{N_y} |c_{N_r,k}^n - \bar{c}^n|^2}.$$

* the normalized Mean Resultant Vector is given by

$$RV^n := \frac{1}{\sum_{i=1}^{N_\theta} c_{N_r,k}^n} \sum_{i=1}^{N_\theta} c_{N_r,k}^n \mathbf{e}_k$$

as defined in circular statistics [45]. Its norm, the Resultant Vector Length (RVL) provides another polarization index.

Figure C1 shows some boundary concentration profiles together with the corresponding RSD and RVL values. In Subfigure C1–(a), the boundary concentration is constant,
both indices are null and the cell is in a resting state. In Subfigure C1-(b), the boundary concentration shows a single pole, both the RSD and the RVL are nonzero, and the cell is in a traveling state. In Subfigure C1-(c), the boundary concentration shows two symmetric poles, such that the cell is in a resting state and the RVL is zero while the RSD is nonzero. Finally, in Subfigure C1-(d), the boundary concentration shows a steep polarization associated with a traveling state. However, both polarity indices are smaller than in Subfigure C1-(b), although nonzero.

<figure><img src="image_11.png" /><figcaption>**Fig. C1** Several boundary concentration profiles with the corresponding polarity indices. (a): Homogenous boundary concentration associated with a resting state, with null Relative Standard Deviation and Resultant Vector Length. (b): Unipolar boundary concentration given by $\theta \mapsto 1-\cos(\theta)$, associated with a traveling state, characterized by nonzero Relative Standard Deviation and Resultant Vector Length. (c): Boundary concentration given by $\theta \mapsto \cos(\theta)$ ($1_{[-\pi/2,\pi/2]} - 1_{[\pi/2,3\pi/2]}$), showing two opposite poles and associated with a resting state. The Relative Standard Deviation is nonzero while the Resultant Vector Length is zero. (d): Boundary concentration given by $\theta \mapsto 21_{[\pi-5\Delta\theta,\pi+5\Delta\theta]}$, with a steep polarization associated with a traveling state. Both Relative Standard Deviation and Resultant Vector Length are nonzero, although smaller than case (b).</figcaption></figure>
### C.3 Computation of the noise term

At each time step, the noise term is computed in cartesian coordinates over all nodes in both radial and angular interfaces to account for spatial correlations. Then, radial and angular components are derived only on the necessary points to account for radial or angular fluxes. More precisely, the set of interface nodes is given by

$$ \{\mathbf{x}_{(j+\frac{1}{2},k)}\}_{\substack{0 \le j \le N_r-1 \\ 1 \le k \le N_\theta}} \cup \{\mathbf{x}_{(j,k+\frac{1}{2})}\}_{\substack{1 \le j \le N_r \\ 0 \le k \le N_\theta-1}} $$

Denote by $A_{\Delta r, \Delta \theta}$ an ordered collection of this set, such that $|A_{\Delta r, \Delta \theta}| = 2N_r N_\theta - (N_r + N_\theta)$. Let us define $Z = (Z_\ell)_{1 \le \ell \le |A_{\Delta r, \Delta \theta}|}$ a Gaussian vector of size $|A_{\Delta r, \Delta \theta}|$ with covariance matrix computed from (46). Then, for $n > 0$, the approximation of the noise term at time $t^n$ and point $\mathbf{x}_i$ ($1 \le i \le |A_{\Delta r, \Delta \theta}|$) is given by

$$ d\mathbf{W}_n^{Q_\sigma}(\mathbf{x}_i) = \sqrt{\Delta t} Z_\ell. $$

### References

[1] Maiuri, P., Rupprecht, J.-F., Wieser, S., Ruprecht, V., Bénichou, O., Carpi, N., Coppey, M., De Beco, S., Gov, N., Heisenberg, C.-P., Lage Crespo, C., Lautenschlaeger, F., Le Berre, M., Lennon-Dumenil, A.-M., Raab, M., Thiam, H.-R., Piel, M., Sixt, M., Voituriez, R.: Actin flows mediate a universal coupling between cell speed and cell persistence. Cell 161(2), 374–86 (2015)

[2] Anon, E., Serra-Picamal, X., Hersen, P., Gauthier, N.C., Sheetz, M.P., Trepat, X., Ladoux, B.: Cell crawling mediates collective cell migration to close undamaged epithelial gaps. Proc Natl Acad Sci U S A 109(27), 10891–6 (2012)

[3] Bravo-Cordero, J.J., Hodgson, L., Condeelis, J.: Directed cell invasion and migration during metastasis. Curr Opin Cell Biol 24(2), 277–83 (2012)

[4] Lauffenburger, D.A., Horwitz, A.F.: Cell migration: a physically integrated molecular process. Cell 84(3), 359–69 (1996)

[5] Maiuri, P., Terriac, E., Paul-Gilloteaux, P., Vignaud, T., McNally, K., Onuffer, J., Thorn, K., Nguyen, P.A., Georgoulia, N., Soong, D., et al.: The First World Cell Race. Current Biology 22(17), 673–675 (2012)

[6] Ben Amar, M., Manyuhina, O., Napoli, G.: Cell motility: a viscous fingering analysis of active gels. The European Physical Journal Plus 126(2), 1–15 (2011)

[7] Blanch-Mercader, C., Casademunt, J.: Spontaneous motility of actin lamellar fragments. Physical Review Letters 110(7) (2013)

[8] Gallinato, O., Ohta, M., Poignard, C., Suzuki, T.: Free boundary problem for cell
protrusion formations: theoretical and numerical aspects. Journal of Mathemati-
cal Biology **75**(2), 263–307 (2017) https://doi.org/10.1007/s00285-016-1080-7

[9] Berlyand, L., Fuhrmann, J., Rybalko, V.: Bifurcation of traveling waves in a keller-segel type free boundary model of cell motility. Comm. Math. Sciences (2018)

[10] Barnhart, E., Lee, K.-C., Allen, G.M., Theriot, J.A., Mogilner, A.: Balance between cell-substrate adhesion and myosin contraction determines the frequency of motility initiation in fish keratocytes. Proc Natl Acad Sci U S A **112**(16), 5045–50 (2015)

[11] Etchegaray, C.: Mathematical and numerical modelling of cell migration. Theses, Université Paris-Saclay (November 2016)

[12] Rybalko, V., Berlyand, L.: Emergence of traveling waves and their stability in a free boundary model of cell motility. Transactions of the American Mathematical Society (376), 1799–1844 (2023)

[13] Lavi, I., Meunier, N., Voituriez, R., Casademunt, J.: Motility and morphodynamics of confined cells. Phys. Rev. E **101**, 022404 (2020) https://doi.org/10.1103/PhysRevE.101.022404

[14] Recho, P., Putelat, T., Truskinovsky, L.: Contraction-driven cell motility. Phys. Rev. Lett. **111**, 108102 (2013) https://doi.org/10.1103/PhysRevLett.111.108102

[15] Recho, P., Putelat, T., Truskinovsky, L.: Mechanics of motility initiation and motility arrest in crawling cells. Journal of the Mechanics and Physics of Solids **84**, 469–505 (2015) https://doi.org/10.1016/j.jmps.2015.08.006

[16] Cucchi, A., Mellet, A., Meunier, N.: Self polarization and traveling wave in a model for cell crawling migration. Discrete and Continuous Dynamical Systems - Series A **42**(5), 2381 (2022) https://doi.org/10.3934/dcds.2021194

[17] Muller, N., Piel, M., Calvez, V., Voituriez, R., Gonçalves-Sá, J., Guo, C.-L., Jiang, X., Murray, A., Meunier, N.: A predictive model for yeast cell polarization in pheromone gradients. PLoS Comput Biol **12**(4), 1004795 (2016)

[18] Cucchi, A., Mellet, A., Meunier, N.: A Cahn-Hilliard Model for Cell Motility. SIAM Journal on Mathematical Analysis **52**(4), 3843–3880 (2020) https://doi.org/10.1137/19M1267969

[19] Calvez, V., Hawkins, R.J., Meunier, N., Voituriez, R.: Analysis of a Nonlocal Model for Spontaneous Cell Polarization. SIAM Journal on Applied Mathematics **72**(2), 594–622 (2012)

[20] Lepoutre, T., Meunier, N., Muller, N.: Cell polarisation model : the 1D case.
Journal des Mathematiques Pures et Appliquees **101**(2), 152-171 (2014)

[21] Lepoutre, T., Meunier, N.: Analysis of a model of cell crawling migration. Communications in Mathematical Sciences **20**(6), 1589-1611 (2022) https://doi.org/10.4310/CMS.2022.v20.n6.a6

[22] Alazard, T., Magliocca, M., Meunier, N.: Mathematical study of a coupled incompressible darcy's free boundary problem with surface tension. Journal de l'École Polytechnique (2025)

[23] Meunier, N., Souplet, P.: Convergence, concentration and critical mass phenomenon in a model of cell motion with boundary signal production (2025)

[24] Hawkins, R.J., Bénichou, O., Piel, M., Voituriez, R.: Rebuilding cytoskeleton roads: active-transport-induced polarization of cells. Phys Rev E Stat Nonlin Soft Matter Phys **80**(4 Pt 1), 040903 (2009)

[25] Calvez, V., Meunier, N., Muller, N., Voituriez, R.: Numerical simulation of the dynamics of molecular markers involved in cell polarisation (2013)

[26] Muller, N.: Études mathématiques et numériques de problèmes non-linéaires et non-locaux issus de la biologie. PhD thesis, Université Paris Descartes (2013)

[27] Etchegaray, C., Meunier, N., Voituriez, R.: Analysis of a non-local and non-linear fokker-planck model for cell crawling migration. arXiv preprint arXiv:1701.06862 (2017)

[28] Peszat, S., Zabczyk, J.: Stochastic Partial Differential Equations with Lévy Noise: An Evolution Equation Approach vol. 113. Cambridge University Press, ??? (2007)

[29] Boulakia, M., Genadot, A., Thieullen, M.: Simulation of spdes for excitable media using finite elements. Journal of Scientific Computing **65**(1), 171-195 (2015)

[30] Callan-Jones, A.C., Voituriez, R.: Active gel model of amoeboid cell motility. New Journal of Physics **15**(2), 025022 (2013)

[31] Etchegaray, C., Meunier, N.: A stochastic model for protrusion activity. ESAIM: Proceedings and Surveys **62**, 56-67 (2018)

[32] De Gennes, P.G., Prost, J.: The Physics of Liquid Crystals. Oxford, ??? (1993)

[33] Kruse, K., Joanny, J.F., Jülicher, F., Prost, J., Sekimoto, K.: Generic theory of active polar gels: a paradigm for cytoskeletal dynamics. Eur Phys J E Soft Matter **16**(1), 5-16 (2005)

[34] Joanny, J.-F., Prost, J.: Active gels as a description of the actin-myosin cytoskeleton. HFSP J **3**(2), 94-104 (2009)
[35] Callan-Jones, A.C., Joanny, J.-F., Prost, J.: Viscous-fingering-like instability of cell fragments. Physical Review Letters **100**(25) (2008)

[36] Kruse, K., Joanny, J.F., Jülicher, F., Prost, J., Sekimoto, K.: Generic theory of active polar gels: a paradigm for cytoskeletal dynamics. Eur Phys J E Soft Matter **16**(1), 5–16 (2005)

[37] Jülicher, F., Kruse, K., Prost, J., Joanny, J.-F.: Active behavior of the cytoskeleton. Physics Reports **449**(1), 3–28 (2007)

[38] Prost, J., Jülicher, F., Joanny, J.: Active gel physics. Nature Physics **11**(2), 111–117 (2015)

[39] Kruse, K., Joanny, J.F., Jülicher, F., Prost, J.: Contractility and retrograde flow in lamellipodium motion. Phys Biol **3**(2), 130–7 (2006)

[40] Dang, I., Gorelik, R., Sousa-Blin, C., Derivery, E., Guérin, C., Linkner, J., Nemethova, M., Dumortier, J.G., Giger, F.A., Chipysheva, T.A., Ermilova, V.D., Vacher, S., Campanacci, V., Herrada, I., Planson, A.-G., Fetics, S., Henriot, V., David, V., Oguievetskaia, K., Lakisic, G., Pierre, F., Steffen, A., Boyreau, A., Peyriéras, N., Rottner, K., Zinn-Justin, S., Cherfils, J., Bièche, I., Alexandrova, A.Y., David, N.B., Small, J.V., Faix, J., Blanchoin, L., Gautreau, A.: Inhibitory signalling to the Arp2/3 complex steers cell migration. Nature **503**(7475), 281–4 (2013)

[41] Krause, M., Gautreau, A.: Steering cell migration: lamellipodium dynamics and the regulation of directional persistence. Nat Rev Mol Cell Biol **15**(9), 577–90 (2014)

[42] Oelz, D., Mogilner, A.: A drift-diffusion model for molecular motor transport in anisotropic filament bundles. Discrete and Continuous Dynamical Systems **36**(8), 4553–4567 (2016) https://doi.org/10.3934/dcds.2016.36.4553

[43] Hawkins, R.J., Piel, M., Faure-Andre, G., Lennon-Dumenil, A., Joanny, J., Prost, J., Voituriez, R.: Pushing off the walls: a mechanism of cell motility in confinement. Physical review letters **102**(5), 058103 (2009)

[44] Enault, S.: Modélisation de la propagation d'une tumeur en milieu faiblement compressible. PhD thesis, Lyon, Ecole normale supérieure (2010)

[45] Jammalamadaka, S.R., SenGupta, A.: Topics in Circular Statistics vol. 5. World Scientific, ??? (2001)