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

HAL
open science

Numerical solutions of a 2D fluid problem coupled to a
nonlinear non-local reaction-advection-diffusion problem
for cell crawling migration in a discoidal domain

Christèle Etchegaray, Nicolas Meunier

► To cite this version:

Christèle Etchegaray, Nicolas Meunier. Numerical solutions of a 2D fluid problem coupled to a non-linear non-local reaction-advection-diffusion problem for cell crawling migration in a discoidal domain. 2017. hal-01760518

HAL Id: hal-01760518

https://hal.science/hal-01760518

Preprint submitted on 11 Apr 2018

**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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical solutions of a 2D fluid problem coupled to a nonlinear non-local reaction-advection-diffusion problem for cell crawling migration in a discoidal domain

Christèle Etchegaray*

Nicolas Meunier†

## Abstract

In this work, we present a numerical scheme for the approximate solutions of a 2D crawling cell migration problem. The model, defined on a non-deformable discoidal domain, consists in a Darcy fluid problem coupled with a Poisson problem and a reaction-advection-diffusion problem. Moreover, the advection velocity depends on boundary values, making the problem nonlinear and non local.

For a discoidal domain, numerical solutions can be obtained using the finite volume method on the polar formulation of the model. Simulations show that different migration behaviours can be captured.

**Keywords:** 2D cell migration, Darcy fluid, nonlinear non-local reaction-advection-diffusion problem, finite volume method.

## 1 Introduction

Cell migration ensures fundamental functions in the body (embryogenesis, immune system), but is also involved in the development of pathologies such as tumor metastasis, arising large research efforts (Bravo-Cordero et al., 2012). However, the responsible intracellular mechanisms involve multiscale interaction in time and space, so that modelling cell migration is challenging and produces interesting problems to study.

For 2D cells crawling on a surface, the motion is friction-based and results from the activity of the so-called **actin cytoskeleton**, that is a dynamics mesh of actin filaments. They are polar: they polymerize at one end and depolymerize at the other end, under the molecular regulation of signalling loops. Overall, the actin mesh can be approximated by an active fluid in our setting (Kruse et al., 2005; Joanny and Prost, 2009).

The mesh grows preferentially at the cell membrane, and shrinks inside the cell body, generating inward actin flows from the membrane. Its mechanical connection to the adhesive substrate generates the friction forces responsible for the cell's displacement.

Modelling this process is a difficult task, because of the large time and space scales, and also because of the large number of effectors of the dynamics. Following key physical ideas of Blanch-Mercader and Casademunt (2013); Maiuri et al. (2015), we proposed in Etchegaray

<footnote>*Institut de Mathématiques de Toulouse, Université Paul Sabatier, 118, route de Narbonne F-31062 Toulouse Cedex 9 (christele.etchegaray@math.univ-toulouse.fr)</footnote>

<footnote>†MAP5, CNRS UMR 8145, Université Paris Descartes, 45 rue des Saints Pères 75006 Paris, France. (nicolas.meunier@parisdescartes.fr)</footnote>
(2016); Etchegaray et al. (2017b,a) a minimal multiscale model for 2D crawling migration, that we recall now.

The cell domain is a non deformable disc $\Omega$, and the problem is formulated in the cell's frame of reference (the domain does not move). First, the actin cytoskeleton is approximated by a Darcy fluid (Blanch-Mercader and Casademunt, 2013), and its dynamics in a crawling situation is modelled by a Poisson problem on the fluid pressure (Etchegaray et al., 2017b). More precisely, for $\mathbf{u}: \mathbb{R}_+ \times \Omega \to \mathbb{R}^2$ the fluid velocity, $p: \mathbb{R}_+ \times \Omega \to \mathbb{R}$ its pressure, we have

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

for $k_d$ the actin depolymerization rate, and $k_p: \mathbb{R}_+ \times \Omega \to \mathbb{R}$ the polymerization function rate at the boundary. The fluid velocity writes

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

with for $\gamma \in \mathbb{R}_+$ the domain velocity

$$ \mathbf{v}(t) = \gamma \int_{\Omega} \nabla p(t, \mathbf{x}) \, d\mathbf{x} = \gamma \int_{\partial\Omega} p(t, \mathbf{x}) \mathbf{n}_{\mathbf{x}} \, d\mathbf{x}, \qquad (3) $$

for $\mathbf{n}_{\mathbf{x}}$ the unit normal vector to $\partial\Omega$ at point $\mathbf{x} \in \partial\Omega_0$.

These equations show that the actin polymerization at the boundary and depolymerization inside the domain may drive a pressure gradient leading to the cell motion. Note that the dynamics therefore arises from the activity at the boundary, and that the cell velocity is nonlocal.

Now, the interaction with the molecular scale consists in studying the dynamics of a molecular inhibitor to polymerization. The molecules diffuse freely inside the cell in an inactive form. They may bind actin filaments and be carried by their flow. Moreover, if activated at the cell membrane, they become able to locally inhibit actin polymerization. Write $c: \mathbb{R}_+ \times \Omega \to \mathbb{R}$ the concentration in an inactive form, and $\mu: \mathbb{R}_+ \times \partial\Omega \to \mathbb{R}$ the activated concentration at the boundary. Then,

→ the inactive molecules follow an advection-diffusion dynamics with advection velocity $\mathbf{u}$,

→ there is an exchange dynamics on $\partial\Omega$ between active and inactive forms,

→ $k_p$ is a decreasing function of $\mu$.

For $D$ the diffusion coefficient, and $k_{\text{on/off}}$ the activation/desactivation rates, the corresponding problem writes

$$ \begin{cases} \partial_t c(t, \mathbf{x}) + \text{div} \left( c(t, \mathbf{x}) \mathbf{u}(t, \mathbf{x}) - D \nabla c(t, \mathbf{x}) \right) = 0 & \text{in } \Omega, \\ (D \nabla c(t, \mathbf{x}) - c(t, \mathbf{x}) \mathbf{u}(t, \mathbf{x})) \cdot \mathbf{n}_{\mathbf{x}} = -k_{\text{on}} c(t, \mathbf{x}) + k_{\text{off}} \mu(t, \mathbf{x}) & \text{on } \partial\Omega, \\ \frac{\partial}{\partial t} \mu(t, \mathbf{x}) = k_{\text{on}} c(t, \mathbf{x}) - k_{\text{off}} \mu(t, \mathbf{x}) & \text{on } \partial\Omega. \end{cases} \qquad (4) $$

Note that the boundary condition ensures mass conservation:

$$ \frac{d}{dt} \left( \int_{\Omega} c(t, \mathbf{x}) \mathrm{d}\mathbf{x} + \int_{\partial\Omega} \mu(t, \mathbf{x}) \mathrm{d}\mathbf{x} \right) = 0. \qquad (5) $$

