# Cell motility in confinement: a computational model for the shape of the cell

Florence Hubert, Meriem Jedouaa, Imene Khames, Julien Olivier, Olivier

Theodoly, Ariane Trescases

► To cite this version:

Florence Hubert, Meriem Jedouaa, Imene Khames, Julien Olivier, Olivier Theodoly, et al. Cell motility in confinement: a computational model for the shape of the cell. ESAIM: Proceedings and Surveys, 2016, 55, pp.148-166. 10.1051/proc/201655166. hal-03352663

HAL Id: hal-03352663

https://hal.science/hal-03352663

Submitted on 24 Sep 2021

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
CELL MOTILITY IN CONFINEMENT: A COMPUTATIONAL MODEL FOR THE SHAPE OF THE CELL

FLORENCE HUBERT<fnref n="1" />, MERIEM JEDOUAA<fnref n="2" />, IMENE KHAMES<fnref n="3" />, JULIEN OLIVIER<fnref n="4" />, OLIVIER THEODOLY<fnref n="5" /> AND ARIANE TRESCASES<fnref n="6" />

**Abstract.** While cells typically tend to spread their cytoplasm in a flat and thin lamellipodium when moving on a flat substrate, it is widely observed that the cytoplasm has a compact shape in microchannels, tending to fulfill the cross-section of the microchannel. We propose a minimal mathematical model for a 2D test case which describes the cell lamellipodium deformations when confined in a channel. We then go through a numerical investigation of this mathematical model and show that it allows to recover qualitatively the physiological characteristics of the confined cell.

**Résumé.** Alors que les cellules ont généralement tendance à présenter un cytoplasme étendu en un large lamellipode extrêmement fin lors du déplacement sur un substrat plat, il est communément observé que le cytoplasme prend une forme compacte lors du déplacement dans des micro-canaux, remplissant au possible le volume contenu dans le micro-channel. Nous proposons un modèle mathématique minimal pour un cas test en 2D qui décrit les déformations du lamellipode en confinement. Nous proposons une exploration numérique de ce modèle mathématique et nous montrons qu'il permet de retrouver qualitativement les caractéristiques physiologiques de la cellule confinée.

INTRODUCTION

Cell migration inside extracellular matrix network is known to play a crucial role in a variety of physiological and pathological processes, as, for example, wound healing, immune surveillance and inflammation [11], or cancer growth [22]. Moreover, the research on cell migration already has come to some biomedical applications, such as the regeneration of tissues [4, 25].

As a toy model to understand cell motility in the extracellular matrix network, many studies have considered the movement of cells in simplified confined environments, and particularly, in micro-channels. Both experimental and theoretical works have been done, questioning in particular the role of adhesion mechanisms for cell displacement in confined settings, while it is known that adhesion plays a crucial role in cell motility on a flat substrate [10–12, 16]. Though it is not intuitive, it was also shown (experimentally and theoretically, [13]) that confined cells could possibly reach velocities significantly larger than the velocities observed on flat substrates.

<footnote n="1">¹ Institut de Mathématiques de Marseille, I2M, UMR 7373, 39 rue F. Joliot Curie, 13453 Marseille, France</footnote>

<footnote n="2">² Laboratoire Jean Kuntzmann, Univ. Grenoble Alpes and CNRS, Grenoble, France</footnote>

<footnote n="3">³ Laboratoire de Mathématiques, INSA de Rouen, Avenue de l'université, F-76801 Saint-Etienne du Rouvray, France</footnote>

<footnote n="4">⁴ Institut de Mathématiques de Marseille, I2M, UMR 7373, 39 rue F. Joliot Curie, 13453 Marseille, France</footnote>

<footnote n="5">⁵ Laboratoire Adhésion & Inflammation, INSERM UMR S 1067, CNRS UMR 7333, Aix-Marseille Université, campus de Luminy, 163 av de Luminy, 13009 Marseille, France</footnote>

<footnote n="6">⁶ DPMMS, Univ. Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK</footnote>
On the other hand, a tremendous effort has been made to understand the motility on flat surfaces. Studies
typically focus on different aspect of cell motility, such as intern mechanisms, interaction with the environment,
communication between cells, polarization... Among them, many aim at describing and explaining the specific
shape of the cell (see for example [14,17,23] and references therein). Indeed, when migrating on a flat substrate,
the typical cell tends to extend its cytoplasm into an extremely thin layer which covers a wide area of substrate.
This flat domain of cytoplasm is called lamellipodium for its remarkable thinness, measured in [1] as 110–160
nm (lamella, lat.: thin layer of material; podos, gr.: foot).

By contrast, experiments show that the shape of the leading edge protrusion in micro-channels is not flat.
The lamellipodium contracts and tends to fulfill a compact volume of the micro-channel. For more details, see
for example [24,26] and references therein for a comparison between migration on flat substrates and in confined
environment. The duality of the shape of the lamellipodium is summarized in figure 1.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Duality of the lamellipodium on a flat substrate (left) and in a micro-channel (right)</figcaption></figure>

In this paper, our aim is to better understand some of the essential mechanisms operating in the shaping of
the lamellipodium in confined domain. We will essentially focus on the influence of the adhesion forces. For that,
we present two toy-computationals models taking into account some of the basic mechanisms of cell migration.
We present numerical simulations associated with these models. This work is a first step in understanding the
impact of the confinement.

We now present the main mechanisms of cell motility, as described in the biology literature. The movement
can typically be decomposed in four steps: in the first step, the lamellipodium extends at the leading edge
(protrusion) by polymerization of actin. More precisely, polymerization in the cytoplasm pushes forward the
membrane. The second step is the adhesion of the membrane at the leading edge to the substrate of the walls.
Adhesion is crucial to acquire momentum (at least on a flat substrate). At this stage, the membrane has adhered
to the wall by its two extremities and the membrane is under tension (as a consequence of the protrusion at the
leading edge and the adhesion at the trailing edge). The third step is deadhesion at the trailing edge, which
allows to release the tension when the cell displaces: this is the fourth step. To summarize, the three main
biological mechanisms are therefore: polymerization at the membrane, adhesion and deadhesion to the wall,
and membrane tension. They are schematically represented (in the case of a micro-channel) on Figure 2.

This paper is organized as follows: In Section 1 are introduced the mathematical toy-models (the *full* model and the *simplified* model) we propose for the cell in a micro-channel with a careful description of the different mechanisms we want to take into account for the motility of the lamellipodium. In Section 2, the numerical treatment and implementation of the full model is presented. In Section 3 is proposed a bench of numerical simulations corresponding to two different configurations of the lamellipodia when entering in the micro channel. In the last section, we sum up the conclusions we have derived from our study and lay down the various questions we wish to work on in further, more detailed work.
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Mechanisms of the cell motility</figcaption></figure>

# 1. MODELLING THE LAMELLIPODIA MOTILITY

## 1.1. Mechanical description of the lamellipodium

The focus of this section is to introduce mathematical models to describe the lamellipodia deformation and movement in a micro-channel. This model inspired by [23] takes into account the lamellipodium immersed in the fluid fulfilling the micro-channel, the polymerization of the actin network, the adhesion of the lamellipodium to the walls and the membrane forces which are here limited to the surface tension. Moreover, it describes the actin dynamics inside the lamellipodia.