Some remarks can be made to highlight the difficulties in the analysis of the model. The fluid velocity rewrites

$$ \mathbf{u}(t, \mathbf{x}) = -\nabla p(t, \mathbf{x}) - \gamma \int_{\partial\Omega} k_p(t, \mathbf{x}) \mathbf{n}_{\mathbf{x}} \, d\mathbf{x}. \qquad (6) $$
We notice that this expression depends on the concentration in activated molecules $\mu$ from the pressure boundary term $k_p$, so that the reaction-advection-diffusion problem is nonlinear. Moreover, the integral term makes it also non-local.

The corresponding 1D model in a special case without activation at the boundary and for a linear $k_p$ writes as a nonlinear non-local advection-diffusion problem, and has been studied in Etchegaray et al. (2017b). The analysis shows the existence of different asymptotic solutions, describing both motile and non motile behaviours. Moreover, for a subcritical mass of molecules, the global weak existence of solutions is established, along with their convergence to a non motile gaussian profile at explicit rate. Finally, under conditions on the initial condition, it is shown that solutions blow up in finite time.

Since the model carries non trivial behaviours, computing numerical solutions is of particular interest. In the following, we develop a finite volume method for the polar formulation of the problem (1)-(2)-(4), since the domain $\Omega$ is a disc.

## 2 Numerical method

We introduce now the finite volume discretization of the model. We assume that the cytokeleton domain is an annulus $\Omega = B(0, R) \setminus B(0, R_{min}) \subset \mathbb{R}^2$, where $B(0, R_{min})$ accounts for the nucleus. We consider a zero-flux boundary condition for the molecular concentration on $C(0, R_{min})$, the circle of center $(0, 0)$ and radius $R_{min}$. As a consequence, it will be natural to study the problem in polar coordinates.

### 2.1 Polar formulation

#### Reaction-advection-diffusion problem

Let $\mathbf{x} = (r \cos(\theta), r \sin(\theta)) \in \Omega$, and $\tilde{c}$ (similarly $\tilde{\mu}$) the polar function such that $\frac{1}{r}\tilde{c}(t, r, \theta) = c(t, \mathbf{x})$ with $(r, \theta) \in [R_{min}, R] \times \mathbb{R}/2\pi\mathbb{Z}$. Then, the problem on the molecular specie writes

$$ \partial_t \tilde{c}(t, r, \theta) = \partial_r \left( D r \partial_r \left( \frac{\tilde{c}(t, r, \theta)}{r} \right) - \tilde{c}(t, r, \theta) \mathbf{u}_r(t, r, \theta) \right) + \partial_\theta \left( \frac{1}{r^2} (D \partial_\theta \tilde{c}(t, r, \theta) - \tilde{c}(t, r, \theta) \mathbf{u}_\theta(t, r, \theta)) \right) \text{ in } \Omega, \quad (7) $$

$$ k_{\text{off}}\tilde{\mu}(t, R, \theta) - k_{\text{on}}\tilde{c}(t, R, \theta) = D R \partial_r \left( \frac{\tilde{c}(t, R, \theta)}{R} \right) - \tilde{c}(t, R, \theta) \mathbf{u}_r(t, R, \theta), \quad \text{on } C(0, R), \quad (8) $$

$$ 0 = D R_{\min} \partial_r \left( \frac{\tilde{c}(t, R_{\min}, \theta)}{R_{\min}} \right) - \tilde{c}(t, R_{\min}, \theta) \mathbf{u}_r(t, R_{\min}, \theta), \quad \text{on } C(0, R_{\min}), \quad (9) $$

$$ \partial_t \tilde{\mu}(t, R, \theta) = k_{\text{on}} \tilde{c}(t, R, \theta) - k_{\text{off}} \tilde{\mu}(t, R, \theta), \quad \text{on } C(0, R). \quad (10) $$

#### Poisson problem on $p$

For the polymerization function $k_p$, let us choose a simple form, that is

$$ k_p(t, \mathbf{x}) = [1 - \delta\mu(t, \mathbf{x})]_+, \quad (11) $$

with $\delta > 0$ and $\mathbf{x} \in \partial\Omega$.
The pressure $p$ is solution of the following problem:

$$-\Delta p(t, \mathbf{x}) = -k_d, \text{ in } \Omega, \qquad (12)$$

$$p(t, \mathbf{x}) = [1 - \delta\mu(t, \mathbf{x})]_{+}, \text{ on } C(0, R), \qquad (13)$$

$$p(t, \mathbf{x}) = 0, \text{ on } C(0, R_{\min}), \qquad (14)$$

where the pressure condition on $C(0, R_{\min})$ is arbitrary, and fix the pressure values. Let us consider these equations in polar coordinates with $\frac{1}{r}\tilde{p}(t, r, \theta) = p(t, \mathbf{x})$ for $(r, \theta) \in [R_{\min}, R] \times \mathbb{R}/2\pi\mathbb{Z}$. We have

$$-\partial_r \left( r \partial_r \left( \frac{\tilde{p}(r, \theta)}{r} \right) \right) - \frac{1}{r^2} \partial_{\theta\theta} \tilde{p}(r, \theta) = -k_d r, \quad \text{in } \Omega, \qquad (15)$$

$$\tilde{p}(t, R, \theta) = R \left[ 1 - \delta \frac{\tilde{\mu}(t, R, \theta)}{R} \right]_{+} \quad \text{on } C(0, R), \qquad (16)$$

$$\tilde{p}(t, R_{\min}, \theta) = 0 \quad \text{on } C(0, R_{\min}). \qquad (17)$$

## 2.2 Discretization

Let $t^n = n \Delta t$ be the time discretization, and $\{r_j = R_{\min} + (j - \frac{1}{2}) \Delta r, j \in \{1, ..., N_r\}\}$ the space discretization of the bounded interval $[R_{\min}, R]$, such that $r_{N_r + \frac{1}{2}} = R$, (therefore $N_R = \frac{R - R_{\min}}{\Delta r}$). Similarly, $\{\theta_k = k \Delta \theta, k \in \{1, ..., N_{\theta}\}\}$ is the space discretization of the periodic interval $\mathbb{R}/2\pi\mathbb{Z}$. We introduce the control volumes $W_{(j,k)} \subset \mathbb{R}^2$ and $V_k \subset \mathbb{R}/2\pi\mathbb{Z}$ with

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

Let $\tilde{c}^n_{(j,k)}$ (resp. $\tilde{\mu}^n_k$) be the approximated value of the exact solution $\tilde{c}(t^n, r_j, \theta_k)$ (resp. $\tilde{\mu}(t^n, \theta_k)$), and $\tilde{p}^n_{(j,k)}$ be the approximated value of the exact solution $\tilde{p}(t^n, r_j, \theta_k)$:

$$\tilde{c}^n_{(j,k)} \approx \frac{1}{\Delta r \Delta \theta} \iint_{W_{(j,k)}} \tilde{c}(t^n, r, \theta) dr d\theta, \\ \tilde{\mu}^n_k \approx \frac{1}{\Delta \theta} \int_{V_k} \tilde{\mu}(t^n, \theta) d\theta, \\ \tilde{p}^n_{(j,k)} \approx \frac{1}{\Delta r \Delta \theta} \iint_{W_{(j,k)}} \tilde{p}(t^n, r, \theta) dr d\theta.$$

Moreover, we write $\mathbf{u}^n$ and $\mathbf{v}^n$ the corresponding discretized velocity functions at time $t^n$.

The resolution is made as follows: for $n \ge 0$, knowing $(\tilde{c}^n, \tilde{\mu}^n)$ allows to compute $\tilde{p}^n$, then $(\mathbf{u}^n, \mathbf{v}^n)$. Finally, $(\tilde{c}^{n+1}, \tilde{\mu}^{n+1})$ is computed using $\mathbf{u}^n$, and so on.

### Problem on $\tilde{p}$

The pressure problem is time-dependent because of the Dirichlet boundary condition. Therefore, it is solved explicitly in time. Write $\mathcal{F}$ for the numerical flux. Then, we have the following scheme for equation (15): for $(j, k) \in \{1, ..., N_r\} \times \{1, ..., N_{\theta}\}$,

$$- \left( \frac{\mathcal{F}_{(j+\frac{1}{2},k)} - \mathcal{F}_{(j-\frac{1}{2},k)}}{\Delta r} + \frac{\mathcal{F}_{(j,k+\frac{1}{2})} - \mathcal{F}_{(j,k-\frac{1}{2})}}{\Delta \theta} \right) = -k_d r_j.$$
The finite volume numerical fluxes are defined by

$$
\begin{align*}
\mathcal{F}_{(j+\frac{1}{2},k)} &= r_{j+\frac{1}{2}} \frac{\frac{\tilde{p}_{(j+1,k)}}{r_{j+1}} - \frac{\tilde{p}_{(j,k)}}{r_j}}{\Delta r}, &
\mathcal{F}_{(j-\frac{1}{2},k)} &= r_{j-\frac{1}{2}} \frac{\frac{\tilde{p}_{(j,k)}}{r_j} - \frac{\tilde{p}_{(j-1,k)}}{r_{j-1}}}{\Delta r}, \\
\mathcal{F}_{(j,k+\frac{1}{2})} &= \frac{1}{r_j^2} \frac{\tilde{p}_{(j,k+1)} - \tilde{p}_{(j,k)}}{\Delta \theta}, &
\mathcal{F}_{(j,k-\frac{1}{2})} &= \frac{1}{r_j^2} \frac{\tilde{p}_{(j,k)} - \tilde{p}_{(j,k-1)}}{\Delta \theta}.
\end{align*}
$$

The Dirichlet boundary conditions (16)-(17) are imposed using ghost values $\tilde{p}_{(0,k)}$ and $\tilde{p}^n_{(N_r+1,k)}$.
For $k \in \{1, \dots, N_\theta\}$,

$$
\mathcal{F}_{(\frac{1}{2},k)} = \frac{r_{\frac{1}{2}}}{r_1} \frac{\tilde{p}_{(1,k)}}{\Delta r},
$$

since $\tilde{p}_{(0,k)} = 0$. Now, $\tilde{p}^n_{(N_r+1,k)} = r_{N_r} \left[ 1 - \delta \frac{\tilde{\mu}^n_{(N_r,k)}}{r_{N_r}} \right]_+$, the corresponding flux writes

$$
\mathcal{F}_{(N_r+\frac{1}{2},k)} = r_{N_r+\frac{1}{2}} \frac{\frac{r_{N_r}}{r_{N_r+1}} \left(1 - \delta \frac{\tilde{\mu}^n_{(N_r,k)}}{r_{N_r}}\right) - \frac{\tilde{p}_{(N_r,k)}}{r_{N_r}}}{\Delta r}.
$$

Therefore, the term $\frac{r_{N_r} r_{N_r+\frac{1}{2}}}{r_{N_r+1} \Delta r} \left[1 - \delta \frac{\tilde{\mu}^n_{(N_r,k)}}{r_{N_r}}\right]_+$ will be included in the right hand side of the
matricial problem.

Similarly, the periodic conditions impose for $j \in \{1, \dots, N_r\}$,

$$
\mathcal{F}_{(j,N_{\theta}+\frac{1}{2})} = \mathcal{F}_{(j,\frac{1}{2})} = \frac{1}{r_j^2} \frac{\tilde{p}_{(j,1)} - \tilde{p}_{(j,N_{\theta})}}{\Delta\theta}.
$$

As a consequence, we write the corresponding system, where the terms in bold account for
the boundary conditions. Note also that each equation is written for $k \in \{1, \dots, N_\theta\}$ with the
convention that if $k = 1$, then $k - 1 = N_\theta$ and $k - \frac{1}{2} = \frac{1}{2}$. Similarly, if $k = N_\theta$, then $k + 1 = 1$
and $k + \frac{1}{2} = \frac{1}{2}$.

For $j = 1$, we write

$$
\frac{1}{\Delta r^2} \left[ r_{\frac{3}{2}} \left( \frac{\tilde{p}_{(2,k)}}{r_2} - \frac{\tilde{p}_{(1,k)}}{r_1} \right) - \frac{\mathbf{r}_{\frac{1}{2}}}{\mathbf{r}_1} \tilde{\mathbf{p}}_{(1,k)} \right] + \frac{\tilde{p}_{(1,k+1)} - 2\tilde{p}_{(1,k)} + \tilde{p}_{(1,k-1)}}{r_1^2 \Delta \theta^2} = k_d r_1. \quad (18)
$$

or equivalently

$$
\frac{1}{\Delta r^2} \left[ \frac{r_{\frac{1}{2}} + r_{\frac{3}{2}}}{r_1} \tilde{p}_{(1,k)} - \frac{r_{\frac{3}{2}}}{r_2} \tilde{p}_{(2,k)} \right] + \frac{-\tilde{p}_{(1,k-1)} + 2\tilde{p}_{(1,k)} - \tilde{p}_{(1,k+1)}}{r_1^2 \Delta \theta^2} = -k_d r_1. \quad (19)
$$

For $j \in \{2, \dots, N_\theta - 1\}$, we write for $k \in \{1, \dots, N_\theta\}$,