We consider the two-dimensional case which corresponds to a toy model of the lamellipodium into the channel and adhering on a two-dimensional substrate.

The basis of the model is to consider the lamellipodium as an “active fluid”: the complex material consisting of a mix of intracellular liquid and actin network will be described by a continuous medium. To describe the properties of this medium, we have to note that the scale of the setting allows us to use low Reynolds number approximation. Now following the discussion in [21] we argue as follows: the intracellular liquid is an incompressible Newtonian fluid and the actin network is an elastic system. We could then go and model our lamellipodium as a viscoelastic material, inhering properties of both its constituents. But we argue that the scale of the outer mechanical forces that would act on the lamellipodium are too large to see the elastic response, driving us to describe the lamellipodium as an incompressible Newtonian fluid with an apparent viscosity $\mu_1$. We are aware that this description is rough and that others have argued that one may take into account the elastic properties of the material by dropping the incompressibility condition, because the medium is slightly compressible (due to the elastic network within) and this is made possible by neglecting the force of pressure [21]. Here for the sake of simplicity we choose not to follow this line of reasoning. Thus our medium will be described by a Stokes law for incompressible fluids with apparent viscosity $\mu_1$. Then the medium would be described by its velocity field $U$ and the pressure field inside it $P$.

## 1.2. Forces driving the lamellipodium

We describe now the forces acting on this medium. In our model we take into account three possible effects on the lamellipodium motion.

First, there is the interaction with the surrounding fluid in the channel. We consider the surrounding fluid to be also a Newtonian fluid at low Reynolds, with viscosity $\mu_f$. The interaction between the two media is taken into account *via* the Level-Set method (described in Section 2).

Then we consider a force acting on the membrane of the lamellipodium called surface tension. This force is a normal force on the membrane proportional to the local curvature of the surface.
Finally, the force that interests us the most is the force coming from the polymerization of actin. Let us describe the mechanism we take. According to experimental observations, the lamellipodium deformation and its movement inside the channel are mainly due to the mechanical properties of the actin network. Hence, actin polymerization/depolymerization is the main mechanism directing the movement of the cell. Moreover, experiments show that the actin assembly is maximal near the membrane and that the actin network polymerizes perpendicularly to the cell/wall interface.

In fact, the polymerization of the actin filaments induces a pressure at the membrane which pushes the plasma membrane outward. The intensity of this pushing force depends on the concentrations of actin and actin binding proteins. The higher the actin concentration, the more the pushing force applied on the membrane.

Note that we may have to distinguish actin under its “free” form whose density we note $\rho_a^{cyt}$ and actin under its polymerized form whose density we note $\rho_a$.

With these observations, we propose a polymerization force which is confined to the cell rim, normal to the membrane and proportional to the polymerized actin density $\rho_a$.

## 1.3. Actin dynamics

What is left to describe is how actin moves around inside the lamellipodium. Note that for the sake of simplicity, we put in our model a fixed amount of total actin to separate what happens inside our lamellipodium piece from the rest. This constant total amount of actin is noted $\rho_a^{tot}$.

We consider the “free” actin to be small molecules that have the ability to diffuse extremely fast inside the lamellipodium. Consequently we will consider that the density of “free” actin is homogeneous in space.

We describe now the dynamics of the polymerized actin. We choose to do it *via* a reaction-transport-diffusion process. The reaction term models the process of depolymerization of the actin. The diffusion is standard at this level of description. Finally, we say that the network itself is transported with the velocity of the medium $U$. Once again we make the assumption that the velocity $U$ describes both the cellular liquid *and* the actin network for the sake of simplicity. Of course, a much refined description could be done by unravelling the effect of the actin network with respect to the intracellular liquid.

When we use the previous dynamics on the actin density we say that we use the *full model*. In our first batch of simulation however we wanted to test the purely mechanical effects on our system. To do so we use the *simplified model*, where we arbitrarily set the polymerized actin density to a fixed value near the plasmic membrane, forgetting every detail of actin dynamics. This assumption, in agreement with the results of the experiments performed in [23], decouples the effect of the protrusion at the membrane from other mechanisms. In the end, this will allow to calibrate the intensity of the protrusion in the complete model.

## 2. MATHEMATICAL FRAMEWORK

Our main interest being the deformation of the lamellipodium, one major goal of this paper is to translate the proposed biophysical ingredients of Section 1 into mathematical equations that are able to efficiently treat the problem of the movement of the membrane.

In order to describe the interaction between the external fluid and the lamellipodium, a purely eulerian formulation is used [8, 9]. To achieve this, the interface is captured thanks to a level set approach.

We first present an overview of the level set method. Then, we introduce the proposed model with a careful description of the polymerization of the actin filaments, the adhesion of the lamellipodium to the walls and the membrane surface tension. Finally, we present the hypothesis made to obtain the simplified model we use to perform numerical simulations.

### 2.1. Outline of the level set method

Introduced in [19] to treat problems involving interfaces, the level set approach is an efficient tool to capture the interface which presents several advantages: it is easy to implement, and topological changes are directly
handled. The general idea of the level set method is to define a scalar function that takes the value zero at the location of the interface to capture.

Let $\Omega$ be a bounded domain in $\mathbb{R}^d$ ($d = 2$ or $d = 3$) partitioned into two subdomains $\Omega_1(t)$ and $\Omega_2(t)$, and let $\Gamma(t)$ be the interface between $\Omega_1(t)$ and $\Omega_2(t)$. The aim is to follow the evolution of the interface $\Gamma(t)$ that is defined as the zero value of a level set function $\phi(\cdot, t)$:

$$ \Gamma(t) = \{x \in \Omega, \phi(x, t) = 0\}. $$

The level set function $\phi(\cdot, t)$ has to be Lipschitz continuous in $x$ on the whole domain $\Omega$ (for all time $t$). The displacement of the interface is deduced from the evolution of the level set function $\phi$. More precisely, consider that the interface is transported with the velocity $V$ of the flow in the whole domain $\Omega$. Then, the level set function is the solution of the scalar transport equation:

$$ \partial_t \phi(t, x) + V(t, x) \cdot \nabla_x \phi(t, x) = 0. \quad (1) $$

We usually define the level set function as a signed distance function that is regular in each corresponding domain,

$$ \phi(x) = \begin{cases} -d(x, \Gamma(t)) & x \in \Omega_1(t), \\ d(x, \Gamma(t)) & x \in \Omega_2(t), \end{cases} \quad (2) $$

where

$$ d(x, \Gamma(t)) = \min_{y \in \Gamma(t)} \|x - y\|. $$

Using this level set function we can define a regularized “Heaviside function $H_\varepsilon$” that is to say a function with value $\approx 1$ in $\Omega_2(t)$ and $\approx 0$ in $\Omega_1(t)$ and transition from 0 to 1 in a zone of size $\approx \varepsilon$ around $\Gamma(t)$: first define

$$ h_{\varepsilon}(s) = \begin{cases} 0 & s \le -\varepsilon, \\ \frac{1}{2}\left(1 + \frac{s}{\varepsilon} + \frac{\sin(\frac{\pi s}{\varepsilon})}{\pi}\right) & |s| \le \varepsilon, \\ 1 & s \ge \varepsilon, \end{cases} \quad (3) $$