$$
\frac{1}{\Delta r^2} \left[ r_{j+\frac{1}{2}} \left( \frac{\tilde{p}_{(j+1,k)}}{r_{j+1}} - \frac{\tilde{p}_{(j,k)}}{r_j} \right) - r_{j-\frac{1}{2}} \left( \frac{\tilde{p}_{(j,k)}}{r_j} - \frac{\tilde{p}_{(j-1,k)}}{r_{j-1}} \right) \right] + \frac{\tilde{p}_{(j,k+1)} - 2\tilde{p}_{(j,k)} + \tilde{p}_{(j,k-1)}}{r_j^2 \Delta \theta^2} = k_d r_j. \quad (20)
$$

or equivalently

$$
\frac{1}{\Delta r^2} \left[ -\frac{r_{j-\frac{1}{2}}}{r_{j-1}} \tilde{p}_{(j-1,k)} + \frac{r_{j+\frac{1}{2}} + r_{j-\frac{1}{2}}}{r_j} \tilde{p}_{(j,k)} - \frac{r_{j+\frac{1}{2}}}{r_{j+1}} \tilde{p}_{(j+1,k)} \right] + \frac{-\tilde{p}_{(j,k-1)} + 2\tilde{p}_{(j,k)} - \tilde{p}_{(j,k+1)}}{r_j^2 \Delta \theta^2} = -k_d r_j. \quad (21)
$$
Finally, for $j = N_r$, we write for $k \in \{1, ..., N_\theta\}$,

$$ \frac{1}{\Delta r^2} \left[ -\frac{\mathbf{r}_{N_r+\frac{1}{2}}}{\mathbf{r}_{N_r}} \tilde{\mathbf{p}}(\mathbf{N_r},k) - r_{N_r-\frac{1}{2}} \left( \frac{\tilde{p}_{(N_r,k)}}{r_{N_r}} - \frac{\tilde{p}_{(N_r-1,k)}}{r_{N_r-1}} \right) \right] + \frac{\tilde{p}_{(N_r,k+1)} - 2\tilde{p}_{(N_r,k)} + \tilde{p}_{(N_r,k-1)}}{r_j^2 \Delta \theta^2} = k_d r_{N_r} \\ - \frac{r_{N_r} r_{N_r+\frac{1}{2}}}{r_{N_r+1} \Delta r} \left[ 1 - \delta \frac{\tilde{\mu}_{N_r,k}^n}{r_{N_r}} \right]_+ . \qquad (22) $$

or equivalently

$$ \frac{1}{\Delta r^2} \left[ -\frac{r_{N_r-\frac{1}{2}}}{r_{N_r-1}} \tilde{p}_{(N_r-1,k)} + \frac{r_{N_r-\frac{1}{2}} + r_{N_r+\frac{1}{2}}}{r_{N_r}} \tilde{p}_{(N_r,k)} \right] + \frac{-\tilde{p}_{(N_r,k-1)} + 2\tilde{p}_{(N_r,k)} - \tilde{p}_{(N_r,k+1)}}{r_j^2 \Delta \theta^2} = -k_d r_{N_r} \\ + \frac{r_{N_r} r_{N_r+\frac{1}{2}}}{r_{N_r+1} \Delta r} \left[ 1 - \delta \frac{\tilde{\mu}_{N_r,k}^n}{r_{N_r}} \right], \qquad (23) $$

Let us now write the corresponding matricial problem. We define the column vector $\mathcal{P}$ by $\mathcal{P}(k + (j - 1)N_\theta) = \tilde{p}_{(j,k)}$ with $(j, k) \in \{1, ..., N_r\} \times \{1, ..., N_\theta\}$:

$$ \mathcal{P} = (\tilde{p}_{(1,1)} \cdots \tilde{p}_{(1,N_\theta)} \tilde{p}_{(2,1)} \cdots \tilde{p}_{(2,N_\theta)} \cdots \tilde{p}_{(N_r,N_\theta)})^T $$

For $\Delta r = \Delta \theta$ the stiffness matrix $A_p$ is defined by

$$ A_p = \begin{pmatrix} \frac{r_{1/2}+r_{1+1/2}}{r_1} I_{N_\theta} & -\frac{r_{1+1/2}}{r_2} I_{N_\theta} & & & & \\ & \ddots & & & & \\ & & \ddots & & & \\ & & & \frac{r_{j-1/2}}{r_{j-1}} I_{N_\theta} & \frac{r_{j-1/2}+r_{j+1/2}}{r_j} I_{N_\theta} & -\frac{r_{j+1/2}}{r_{j+1}} I_{N_\theta} & \\ & & & & \ddots & \\ & & & & & \ddots & \\ & & & & & & \frac{r_{N_r-1/2}}{r_{N_r-1}} I_{N_\theta} & \frac{r_{N_r-1/2}+r_{N_r+1/2}}{r_{N_r}} I_{N_\theta} \\ & & & & & & & \frac{1}{r_1^2} A & \\ & & & & & & & \frac{1}{r_2^2} A & \\ & & & & & & & & \ddots & \\ & & & & & & & & & \frac{1}{r_{N_r-1}^2} A & \\ & & & & & & & & & & \frac{1}{r_{N_r}^2} A \end{pmatrix}, \quad (24) $$

where the second matrix accounts for the angular diffusion, with $A \in M_{N_\theta}(\mathbb{R})$ the classical diffusion matrix with periodic flux boundary conditions:

$$ A = \begin{pmatrix} 2 & -1 & & & -1 \\ -1 & 2 & \ddots & & \\ & \ddots & \ddots & & \\ & & 2 & -1 & \\ -1 & & -1 & 2 & \end{pmatrix}. \qquad (25) $$
The right-hand side in (15) and the flux boundary condition (16) on C(0, $R_{max}$) imposes this right hand side column vector of length $N_r$ $N_\theta$:

$$ \mathcal{R}_p^n = -k_d \begin{pmatrix} (r_1)_k \\ \vdots \\ (r_j)_k \\ \vdots \\ (r_{N_r})_k \end{pmatrix} + \frac{r_{N_r} r_{N_r+\frac{1}{2}}}{r_{N_r+1} \Delta r} \begin{pmatrix} 0 \\ \vdots \\ 0 \\ \left[1 - \delta^{\frac{\tilde{\mu}_{N_r,k}^n}{r_{N_r}}}\right]_+ \end{pmatrix}. $$

We use a standard numerical method to invert the symmetric positive definite matrix $\frac{1}{\Delta r^2} \mathcal{A}_p$ and then resolve at each time step

$$ \mathcal{P} = \left( \frac{1}{\Delta r^2} \mathcal{A}_p \right)^{-1} \mathcal{R}_p^n. $$