and then let $H_\varepsilon(t, x) = h_\varepsilon(\phi(t, x))$. We will also use a regularized “Dirac” function by setting

$$ z_{\varepsilon}(s) = \begin{cases} 0 & s \le -\varepsilon, \\ \frac{1}{2\varepsilon}\left(1 + \cos\left(\frac{\pi s}{\varepsilon}\right)\right) & |s| \le \varepsilon, \\ 0 & s \ge \varepsilon, \end{cases} \quad (4) $$

so that $\zeta_\varepsilon(t, x) = z_\varepsilon(\phi(t, x))$ is a regularized “Dirac” function with support at the interface. In this case, $\varepsilon$ represents half of the interface thickness.

Moreover geometrical characteristics of the curve $\Gamma(t)$ such as normal vectors $n$ and curvature $\kappa$ are obtained explicitly using the level set function:

$$ n = \frac{\nabla \phi}{|\nabla \phi|}, \quad \kappa = \nabla \cdot n. $$

## 2.2. Equations of the models

The spatial domain of interest is noted $\Omega$, and is taken rectangular to model a portion of the micro-canal. The upper and lower boundaries of the rectangular $\Omega$, called $\Gamma_d$ and resp. $\Gamma_b$, are the physical boundaries of the micro-canal, while the left and right boundaries of $\Omega$, called $\Gamma_a$ and resp. $\Gamma_c$, are artificial/numerical boundaries (see Figure 3). We study the evolution up to a time $T > 0$.

As previously said, both the lamellipodium and the external fluids are Stokes fluid with respective viscosities $\mu_1$ and $\mu_f$. They do not mix and they occupy the full domain $\Omega$. To treat their interaction in a level-set method,
we define the velocity field $U(x,t)$ and the pressure field $P(x,t)$ in the whole $\Omega_T = \Omega \times [0, T]$ to be the velocity and pressure of the lamellipodium fluid when $x$ is in the space occupied by the lamellipodium and those of the external fluid elsewhere. This allows to write a single equation in the whole domain rather than equations on separate moving domains. Defining $\mu(\phi(x,t)) = \mu_1 + (\mu_f - \mu_1)h_{\varepsilon}(\phi(x,t))$ one can then write

$$ \begin{cases} -\nabla \cdot (2\mu(\phi)D(U)) + \nabla P = F_{\text{ext}} & \text{in } \Omega_T, \\ \nabla \cdot U = 0 & \text{in } \Omega_T, \end{cases} $$

where $D(U) = \frac{\nabla U + (\nabla U)^T}{2}$ is the symmetrized strain rate tensor and $F_{\text{ext}}$ is the sum of the external forces acting on the lamellipodium and the external fluid. Since the lamellipodium/external fluid interaction is already taken into account, $F_{\text{ext}}$ is the sum of two contributions that we describe now.

### Polymerization of the actin network.

The first contribution is the effect of actin polymerization. Following [23] we describe this effect by a stress tensor which is localized on the membrane and pushing towards the exterior of the lamellipodium with an intensity proportional to the polymerization activity, measured in our model through the density of polymerized actin. Generally speaking, the stress tensor would be given by

$$ \sigma_{\text{poly}} = -\eta_a^0 \rho_a \delta_{\partial\Omega_{1,t}} n \otimes n \quad (5) $$

where $n$ is the outer normal at the membrane and $\delta$ is the Dirac function. In the context of level-sets methods, it can be extended to the whole domain via a regularization of the Dirac function through

$$ \sigma_{\text{poly}} = -\eta_a^0 \rho_a z_{\varepsilon}(\phi)n(\phi) \otimes n(\phi), \quad \text{in } \Omega_T. \quad (6) $$

Here $\eta_a^0$ is the actin protrusion coefficient and $n(\phi)$ is the normal to the membrane pointing outward. The actin density $\rho_a$ is also defined on the whole $\Omega_T$ even if their is no meaning to nonzero values of $\rho_a$ outside of the lamellipodium. Note that the force is the divergence of the tensor $\nabla \cdot \sigma_{\text{poly}}$.

### The membrane force

The second contribution to the external forces is the membrane force. We denote as $F_{\text{mem}}$ the surface tension at the membrane, which can be expressed directly by using the level set function:

$$ F_{\text{mem}} = \gamma\kappa(\phi)n(\phi), \quad \text{in } \Omega_T, \quad (7) $$

where $\gamma$ is the tension coefficient.

### The displacement of the membrane

The displacement of the membrane is specified by the transport of the level set function with the flow velocity $U$.

$$ \partial_t \phi + U \cdot \nabla \phi = 0, \quad \text{in } \Omega_T. \quad (8) $$

### Reaction-transport-diffusion equation for the actin network

In our “simplified” model, $\rho_a$ is prescribed so we do not need more equations. It is given by

$$ \rho_a = C \zeta_{\varepsilon_2}(\phi), $$

where $C > 0$. In this manner, $\rho_a$ is positive near the interface at a distance of $\approx \varepsilon_2$. There is no particular reason why $\varepsilon$ and $\varepsilon_2$ should be the same and in numerical simulations we take we take $\varepsilon_2 = 3 \times \varepsilon$.
On the other hand, for the “full” model, we have to write an equation for the reaction-transport-diffusion process. We use the model of [23]:

$$ \partial_t \rho_a = -\nabla \cdot (\rho_a U) + \nabla \cdot (D_a \nabla \rho_a) + f(\rho_a, \rho_a^{cyt}), \quad \text{in } \Omega_{1,T} = \{(x, t), t \in [0, T], x \in \Omega_1(t)\}. \quad (9) $$

Here $D_a$ is an effective diffusion coefficient of the actin network. This constant effective diffusion arises from random events related to various chemical processes. The function $f(\rho_a, \rho_a^{cyt})$ describes the reaction kinetics of the actin filaments : polymerization and depolymerization. The specific form of $f$ we use is again taken from [23] and is classical in the modelling of positive feedback reactions see [18]:

$$ f(\rho_a, \rho_a^{cyt}) = k_b \left( \frac{\rho_a^2}{K_a^2 + \rho_a^2} + k_a \right) \rho_a^{cyt} - k_c \rho_a. $$

The second term is a simple depolymerization term with constant rate $k_c$. The first term models polymerization, which is taken proportional to the cytosolic actin monomer concentration $\rho_a^{cyt}$, and is increasing in $\rho_a$ with saturation. Here, $k_a$ is the base polymerization rate, $k_b$ is the actin polymerization rate, $K_a$ is a positive feedback threshold.

Note that given the hypotheses we made on the dynamics of “free” actin monomers, we can also write

$$ \rho_a^{cyt}(t) = \frac{\left( \rho_a^{tot} - \int_{\Omega_1} \rho_a \right)}{A(t)}. $$

Recall that $\rho_a^{tot}$ is the total amount of actin. The term $A(t)$ is the total area of the lamellipodium.

### Boundary conditions

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Computational domain $\Omega$</figcaption></figure>

To complete our system we have to provide boundary conditions on each of the four sides of the rectangular $\Omega$ and initial conditions.