## Equations for u and v

The velocities **u** and **v** depend on the pressure $p$, which is obtained through a time explicit scheme. Therefore, they also depend explicitly in time on the concentration.

The equation on **v** in polar coordinates writes

$$ \mathbf{v}(t) = \gamma \int_{0}^{2\pi} \left[ 1 - \delta \frac{\tilde{\mu}(t, R, \theta)}{R} \right]_{+} \mathbf{n} d\theta. \quad (26) $$

We compute numerically the velocity in cartesian coordinates $\mathbf{v}_{\text{cart}}^n := (v_x^n, v_y^n)^T$:

$$ v_x^n = \gamma \Delta\theta \sum_{k=1}^{N_\theta} \left[ 1 - \delta \frac{\tilde{\mu}_{(N_r,k)}^n}{R} \right]_+ \cos(\theta_k), \quad (27) $$

$$ v_y^n = \gamma \Delta\theta \sum_{k=1}^{N_\theta} \left[ 1 - \delta \frac{\tilde{\mu}_{(N_r,k)}^n}{R} \right]_+ \sin(\theta_k). \quad (28) $$

Then, a polar change of coordinates leads to $\mathbf{v}^n := (v_r^n, v_\theta^n)^T$.

For the fluid velocity, we have

$$ \mathbf{u}(t, \mathbf{x}) = -\nabla p(t, \mathbf{x}) - \mathbf{v}(t), \quad (29) $$

that rewrites

$$ \mathbf{u}(t, r, \theta) = - \left( \partial_r \left( \frac{\tilde{p}(t, r, \theta)}{r} \right) + v_r \right) \mathbf{e}_r - \left( \frac{1}{r} \partial_\theta \tilde{p}(t, r, \theta) + v_\theta \right) \mathbf{e}_\theta, \quad (30) $$

since $u_r = u(t, \mathbf{x}) \cdot \mathbf{e}_r$ and $u_\theta = r u(t, \mathbf{x}) \cdot \mathbf{e}_\theta$. We define at time $t^n$

$$ u_{(j+\frac{1}{2},k)}^n = -\frac{\frac{\tilde{p}_{(j+1,k)}}{r_{j+1}} - \frac{\tilde{p}_{(j,k)}}{r_j}}{\Delta r} - v_r^n, $$

$$ u_{(j-\frac{1}{2},k)}^n = -\frac{\frac{\tilde{p}_{(j,k)}}{r_j} - \frac{\tilde{p}_{(j-1,k)}}{r_{j-1}}}{\Delta r} - v_r^n, $$

$$ u_{(j,k+\frac{1}{2})}^n = -\frac{1}{r_j} \frac{\tilde{p}_{(j,k+1)} - \tilde{p}_{(j,k)}}{\Delta\theta} - v_\theta^n, $$

$$ u_{(j,k-\frac{1}{2})}^n = -\frac{1}{r_j} \frac{\tilde{p}_{(j,k)} - \tilde{p}_{(j,k-1)}}{\Delta\theta} - v_\theta^n. $$
### 2.2.1 Problem for $\tilde{c}$ and $\tilde{\mu}$

For simplicity, we call again $\mathcal{F}$ the numerical fluxes. We can write the following scheme for equation (7): for $(j, k) \in \{1, ..., N_r\} \times \{1, ..., N_\theta\}$,

$$ \frac{\tilde{c}_{(j,k)}^{n+1} - \tilde{c}_{(j,k)}^{n}}{\Delta t} = \frac{\mathcal{F}_{(j+\frac{1}{2},k)} - \mathcal{F}_{(j-\frac{1}{2},k)}}{\Delta r} + \frac{\mathcal{F}_{(j,k+\frac{1}{2})} - \mathcal{F}_{(j,k-\frac{1}{2})}}{\Delta \theta}. $$

We define now the numerical fluxes. The diffusion part is implicit, so that no CFL condition is needed (Allaire, 2005), while the advection is explicit due to the nonlinearity in the expression of $v$.

We have:

$$ \mathcal{F}_{(j+\frac{1}{2},k)} = D r_{j+\frac{1}{2}} \frac{\frac{\tilde{c}_{(j+1,k)}^{n+1}}{r_{j+1}} - \frac{\tilde{c}_{(j,k)}^{n+1}}{r_j}}{\Delta r} - A^{up} \left( u_{(j+\frac{1}{2},k)}^{n}, \tilde{c}_{(j,k)}^{n}, \tilde{c}_{(j+1,k)}^{n} \right), $$

$$ \mathcal{F}_{(j-\frac{1}{2},k)} = D r_{j-\frac{1}{2}} \frac{\frac{\tilde{c}_{(j,k)}^{n+1}}{r_j} - \frac{\tilde{c}_{(j-1,k)}^{n+1}}{r_{j-1}}}{\Delta r} - A^{up} \left( u_{(j-\frac{1}{2},k)}^{n}, \tilde{c}_{(j-1,k)}^{n}, \tilde{c}_{(j,k)}^{n} \right), $$

$$ \mathcal{F}_{(j,k+\frac{1}{2})} = \frac{1}{r_j^2} \left( D \frac{\tilde{c}_{(j,k+1)}^{n+1} - \tilde{c}_{(j,k)}^{n+1}}{\Delta \theta} - A^{up} \left( u_{(j,k+\frac{1}{2})}^{n}, \tilde{c}_{(j,k)}^{n}, \tilde{c}_{(j,k+1)}^{n} \right) \right), $$

$$ \mathcal{F}_{(j,k+\frac{1}{2})} = \frac{1}{r_j^2} \left( D \frac{\tilde{c}_{(j,k)}^{n+1} - \tilde{c}_{(j,k-1)}^{n+1}}{\Delta \theta} - A^{up} \left( u_{(j,k-\frac{1}{2})}^{n}, \tilde{c}_{(j,k-1)}^{n}, \tilde{c}_{(j,k)}^{n} \right) \right), $$

where $A^{up}$ is the advection term expressed by

$$ A^{up}(\mathbf{u}, x_{-}, x_{+}) = \begin{cases} \mathbf{u} x_{-}, & \text{if } u > 0, \\ \mathbf{u} x_{+}, & \text{if } u < 0. \end{cases} \qquad (31) $$

The external boundary condition (8) yields

$$ \mathcal{F}_{(N_r+\frac{1}{2},k)} = k_{\text{off}} \mu_k^{n+1} - k_{\text{on}} \tilde{c}_{N_r,k}^{n+1} $$