On the inner and outer boundary $\Gamma_a$ and $\Gamma_c$ (see Figure 3 for the definition of these) we disregard the phenomena happening outside the domain, in the sense that, in this prospective toy model, we pretend that the behaviour of the quantity $U$ is locally (close to the boundary) homogeneous in the horizontal variable, and equal to its value at the (inner) boundary. Therefore, no new force is input in the model at the boundary. For this we choose the homogeneous Neumann boundary conditions. Of course more sophisticated boundary conditions could be considered, such as (4) in [19]. Here our choice is driven by simplicity of formulation and implementation. Therefore, we have for the flow velocity:

$$ \nabla U \cdot n = 0 \quad \text{on } \Gamma_a \cup \Gamma_c. $$

The lamellipodium is in contact with the bottom wall $\Gamma_b$ and the top wall $\Gamma_d$ and adheres to these walls. We assume non permeability at the walls. For the “simplified” model, we assume the adhesion to be very strong, so
that we chose to model both lamellipodium/wall and external fluid/wall interaction by perfect adhesion. This thus reads:

$$U = 0 \text{ on } \Gamma_d \cup \Gamma_b.$$

For the “full” model, we model the adhesion by a friction force. What we call adhesion here is therefore a force which is specifically tangent to the wall of contact. We could argue that adhesion may also involve normal forces. Here, for the sake of simplicity, we consider the polymerization force as the main active mechanism in the direction normal to the boundary. If we denote as $U = (U_x, U_y)^t$ the velocity and as $\sigma(U, P) = 2\mu D(U) - P\text{Id}$ we obtain the following boundary conditions:

$$\begin{cases} U \cdot n = 0 & \text{on } \Gamma_d \cup \Gamma_b, \\ [(\sigma(U, P) + \sigma_{\text{poly}}(\phi))n]_{\text{tan}} = \sigma_{xy} = \mu(\phi)(\partial_x U_y + \partial_y U_x) = -\lambda(\phi)U_{\text{tan}} = -\lambda(\phi)U_x & \text{on } \Gamma_d \cup \Gamma_b. \end{cases}$$

Denoting by $\lambda_f$ and $\lambda_1$ the friction coefficients inside the external fluid and the lamellipodium, the friction field $\lambda$ in the whole domain $\Omega$ is given by:

$$\lambda(\phi) = \lambda_1 + (\lambda_f - \lambda_1)H_{\varepsilon}(\phi). \quad (10)$$

For the density of actin, for the same reasons as for the flow velocity $U$, we also take the homogeneous Neumann boundary conditions:

$$\nabla \rho_a \cdot n = 0, \quad \text{on } \Gamma_a \cup \Gamma_b \cup \Gamma_c \cup \Gamma_d.$$

Again, note that this assumption is not meant to be realistic, since there may be significant exchanges at the (fictive) boundary $\Gamma_a$. One should rather see this assumption as an artificial way to theoretically isolate the mechanisms happening at the edge of the lamellipodium from the rest. Alternatively said, all mechanisms not happening in the fictive domain are neglected.

## The two models

Let us sum up the equations considered. Firstly the *Complete model* writes : find $(U, P, \rho_a, \phi)$ satisfying