for $k \in \{1, ..., N_\theta\}$. The zero flux boundary condition (9) imposes that $\mathcal{F}_{(\frac{1}{2},k)} = 0$ for $k \in \{1, ..., N_\theta\}$. Similarly, the periodic conditions impose for $j \in \{1, ..., N_r\}$

$$ \mathcal{F}_{(j,N_{\theta}+\frac{1}{2})} = \mathcal{F}_{(j,\frac{1}{2})} = \frac{1}{r_j^2} \left( D \frac{\tilde{c}_{(j,1)}^{n+1} - \tilde{c}_{(j,N_\theta)}^{n+1}}{\Delta \theta} - A^{up} \left( u_{(j,\frac{1}{2})}^{n}, \tilde{c}_{(j,N_\theta)}^{n}, \tilde{c}_{(j,1)}^{n} \right) \right). $$

We write the corresponding scheme and group the implicit (resp. explicit) terms on the left-hand-side (resp. right-hand-side). Note also that each equation is written for $k \in \{1, ..., N_\theta\}$ with the convention that if $k = 1$, then $k - 1 = N_\theta$ and $k - \frac{1}{2} = \frac{1}{2}$. Similarly, if $k = N_\theta$, then $k + 1 = 1$ and $k + \frac{1}{2} = \frac{1}{2}$.

For $j = 1$, we have

$$ \left( 1 + \frac{D \Delta t}{\Delta r^2} \frac{r_{1+\frac{1}{2}}}{r_1} \right) \tilde{c}_{(1,k)}^{n+1} - \frac{D \Delta t}{\Delta r^2} \frac{r_{1+\frac{1}{2}}}{r_2} \tilde{c}_{(2,k)}^{n+1} + \frac{D \Delta t}{r_1^2 \Delta \theta^2} \left( -\tilde{c}_{(1,k-1)}^{n+1} + 2\tilde{c}_{(1,k)}^{n+1} - \tilde{c}_{(1,k+1)}^{n+1} \right) \\ = \tilde{c}_{(1,k)}^{n} - \frac{\Delta t}{\Delta r} A^{up} \left( u_{(1+\frac{1}{2},k)}^{n}, \tilde{c}_{(1,k)}^{n}, \tilde{c}_{(2,k)}^{n} \right) \quad (32) \\ + \frac{\Delta t}{r_1^2 \Delta \theta} \left( A^{up} \left( u_{(1,k-\frac{1}{2})}^{n}, \tilde{c}_{(1,k-1)}^{n}, \tilde{c}_{(1,k)}^{n} \right) - A^{up} \left( u_{(1,k+\frac{1}{2})}^{n}, \tilde{c}_{(1,k)}^{n}, \tilde{c}_{(1,k+1)}^{n} \right) \right). $$
Now, for $j \in 2, ..., N_r - 1$, we obtain

$$ 
\begin{aligned}
& \tilde{c}_{(j,k)}^{n+1} + \frac{D\Delta t}{\Delta r^2} \left( -\frac{r_{j-\frac{1}{2}}}{r_{j-1}} \tilde{c}_{(j-1,k)}^{n+1} + \frac{r_{j-\frac{1}{2}} + r_{j+\frac{1}{2}}}{r_j} \tilde{c}_{(j,k)}^{n+1} - \frac{r_{j+\frac{1}{2}}}{r_{j+1}} \tilde{c}_{(j+1,k)}^{n+1} \right) + \frac{D\Delta t}{r_j^2 \Delta \theta^2} \left( -\tilde{c}_{(j,k-1)}^{n+1} + 2\tilde{c}_{(j,k)}^{n+1} - \tilde{c}_{(j,k+1)}^{n+1} \right) \\
& = \tilde{c}_{(j,k)}^n + \frac{\Delta t}{\Delta r} \left( A^{up} \left( u_{(j-\frac{1}{2},k)}^n, \tilde{c}_{(j-1,k)}^n, \tilde{c}_{(j,k)}^n \right) - A^{up} \left( u_{(j+\frac{1}{2},k)}^n, \tilde{c}_{(j,k)}^n, \tilde{c}_{(j+1,k)}^n \right) \right) \\
& \quad + \frac{\Delta t}{r_j^2 \Delta \theta} \left( A^{up} \left( u_{(j,k-\frac{1}{2})}^n, \tilde{c}_{(j,k-1)}^n, \tilde{c}_{(j,k)}^n \right) - A^{up} \left( u_{(j,k+\frac{1}{2})}^n, \tilde{c}_{(j,k)}^n, \tilde{c}_{(j,k+1)}^n \right) \right).
\end{aligned}
\quad (33) $$

Finally, for $j = N_r$, we get

$$ 
\begin{aligned}
& -\frac{k_{\text{off}}\Delta t}{\Delta r}\mu_k^{n+1}\left(1+\frac{k_{\text{on}}\Delta t}{\Delta r}+\frac{D\Delta t}{\Delta r^2}\frac{r_{N_r-\frac{1}{2}}}{r_{N_r}}\right)\tilde{c}_{(N_r,k)}^{n+1} - \frac{D\Delta t}{\Delta r^2}\frac{r_{N_r-\frac{1}{2}}}{r_{N_r-1}}\tilde{c}_{(N_r-1,k)}^{n+1} \\
& +\frac{D\Delta t}{r_{N_r}^2\Delta \theta^2}\left(-\tilde{c}_{(N_r,k-1)}^{n+1}+2\tilde{c}_{(N_r,k)}^{n+1}-\tilde{c}_{(N_r,k+1)}^{n+1}\right) = \tilde{c}_{(N_r,k)}^n + \frac{\Delta t}{\Delta r}\left(A^{up}\left(u_{(N_r-\frac{1}{2},k)}^n, \tilde{c}_{(N_r-1,k)}^n, \tilde{c}_{(N_r,k)}^n\right)\right) \\
& \quad +\frac{\Delta t}{r_{N_r}^2\Delta \theta}\left(A^{up}\left(u_{(N_r,k-\frac{1}{2})}^n, \tilde{c}_{(N_r,k-1)}^n, \tilde{c}_{(N_r,k)}^n\right) - A^{up}\left(u_{(N_r,k+\frac{1}{2})}^n, \tilde{c}_{(N_r,k)}^n, \tilde{c}_{(N_r,k+1)}^n\right)\right).
\end{aligned}
\quad (34) $$

The membrane activation equation writes in polar coordinates

$$ \partial_t \mu(t, R, \theta) = k_{\text{on}}c(t, R, \theta) - k_{\text{off}}\mu(t, R, \theta), \quad \text{on } C(0, R). \quad (35) $$

At each time step, the implicit discretization of equation (35) for $k \in \{1, ..., N_\theta\}$ writes

$$ - \Delta t k_{\text{on}} \tilde{c}_k^{n+1} + (1 + \Delta t k_{\text{off}}) \mu_k^{n+1} = \mu_k^n. \quad (36) $$

For simplicity, we treat both the free and activated concentrations in the same linear problem of unknown

$$ \mathcal{E}^n = (\tilde{c}_{(1,1)}^n \dots \tilde{c}_{(1,N_\theta)}^n \tilde{c}_{(2,1)}^n \dots \tilde{c}_{(2,N_\theta)}^n \dots \tilde{c}_{(N_r,N_\theta)}^n \tilde{\mu}_1^n \dots \tilde{\mu}_{N_\theta}^n)^T. $$

We have at each time step

$$ \left( I_{(N_r+1)N_\theta} + \frac{\Delta t}{\Delta r^2} \mathcal{A} \right) \mathcal{E}^{n+1} = \left( I_{(N_r+1)N_\theta} - \frac{\Delta t}{\Delta r} \mathcal{B}^n \right) \mathcal{E}^n, $$

where $\mathcal{A}$ is the diffusion matrix, and $\mathcal{B}^n$ the advection matrix.
The diffusion matrix writes

$$
A = \begin{pmatrix}
\frac{r_{1/2}}{r_1} I_{N_\theta} & -\frac{r_{1+1/2}}{r_2} I_{N_\theta} & & & & \\
& \ddots & & & & \\
& & \frac{r_{j-1/2}+r_{j+1/2}}{r_j} I_{N_\theta} & -\frac{r_{j+1/2}}{r_{j+1}} I_{N_\theta} & & \\
& & & \ddots & & \\
& & & & \ddots & \\
& & & & & \frac{r_{N_r-1/2}}{r_{N_r-1}} I_{N_\theta} \left( \begin{array}{cc} \frac{r_{N_r+1/2}}{r_{N_r}} + \Delta r k_{\text{on}} & I_{N_\theta} \\ -\Delta r^2 k_{\text{on}} & I_{N_\theta} \end{array} \right) & -\Delta r k_{\text{off}} I_{N_\theta} & \\
& & & & \Delta r^2 k_{\text{off}} I_{N_\theta} & \\
& & & & & \left( \begin{array}{cc} \frac{1}{r_1^2} A & \\ & \ddots \\ & & \ddots \\ & & \frac{1}{r_{N_r-1}^2} A & \\ & & & \frac{1}{r_{N_r}^2} A & \\ & & & & 0_{N_\theta} \end{array} \right) & , (37)
\end{pmatrix}
$$

where the second matrix accounts for the angular diffusion, with $A \in M_{N_\theta}(\mathbb{R})$ the classical diffusion matrix with periodic flux boundary conditions:

$$
A = \begin{pmatrix}
2 & -1 & & & -1 \\
-1 & 2 & \ddots & & \\
& \ddots & \ddots & \ddots & \\
& & 2 & -1 & \\
-1 & & -1 & 2 &
\end{pmatrix} . \qquad (38)
$$

Now, for the advection term $A^{up}$ defined by equation (31), we write $(u)^+ = \max(u, 0)$ and $(u)^- = \min(u, 0)$ so that $A^{up}(u, \tilde{c}_{j,k}, \tilde{c}_{j,k+1}) = (u)^+ c_{j,k} + (u)^- c_{j,k+1}$. Therefore, we introduce the following diagonal matrices for $j \in \{1, \dots, N_r\}$, $U_{j+\frac{1}{2}}^+ \in M_{N_\theta}(\mathbb{R})$ and $U_{j+\frac{1}{2}}^- \in M_{N_\theta}(\mathbb{R})$:

$$
U_{j+\frac{1}{2}}^{\pm} = \begin{pmatrix} \cdots \\ (u_{(j+\frac{1}{2},k-1)}^{n})^{\pm} & (u_{(j+\frac{1}{2},k)}^{n})^{\pm} & (u_{(j+\frac{1}{2},k+1)}^{n})^{\pm} \\ & \ddots & \vdots \\ & & \ddots \end{pmatrix} .
$$
Thus we can define the advection matrix:

$$
\mathcal{B}^n =
\begin{pmatrix}
\begin{pmatrix} U_{\frac{3}{2}}^{+} & U_{\frac{3}{2}}^{-} \\ \vdots & \vdots \\ U_{j+\frac{1}{2}}^{+} & U_{j+\frac{1}{2}}^{-} \\ \vdots & \vdots \\ U_{N_r-\frac{1}{2}}^{-} & 0 \\ 0 & 0_{N_\theta} \end{pmatrix} & -
\begin{pmatrix}
0 \\ U_{\frac{3}{2}}^{+} & \cdots & U_{j-\frac{1}{2}}^{+} & U_{j-\frac{1}{2}}^{-} \\ \vdots & \ddots & \vdots & \vdots \\ U_{N_r-\frac{1}{2}}^{+} & U_{N_r-\frac{1}{2}}^{-} & 0 & 0_{N_\theta} \end{pmatrix} \\
+ &
\begin{pmatrix}
\frac{1}{r_1^2} B^n \\ \vdots \\ \frac{1}{r_{N_r-1}^2} B^n \\ \vdots \\ \frac{1}{r_{N_r}^2} B^n \\ 0 \end{pmatrix}
\end{pmatrix}
\quad (39)
$$

where the discrete advection matrix $B^n \in M_{N_\theta}(\mathbb{R})$ with periodic flux condition on the boundary is defined as in Allaire (2005)

$$
\mathcal{B}^n =
\begin{pmatrix}
(u_{\frac{3}{2}}^n)^+ & (u_{\frac{3}{2}}^n)^- \\
\vdots & \vdots \\
(u_{j+\frac{1}{2}}^n)^+ & (u_{j+\frac{1}{2}}^n)^- \\
\vdots & \vdots \\
(u_{N_\theta+\frac{1}{2}}^n)^- & (u_{N_\theta+\frac{1}{2}}^n)^+
\end{pmatrix}
\tag{40}
\\[\bigskipamount]
-
\begin{pmatrix}
(u_{\frac{1}{2}}^n)^- \\
(u_{\frac{3}{2}}^n)^+ \\
\vdots \\
(u_{j-\frac{1}{2}}^n)^+ & (u_{j-\frac{1}{2}}^n)^- \\
\vdots & \vdots \\
(u_{N_\theta-\frac{1}{2}}^n)^+ & (u_{N_\theta-\frac{1}{2}}^n)^-
\end{pmatrix}.
$$

3 Results

The discretization scheme was implemented using MATLAB. We performed some numerical
simulations to test the scheme’s numerical convergence. Since the problem has a boundary
nonlinearity, comparing a numerical solution to an exact one is out of reach in general.