$$\left\{
\begin{array}{ll}
-\nabla \cdot (2\mu(\phi)D(U)) + \nabla P = F_{\text{mem}} + \nabla \cdot \sigma_{\text{poly}}, & \text{in } \Omega_T, \\
\nabla \cdot U = 0, & \text{in } \Omega_T, \\
U \cdot n = 0, \quad \mu(\phi)(\partial_x U_y + \partial_y U_x) = -\lambda(\phi)U_x, & \text{on } \Gamma_d \cup \Gamma_b \times [0, T], \\
\nabla U \cdot n = 0, & \text{on } \Gamma_a \cup \Gamma_b \times [0, T], \\
\\
\partial_t \rho_a = -\nabla \cdot (\rho_a U) + D_a \Delta \rho_a + f(\rho_a, \rho_a^{\text{cyt}}), & \text{in } \Omega_{1,T}, \\
\nabla \rho_a \cdot n = 0, & \text{on } \Gamma_a \cup \Gamma_b \cup \Gamma_c \cup \Gamma_d \times [0, T], \\
\\
\rho_a^{\text{cyt}}(t) = \frac{(\rho_a^{\text{tot}} - \int_{\Omega_1} \rho_a)}{A(t)}, & \text{in } [0, T], \\
\\
\partial_t \phi + U \cdot \nabla \phi = 0, & \text{in } \Omega_T, \\
(U \cdot n)_+ \phi & \text{prescribed on } \Gamma_a \cup \Gamma_c.
\end{array}
\right.
\quad (11)$$
Secondly, the Simplified model writes: find $(U, P, \phi)$ solution of the system

$$ 
\left\{
\begin{array}{ll}
-\nabla \cdot (2\mu(\phi)D(U)) + \nabla P = \nabla \cdot \sigma_{\text{poly}} + F_{\text{mem}}(\phi) & \text{in } \Omega_T, \\
\nabla \cdot U = 0 & \text{in } \Omega_T, \\
U = 0 & \text{on } \Gamma_b \cup \Gamma_d \times [0, T], \\
\nabla U \cdot n = 0 & \text{on } \Gamma_a \cup \Gamma_d \times [0, T], \\
\rho_a = C\zeta_{\varepsilon_2}(\phi) & \text{in } \Omega_T, \\
\partial_t \phi + U \cdot \nabla \phi = 0 & \text{in } \Omega_T, \\
(U \cdot n)_+ \phi & \text{prescribed on } \Gamma_a \cup \Gamma_c.
\end{array}
\right.
\tag{12}
$$

Both these systems are complemented with initial conditions.

### 3. NUMERICAL IMPLEMENTATION OF SYSTEM (12)

This section is devoted to the numerical implementation of the system (12).

The system is discretized by a finite difference method on a staggered grid where the pressure and the level set function are located at the center of the mesh cells and the velocity at the center of the sides (see figure 4). The divergence free is computed at the pressure point which enforces the volume constraint accurately. Following the ideas introduced in [5,6], we solve the steady Stokes equations by using a variant of the projection method of Chorin [7] and a fixed point method.

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Staggered grid</figcaption></figure>

#### 3.1. General algorithm

Given a time step $\Delta t$, we set $t^n = n\Delta t$, $U^n \approx U(\cdot, t^n)$, and $\phi^n \approx \phi(\cdot, t^n)$.

At each time step, the general algorithm performs the following steps:

(1) Using a projection method of Chorin type and fixed point method we resolve the Stationnary Stokes equation on $U$.

(2) Resolution of the transport equation on the level set field using a WENO 5 scheme and an Euler explicit scheme in time:

$$ \phi^{n+1} = \phi^n - \Delta t \times U \cdot \nabla \phi^n $$

(3) Perform a fast marching method by resolving the following eikonal equation in the entire computational domain $\Omega$:

$$ | \nabla \phi^n | = 1, \quad \text{in } \Omega $$
using a first order numerical scheme [20] with an horizontal and a vertical space steps $\Delta x$ and $\Delta y$ in two-dimensions we solve:

$$ \max(\max(D_x^{-}\phi_{ij}^n, 0)^2, \min(D_x^{+}\phi_{ij}^n, 0)^2) + \max(\max(D_y^{-}\phi_{ij}^n, 0)^2, \min(D_y^{+}\phi_{ij}^n, 0)^2) = 1 \quad (13) $$

$$ \text{where } D_x^{-}\phi_{ij}^n = \frac{\phi_{ij}^n - \phi_{i-1,j}^n}{\Delta x}, \quad D_x^{+}\phi_{ij}^n = \frac{\phi_{i+1,j}^n - \phi_{i,j}^n}{\Delta x}, \quad D_y^{-}\phi_{ij}^n = \frac{\phi_{ij}^n - \phi_{i,j-1}^n}{\Delta y}, \quad D_y^{+}\phi_{ij}^n = \frac{\phi_{i,j+1}^n - \phi_{i,j}^n}{\Delta y}. $$

## 3.2. Resolution of the stationary Stokes equation

We use an adapted projection method of Chorin type and a fixed point method to solve the Steady Stokes equations. This method is taken from [5, 6]. The idea of the projection method is to solve the conservation of momentum without the pressure term for an intermediate state $U^*$. Then, this solution is projected on the space free-divergence fields $\mathbb{P}$:

$$ U = \mathbb{P}(U^*) = U^* - \nabla\xi. \quad (14) $$

where the projector $\xi$ is solution of the problem:

$$ \left\{ \begin{array}{ll} \Delta\xi = \operatorname{div}(U^*) & \text{in } \Omega, \\ \partial_n\xi = 0 & \text{on } \Gamma_a \cup \Gamma_b, \cup \Gamma_c \cup \Gamma_d. \end{array} \right. $$

To treat the non homogeneous viscosity term we decompose it as follows:

$$ \nabla \cdot (2\mu D(U)) = \mu \Delta U + 2D(U)\nabla\mu + \mu \nabla \operatorname{div} U. $$

Combining with the incompressibility condition $\operatorname{div} U = 0$ and $U = U^* - \nabla\xi$, it comes:

$$ \nabla \cdot (2\mu D(U)) = \mu \Delta U^* - \mu \nabla \Delta \xi + 2D(U)\nabla\mu. $$

Then, setting $p = -\mu(\phi)\Delta\xi$, the intermediate state $U^*$ is solution of (15) where $U^* = (u^*, v^*)^t$.

$$ \left\{ \begin{array}{ll} -\mu\Delta U^* = -\operatorname{div}(U^*)\nabla\mu - 2D(u)\nabla\mu + F_{\text{mem}}(\phi) + \nabla \cdot \sigma_{\text{poly}} & \text{in } \Omega, \\ \partial_x u^* = \partial_{xx}\xi & \text{on } \Gamma_a \cup \Gamma_c, \\ \partial_x v^* = 0 & \text{on } \Gamma_a \cup \Gamma_c, \\ u^* = \partial_x \xi & \text{on } \Gamma_b \cup \Gamma_d, \\ v^* = 0 & \text{on } \Gamma_a \cup \Gamma_c. \end{array} \right. \quad (15) $$

The term depending on $U^*$ in the right hand side of (15) is treated explicitly. In fact, 15 is transformed into the fixed point problem:

$$ \left\{ \begin{array}{ll} -\mu\Delta U_{k+1}^* = -(2D(U_k) + \operatorname{div}(U_k^*)Id)\nabla\mu + F_{\text{mem}}(\phi) + \nabla \cdot \sigma_{\text{poly}} & \text{in } \Omega_T, \\ +U_{k+1}^* = \nabla\xi_k & \text{on } \Gamma_b \cup \Gamma_d, \\ u_{k+1}^* = \partial_{xx}\xi_k & \text{on } \Gamma_a \cup \Gamma_c, \\ v_{k+1}^* = 0 & \text{on } \Gamma_a \cup \Gamma_c. \end{array} \right. \quad (16) $$

If one wants to apply friction boundary condition then the boundary conditions of $U_{k+1}^*$ on $\Gamma_b \cup \Gamma_d$ are changed by:

$$ \left\{ \begin{array}{ll} u_{k+1}^* = -\partial_x \xi^k + \frac{\mu}{\lambda} \partial_y u_k & \text{on } \Gamma_b \cup \Gamma_d, \\ v_{k+1}^* = 0 & \text{on } \Gamma_b \cup \Gamma_d. \end{array} \right. $$

The last step of the projection method is the correction of the velocity $U$:

$$ U = U^* - \nabla\xi. $$
For the numerical simulations, the error parameter $\epsilon_m$ of the fixed point algorithm (see algorithm 1) is fixed to $10^{-2} \min(\Delta_x, \Delta_y)$. The proposed algorithm is very fast as all the equations are solved using the fast Poisson solver FISHPACK.

**Algorithm 1 Fixed point procedure**

$u_0^* = u_0$ is known at initialisation $k = 0$

while $((|| D(U_k) - D(U_{k-1}) || \ge \epsilon_m) \text{ OR } (|| U_k ||_{\Gamma_b \cup \Gamma_d} \ge \epsilon_m) \text{ OR } (|| \partial_n U_k ||_{\Gamma_a \cup \Gamma_c} \ge \epsilon_m))) \text{ do}$
$U_k^*, U_k = \mathbb{P}(U_k^*)$ and $\nabla \xi_k = U_k^* - \mathbb{P}(U_k^*)$ are known

Find $U_{k+1}^*$ solution of the following system:

$$ \left\{ \begin{array}{ll} -\mu \Delta U_{k+1}^* = -(2D(U_k) + \text{div}(U_k^*) Id) \nabla \mu + F_{\text{mem}}(\phi) + \nabla \cdot \sigma_{\text{poly}}, & \text{in } \Omega_T, \\ +U_{k+1}^* = \nabla \xi_k, & \text{on } \Gamma_b \cup \Gamma_d, \\ u_{k+1}^* = \partial_{xx} \xi_k, & \text{on } \Gamma_a \cup \Gamma_c, \\ v_{k+1}^* = 0, & \text{on } \Gamma_a \cup \Gamma_c. \end{array} \right. \qquad (17) $$

Resolve the equation on $\xi_{k+1}$:

$$ \left\{ \begin{array}{ll} \Delta \xi_{k+1} = \text{div} U_{k+1}^* & \text{in } \Omega, \\ \partial_n \xi_{k+1} = 0 & \text{on } \Gamma_a \cup \Gamma_b, \cup \Gamma_c \cup \Gamma_d. \end{array} \right. \qquad (18) $$

Correct the velocity:

$$ U_{k+1} = U_{k+1}^* - \nabla \xi_{k+1}. $$

$k = k + 1$

end while

# 4. NUMERICAL ILLUSTRATIONS

In this section, we present the numerical results obtained with the proposed models. We first present the simulations obtained with the simplified model (12) for two different configurations of the lamellipodium, using the algorithm described in Section 3.

Then a Freefem++ implementation of the full model (11) is performed.

## 4.1. Choice of the parameters

Our choice of parameters mainly relies on the parameters used in [23]. In [23], the following quantities measured by experiments are given: $\eta_a^0, k_a, k_b, K_a, k_c, D_a, \gamma$. Then, the total density of actin $\rho_a^{\text{tot}}$ is simply computed as the ratio of the total amount of actin $m_a^{\text{tot}}$ (also taken from [23]) and the average area of the cell $A_{\text{cell}}$, that is, $\rho_a^{\text{tot}} = m_a^{\text{tot}}/A_{\text{cell}} = 800/200 = 4 \mu\text{m}^{-2}$. The viscosity of the surrounding fluid, mainly composed of water, is taken equal to the viscosity of water. To simplify, we first consider (in the simplified model) that the viscosity of the actin flow and the viscosity of the surrounding fluid are the same, while in the full model, we consider that the actin flow is less viscous than the surrounding fluid: $\mu_1 = 0.5\mu_f$.

In the case where our model prescribes the actin density close to the membrane, we take the value from [2] (see Table 1).

Finally, in the simplified model, we model the adhesion of the membrane at the contact with the channel by homogenous Dirichlet boundary conditions. In the full model, the adhesion is modeled by a friction force. Up to our knowledge, there is no data available for the (adhesion) friction coefficient in this context. However, what is known is an effective friction coefficient (eg. [23]), which does not take into account the adhesion process.
<table><caption>TABLE 1. Numerical values of the parameters for the model (12)</caption><thead><tr><th></th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td>μ<sub>1</sub></td><td>effective viscosity of actin flow</td><td>0.5·10<sup>3</sup>–10<sup>3</sup> pN.s/μm</td></tr><tr><td>μ<sub>f</sub></td><td>viscosity of the surrounding fluid</td><td>10<sup>3</sup> pN.s/μm</td></tr><tr><td>η<sup>0</sup><sub>a</sub></td><td>F-actin protrusion coefficient</td><td>560 pN.μm<sup>2</sup></td></tr><tr><td>k<sub>a</sub></td><td>base polymerization rate</td><td>0.01 s<sup>-1</sup></td></tr><tr><td>k<sub>b</sub></td><td>F-actin polymerization rate</td><td>10 s<sup>-1</sup></td></tr><tr><td>K<sub>a</sub></td><td>positive feedback threshold</td><td>1 μm<sup>-2</sup></td></tr><tr><td>k<sub>c</sub></td><td>F-actin depolymerization rate</td><td>10 s<sup>-1</sup></td></tr><tr><td>D<sub>a</sub></td><td>actin network diffusion coefficient</td><td>0.8 μm<sup>2</sup>/s</td></tr><tr><td>γ</td><td>tension coefficient</td><td>20 pN</td></tr><tr><td>ρ<sup>tot</sup><sub>a</sub></td><td>total density of actin</td><td>4 μm<sup>-2</sup></td></tr><tr><td>ρ<sub>a</sub></td><td>density of polymerized actin at the membrane</td><td>0.28 μm<sup>-2</sup></td></tr><tr><td>λ</td><td>effective (adhesion) friction coefficient</td><td>1–100 Pa. s/μm</td></tr></tbody></table>

Our (adhesion) friction coefficient should therefore be higher than this effective (pure) friction coefficient. We tried different orders of magnitude, and kept the most significant ones.

We end up with the parameters described in Table 1.

## 4.2. Choice of the initial shape

To run the simulations, we have to decide an initial shape for the membrane. Recall that we want to observe the duality of the shape of the lamellipodium, that is, we expect the shape to evolve into a compact shape which fulfills the micro-channel.

We first choose a simple initial shape, where the free membrane (the membrane which is not in contact with the walls) is a semi-ellipse. In addition to its simple implementation, this shape has the advantage of being already quite compact, therefore it is a good departure to observe whether compact shapes are stable, or whether the system will evolve to a less compact shape.

The second initial shape we choose is more realistic: it models the lamellipodium at its entrance in the micro-channel. A long part of the thin lamellipodium has already progressed on one wall (say, the lower wall $Γ_b$), dragging the cell body behind itself. The situation represented is when the cell body has just attained the entrance, and a tiny portion of the membrane has just adhered to the upper wall $Γ_d$. Therefore, the membrane has the shape of a stiff stair (or a letter "L"): the lower, long and thin stair is the thin lamellipodium, and the upper, short stair is the point where the membrane has newly adhered.

## 4.3. Numerical results for the simplified model (12)

We present in figures 5 and 6 six snapshots of the results of the simulations. What is represented is the value of the level set function, and the white line is its zero-value, that is, the membrane.

The first simulation represented in figure 5, deals with the movement of the lamellipodium into the micro channel with an initial shape corresponding to the first configuration (semi-ellipse). The associated computational domain $Ω$ is a rectangle of size $[0, 4μm]×[0, 1μm]$. The simulation is performed on a grid of size (512×128) corresponding to a space step $Δ_x = 7.8125·10^{-3}μm$.

The second simulation, illustrated in figure 6, represents the movement of the lamellipodium which just entered into the channel (as described in the last paragraph of Subsection 4.2). The computational domain is a rectangle of size $[0, 2μm] × [0, 1μm]$, and the corresponding grid resolution is (256 × 128). The interface thickness $ε$ is set to $2Δ_x$.
The first image of each sequence is the initial datum. In both cases, we observe that the cytoplasm tends to fulfill the micro-channels. For the first case, we observe an immediate progression of the membrane in the micro-channel. For the second case, the process can be decomposed in two phases: the cytoplasm first tends to change its shape into a convex shape, as compact as possible, then it starts to progress in the micro-channel. From what we observed of the numerical simulations performed for different sets of data (initial conditions and parameters), we can argue that this behaviour in two phases is quite typical of our model. For the case of the semi-circular initial datum, the initial shape being already compact, it is natural to expect that the system directly starts with the second phase.

Note that no progression is observed precisely at the wall: this is a direct consequence of the assumption of perfect adhesion, which prevents the fluid from any displacement at the contact with the wall. In particular, the membrane in contact with the wall cannot progress (or deadhere). However, the general movement of the free membrane indicates its tendency to progress in the micro-channel. Therefore, on one hand, we seem to observe a tendency to progress; on the other hand, one of the main limitations of our model is the assumption of perfect adhesion. In conclusion, we expect that removing the assumption of perfect adhesion would lead to observe exactly what is aimed, that is, a tendency of the membrane to remodel in a compact shape, followed by a progression of the cytoplasm in the micro-channel, including at the contact with the wall.

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Behaviour of the lamellipodium with having an initial shape of semi-ellipse. The background color shows the level set amplitude.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Behaviour of the lamellipodium having an initial shape of a letter "L". The background color shows the level set amplitude.</figcaption></figure>

## 4.4. Numerical simulations for the full model

We implement the model (11) with the finite element solver Freefem++. We chose the FE-spaces $\mathbb{P}^1$ for the pressure and $\mathbb{P}^1$-bubble for all the other functions. These spaces are classical in the implementation of Stokes equations in Finite Elements methods [3].

Firstly the friction boundary condition is classically taken into account in the variational formula of the Stokes problem.

Secondly, we now have to treat the reaction-transport-diffusion equation on $\rho_a$ on the moving domain $\Omega_{1,T}$. To circumvent this difficulty we approximate the equation (9) by an equation defined on the whole domain $\Omega_T$ while penalizing the Neumann boundary condition $\nabla\rho_a \cdot n = 0$ on the free membrane $\Gamma_1$. We still note $\rho_a$ the solution of this new equation:

$$ \partial_t \rho_a = -\nabla \cdot (\rho_a U) - \nabla \cdot (D_a \nabla \rho_a) + f(\rho_a, \rho_a^{cyt}) + \nabla \cdot \left( h_{\varepsilon}(\phi) \frac{\nabla \phi \otimes \nabla \phi}{\varepsilon |\nabla \phi|^2} \nabla \rho_a \right). $$

For the numerical simulations, we consider the test case of the lamellipodium which just enter into the channel (having an initial shape of a letter "L") and we test our model for two different friction values $\lambda_1 = 1$ and $\lambda_1 = 100$. The friction forces exered by the wall on the surrounding fluid are neglected ($\lambda_f = 0$), giving the following friction coefficient :

$$ \lambda(\phi) = \lambda_1(1 - H_{\varepsilon}(\phi)). $$

The computational domain $\Omega$ is a rectangle of size $[0, 1\mu m] \times [0, 4\mu m]$. We use the parameters of Table 1, except for the parameters $D_a$ and $\eta_a$ which are taken as follow: $D_a = 1$, $\eta_a = 100$. In these simulations, the space steps are $\Delta_x = \Delta_y = 6.25 \cdot 10^{-2}\mu m$, and the time step is $\Delta_t = 10^{-3}\text{s}$. We also chose the regularization parameter $\varepsilon$ to be $3\Delta_x$.

For completeness, the initial shape is obtained by the following FreeFem++ instructions so as to encode the "L" shape:

```matlab
func real heav0(real t){
    real out=0.5*(1+tanh(t));
    return out;
}

func real heavy(real t){
    real out=heav0(t/e);
```
```c
return out;
}
func real hauty(real t){
    real out=0.25*(heavy(1-t)-0.5)+0.375*(heavy(0.2-t));
    return out;
}
func phi0=tanh(hauty(x)-abs(y-hauty(x)));
```

Note that the variable *e* stands here for *ε*. These instructions will project onto the finite element space a function that is theoretically equal to 0 on the segments [(0, 1), (0.2, 1)], [(0.2, 1), (0.2, 0.25)], [(0.2, 0.25), (1, 0.25)], [(1, 0.25), (1, 0)] and [(1, 0), (0, 0)]. The actual shape is regularized by the use of regularized Heaviside functions heav0 and hauty so as to give in the end the shape that can be seen on Figures 7 and 8.

Figure 7 shows the simulations obtained in the case of a small friction force $\lambda_1 = 1$. As expected, the membrane slips along, attaches to and breaks away from the wall while moving forward into the channel. As the friction coefficient is small, the adhesion is quite weak. Consequently, the progression in the channel is quite fast.

Figure 8 shows the simulations obtained in the case of a large friction force $\lambda_1 = 100$, which is to be more realistic. Again, we observe a progression of the lamellipodium in the channel, including at the contact with the wall (through sliding and adhesion/de-adhesion). As the friction coefficient is high, the adhesion is strong. Consequently, in comparison with the previous case, the progression in the channel is slowed down. Notice indeed that the lengths of time required to fulfill the channel are not the same in the two cases.

The simulations obtained here are quite promising: as in the case with perfect adhesion we observe that the lamellipodium tends to remodel in a compact shape and to progress into the channel; but contrarily to the case with perfect adhesion, in this more realistic case the cell progresses also along the wall. These results motivate the enrichment of the main code we developed with realistic boundary conditions.

<figure><img src="image_8.png" /><figcaption>FIGURE 7. Behaviour of the lamellipodium having an initial shape of a letter "L". The friction coefficient $\lambda_1$ is set to 1. The background color represents the values of the regularized heaviside function $H_\epsilon$ applied to the level-set function: the red color corresponds to the lamellipodium and the blue object represents the surrounding fluid.</figcaption></figure>

## 4.5. Comparison with an unconfined case

In order to test the biological hypothesis that the shape of the lamellipodium is prescribed by the confinement, it is interesting to test our model in the case of a non confined domain. To do so, we use the full model implemented in FreeFem with a larger height for the computational domain $\Omega$. This models the situation where the lamellipodium moves in an environment whose dimensions are wider than its own size. We take $\Omega = [0, 4\mu m] \times [0, 4\mu m]$. All other parameters are the same as in the previous subsection (with $\lambda_1 = 100$). The initial shape for this case is given by the following set of instructions:
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Behaviour of the lamellipodium having an initial shape of a letter "L". The friction coefficient $\lambda_1$ is set to 100. The background color represents the values of the regularized heaviside function $H_\varepsilon$ applied to the level-set function: the red color corresponds to the lamellipodium and the blue object represents the surrounding fluid.</figcaption></figure>

```matlab
real x0=1;
real x1=0.25;
func prof=x1+(x0-x1)*(1-y);
func phi0=heav0((prof-x)/(hx))*(-100*(x+1)*(prof-x)*y*(1-y))+2*(1-heav0((prof-x)/(hx)));
```

See Figure 9. This initial condition is an analog of the "L" shape of the confined case in the sense that it also models the instant where the body of the cell has just started to enter the canal. While in the confined case the simulations start when the body has just touched both walls of the narrow canal, this situation is unlikely to happen in the wide canal. Note that for simulations we carried on with an initial "L" shape in this unconfined domain, the cell quickly moves back out of the canal, therefore we can not observe the cell deformation during a meaningful time.

We observe that the lamellipodium mainly displaces in the horizontal direction, going forward and backward. Even though after a long time, the lamellipodium has a light tendency to grow in the vertical direction, this growth cannot be considered as significant and may be a numerical artefact. Indeed, the vertical growth is around 15%, while the horizontal expansion reaches far above 100%. This is definitely not comparable to the clear and immediate tendency of fullfilling the domain that we observe in the confined case. However, its shape is ambiguous: on one hand, it does not fullfill the canal; but on the other hand, it does not reach the expected flat shape.

Therefore, while these tests tend to confirm our hypothesis that the ingredients chosen allow to recover the compact shape in a confined domain, the model does not seem to be adapted to represent an unconfined situation.

It is also interesting to notice that, as observed in some biological experiments, see [15], the lamellipodium may reach higher velocities under some experimental conditions in a confined environment.

## 5. CONCLUSION AND OUTLOOK

This proceeding is a report on a point of entry to study the effect of confinement on cell motility. It presents two prospective toy-models and associated numerical simulations that show that to understand the shape taken by the lamellipodium in a micro channel is a highly complex and challenging question. In this work, we wish to report that a few ingredients can lead, at least at the numerical level, to a behaviour resembling the “plug” form of cells in microchannels, reported in the experimental literature [24, 26]. This minimal model allows to recover qualitatively two main properties of the confined cell: its tendency to reshape compactly, and its actual progression inside the domain. It therefore suggests that the specificities of the cell motility in confined
<figure><img src="image_10.png" /><figcaption>FIGURE 9. Behaviour of the lamellipodium in the unconfined case. The background color represents the values of the regularized heaviside function $H_{\varepsilon}$ applied to the level-set function: the red color corresponds to the lamellipodium and the blue object represents the surrounding fluid.</figcaption></figure>

environment could originate from the same active mechanisms than the cell motility on a flat substrate, and not necessarily mechanisms that are specific to the confined domain.

However, the ambiguous behaviour observed for the model in an unconfined domain does not allow to conclude that the model in its present form explains the duality between flat and "plug" form. Therefore, it seems that at this stage the model is not appropriate to model unconfined situations. There are many possible reasons for that, lying in the several approximations we made in this extremely simplified model. We discuss some of them below.

As we wanted to keep the description at a macroscopic level, we chose to describe the lamellipodium as essentially its cellular liquid. The effect of the actin polymerization thus appears only in a protrusive force acting on this liquid. However it is a subject of debate (we have given a few of its terms in Section 1) whether this is a relevant assumption or whether the effect of the elasticity of the actin network should be taken a lot more into account in a more detailed mechanical model of the lamellipodium (again see for instance [21]). This is a really difficult question even for “passive” material (like soil for instance) and it leads to intricate models. As further developments, it will be crucial to determine if this description is sufficient to see the duality of flatness/“plugness” we wish to show.

The description of “adhesion” as being a friction of the water molecules on the solid boundaries can largely be discussed. In particular, the impact of adhesion in terms of normal forces (to the boundary) was here deliberately neglected. Of course, more sophisticated modelling will be investigated in further work. Overall, boundary conditions at the walls or at the back of the lamellipodium can also be discussed and rework in general to accomplish a finer description of the system.

Finally, remark that the well-posedness of both systems (complete or simplified) is a non trivial issue that has to be addressed in the future. Similarly, more relevant numerical methods should be tested for the complete
model. This would allow to test a larger set of boundary conditions. Even at the level of the current code, the (numerical) stability would need to be studied in depth to check what is the domain of parameters (numerical or from modelling) that can be used, and also the stability of our preliminary results with respect to initial shape.

## REFERENCES

[1] M Abercrombie, JEM Heaysman, and SM Pegrum. The locomotion of fibroblasts in culture: II. movements of particles on the dorsal surface of the leading lamella. *Experimental cell research*, 62(2):389-398, 1970.

[2] VC Abraham, V Krishnamurthi, DL Taylor, and F Lanni. The actin-based nanomachine at the leading edge of migrating cells. *Biophysical journal*, 77(3):1721-1732, 1999.

[3] F Brezzi and J Pitkäranta. *On the Stabilization of Finite Element Approximations of the Stokes Equations*, pages 11–19. Vieweg+Teubner Verlag, Wiesbaden, 1984.

[4] RM Capito and M Spector. Scaffold-based articular cartilage repair. *Engineering in Medicine and Biology Magazine, IEEE*, 22(5):42–50, 2003.

[5] R Chatelin. *Méthodes numériques pour l'écoulement de Stokes 3D: fluides à viscosité variable en géométrie complexe mobile; application aux fluides biologiques*. PhD thesis, Université Paul Sabatier-Toulouse III, 2013.

[6] R Chatelin and P Poncet. Hybrid grid-particle methods and penalization: A sherman-morrison-woodbury approach to compute 3d viscous flows using fft. *Journal of Computational Physics*, 269:314-328, 2014.

[7] AJ Chorin. Numerical solution of the navier-stokes equations. *Mathematics of computation*, 22(104):745-762, 1968.

[8] GH Cottet and E Maitre. A level-set formulation of immersed boundary methods for fluid-structure interaction problems. *Comptes Rendus Mathematique*, 338(7):581-586, 2004.

[9] GH Cottet and E Maitre. A level set method for fluid-structure interactions with immersed surfaces. *Mathematical models and methods in applied sciences*, 16(03):415-438, 2006.

[10] P Friedl, F Entschladen, C Conrad, B Niggemann, and KS Zänker. Cd4+ t lymphocytes migrating in three-dimensional collagen lattices lack focal adhesions and utilize β1 integrin-independent strategies for polarization, interaction with collagen fibers and locomotion. *European journal of immunology*, 28(8):2331-2343, 1998.

[11] P Friedl and B Weigelin. Interstitial leukocyte migration and immune function. *Nature immunology*, 9(9):960-969, 2008.

[12] WS Haston, JM Shields, and PC Wilkinson. Lymphocyte locomotion and attachment on two-dimensional surfaces and in three-dimensional matrices. *The Journal of cell biology*, 92(3):747-752, 1982.

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

[14] M Herant and M Dembo. Form and function in cell motility: from fibroblasts to keratocytes. *Biophysical Journal*, 98(8):1408-1417, 2010.

[15] J Jacobelli, RS Friedman, MA Conti, AM Lennon-Dumenil, M Piel, CM Sorensen, RS Adelstein, and MF Krummel. Confinement-optimized three-dimensional t cell amoeboid motility is modulated via myosin iia-regulated adhesions. *Nature immunology*, 11(10):953-961, 2010.

[16] T Lämmermann, BL Bader, SJ Monkley, T Worbs, R Wedlich-Söldner, K Hirsch, M Keller, R Förster, DR Critchley, R Fässler, et al. Rapid leukocyte migration by integrin-independent flowing and squeezing. *Nature*, 453(7191):51-55, 2008.

[17] A Manhart, C Schmeiser, N Sfakianakis, and D Oelz. An extended filament based lamellipodium model produces various moving cell shapes in the presence of chemotactic signals. *arXiv preprint arXiv:1502.02442*, 2015.

[18] Y Mori, A Jilkine, and L Edelstein-Keshet. Wave-pinning and cell polarity from a bistable reaction-diffusion system. *Biophysical Journal*, 94:3684-3697.

[19] S Osher and JA Sethian. Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations. *Journal of computational physics*, 79(1):12-49, 1988.

[20] E Rouy and A Tourin. A viscosity solutions approach to shape-from-shading. *SIAM Journal on Numerical Analysis*, 29(3):867-884, 1992.

[21] B Rubinstein, MF Fournier, K Jacobson, AB Verkhovsky, and A Mogilner. Actin-myosin viscoelastic flow in the keratocyte lamellipod. *Biophys. J.*, 97:1853-1863, 2009.

[22] E Sahai. Illuminating the metastatic process. *Nature Reviews Cancer*, 7(10):737-749, 2007.

[23] D Shao, H Levine, and WJ Rappel. Coupling actin flow, adhesion, and morphology in a computational cell motility model. *Proceedings of the National Academy of Sciences*, 109(18):6851-6856, 2012.

[24] K Wilson, A Lewalle, M Fritzsche, R Thorogate, T Duke, and G Charras. Mechanisms of leading edge protrusion in interstitial migration. *Nature communications*, 4, 2013.

[25] IV Yannas, E Lee, DP Orgill, EM Skrabut, and GF Murphy. Synthesis and characterization of a model extracellular matrix that induces partial regeneration of adult mammalian skin. *Proceedings of the National Academy of Sciences*, 86(3):933-937, 1989.
[26] AK Yip, KH Chiam, and P Matsudaira. Traction stress analysis and modeling reveal amoeboid migration in confined spaces is accompanied by expansive forces and requires the structural integrity of the membrane-cortex interactions. *Integrative Biology*, 2015.