We consider the case of a polarised initial condition that writes $\tilde{c}(0, r, \theta) = \cos(\theta - \pi) + 1$, $\tilde{\mu}(0, \theta) = 0.5\tilde{c}(0, R, \theta)$. We fix some parameters: $k_{\text{off}} = D = k_d = 1$; $\delta = \gamma = 2$.

Illustrative example

By changing the value of the activation rate $k_{\text{on}}$, we observe qualitatively different stationary solutions. The figure 1 shows the time evolution in the molecular concentration in the cell body
in two characteristic cases. The left figure shows a typically non motile profile, while the right figure displays a polarisation situation.

<figure><img src="image_3.png" /><figcaption>Figure 1: Numerical simulation of the evolution in the molecular concentration over time. Discretization parameters: $\Delta t = 10^{-2}$, $\Delta r = \Delta \theta = 2\pi/120 \approx 5.3 * 10^{-2}$. Parameters: $R = 1.5$, $k_{\text{off}} = D = k_d = 1$; $\delta = \gamma = 2$. Initial condition: $\tilde{c}(0, r, \theta) = \cos(\theta - \pi) + 1$, $\tilde{\mu}(0, \theta) = 0.5\tilde{c}(0, R, \theta)$.</figcaption></figure>

Note that for $k_{\text{on}} = 0$, it is clear that the stationary concentration in activated molecules is 0 everywhere on the boundary, so that $v = 0$ and the nonlinearity disappears. In the following, we rather focus on cases where the stationary state can be asymmetric.

## Non-zero stationary polarisation

We also performed the same simulation for $k_{\text{on}} = 0.3$ and $R = 1$, with an angular discretization step $\Delta \theta = 2\pi/160 \approx 3.9 * 10^{-2}$, and varying values for $\Delta r$ and $\Delta t$. The system is considered at numerical steady state when the concentration $\mu$ in activated molecules has stayed unchanged for 1 numerical hour. Then, we obtain the time to attain the steady state, as well as the stationary polarisation of the system, quantified up to a constant by the cell velocity vector.

We consider the following parameter values:

$$ \begin{array}{c|ccccc} \Delta r & 5 * 10^{-3} & 10^{-2} & 2 * 10^{-2} & 2.5 * 10^{-2} & \\ \Delta t & 10^{-3} & 2 * 10^{-3} & 5 * 10^{-3} & 10^{-2} & 1.5 * 10^{-2} \end{array} $$

The figure 2 shows that the time to attain the stationary state is a consequence of the polarisation state of the stationary solution. Indeed, considering that the initial condition is polarised, and the simulation shows the system's depolarisation towards a less polarised steady state, then the more polarised it is, the sooner it is attained.

We also clearly notice that the smaller $\Delta t$ gets, the more polarised the stationary state gets. This trend can be distinguished independently of $\Delta r$. The radial step has a smaller effect, but more precise grids are correlated with lower polarised solutions.

Overall, the simulations show that a fine time step is fundamental to catch the relevant stationary type of behaviour. The polarisation module being decreasing with the time step, we can infer that the numerical solution approach a polarised state. However, this also shows that
<figure><img src="image_4.png" /><figcaption>Figure 2: Left: Time to attain the stationary state and Right: stationary polarisation module for varying $\Delta r$ and $\Delta t$, and $k_{\text{on}} = 0.3$.</figcaption></figure>

the previous numerical simulations are better understood in a qualitative sense rather than for
quantitative purposes.

Low stationary polarisation

Finally, we perform the same numerical test with $k_{\text{on}} = 0.1$ to check how the discretization steps
may generate an error between a polarised and a non polarised stationary solution. We obtain
the values depicted in figure 3.

We can see that the stationary polarisation levels are very low compared to the case where $k_{\text{on}} = 0.3$. The same trend appears for the effect of $\Delta t$, while the effect of $\Delta r$ is less visible since the nonlinear term is small.

To determine if these polarisation levels could be comparable to a true symmetric stationary
state, we performed the same simulation for $k_{\text{on}} = 0$ and for the most precise time-space grid.
The stationary polarisation module was approximately equal to $1.63 * 10^{-4}$, so that no clear
distinction can be made between these cases.

4 Conclusion

In this work, we have presented the finite volume discretization of a multiscale model for 2D cell
crawling migration consisting in a Darcy fluid dynamics coupled to a Poisson problem and to a
nonlinear and non-local reaction-advection-diffusion problem for the concentration in a molecular
specie. The simulation of numerical solutions of this type of problems is very useful since the
mathematical analysis is necessarily limited, whereas the model show varied behaviours.

The discretization method showed good qualitative numerical result. In particular, the
molecular mass is preserved numerically. However, in critical cases, the scheme seems not able
to make the distinction between polarised and unpolarised states. A natural continuation will
consist in taking into account the interaction with the environment (chemical signal, mechanical
obstacle) as a bias for motion. Finally, further studies should be based on an implicit treatment
of the nonlinearity Cancès et al. (2017b,a).
<figure><img src="image_5.png" /><figcaption>Figure 3: Left: Time to attain the stationary state and Right: stationary polarisation module for varying $\Delta r$ and $\Delta t$ and $k_{\text{on}} = 0.1$.</figcaption></figure>

References

Allaire, G. (2005). *Analyse numérique et optimisation: Une introduction à la modélisation mathématique et à la simulation numérique*. Editions Ecole Polytechnique.

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

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

Cancès, C., Chainais-Hillairet, C., and Krell, S. (2017a). A nonlinear discrete duality finite volume scheme for convection-diffusion equations. In *International Conference on Finite Volumes for Complex Applications*, pages 439-447. Springer.

Cancès, C., Chainais-Hillairet, C., and Krell, S. (2017b). Numerical analysis of a nonlinear free-energy diminishing discrete duality finite volume scheme for convection diffusion equations. *arXiv preprint arXiv:1705.10558*.

Etchegaray, C. (2016). *Mathematical and numerical modelling of cell migration*. Theses, Université Paris-Saclay.

Etchegaray, C., Meunier, N., and Voituriez, R. (2017a). A 2D deterministic model for cell crawling, a minimal multiscale approach. Submitted.

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

Joanny, J.-F. and Prost, J. (2009). Active gels as a description of the actin-myosin cytoskeleton. *HFSP J*, 3(2):94-104.

Kruse, K., Joanny, J. F., Jülicher, F., Prost, J., and Sekimoto, K. (2005). Generic theory of active polar gels: a paradigm for cytoskeletal dynamics. *Eur Phys J E Soft Matter*, 16(1):5-16.
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., and Voituriez, R. (2015). Actin flows mediate a universal coupling between cell speed and cell persistence. *Cell*, 161(2):374-86.