Analysis of a non-local and non-linear Fokker-Planck model for cell crawling migration

TL;DR

This paper introduces and analyzes a non-local, non-linear Fokker-Planck model for cell crawling migration, capturing the universal coupling between cell speed and persistence, with results on solution behavior and convergence.

Contribution

It presents a novel non-local, non-linear Fokker-Planck model for cell motility, including analysis of solution existence, blow-up conditions, and convergence properties.

Findings

Solutions are global below critical mass

Solutions blow up above critical mass

Quantitative convergence to equilibrium

Abstract

Cell movement has essential functions in development, immunity and cancer. Various cell migration patterns have been reported and a general rule has recently emerged, the so-called UCSP (Universal Coupling between cell Speed and cell Persistence), [30]. This rule says that cell persistence, which quantifies the straightness of trajectories, is robustly coupled to migration speed. In [30], the advection of polarity cues by a dynamic actin cytoskeleton undergoing flows at the cellular scale was proposed as a first explanation of this universal coupling. Here, following ideas proposed in [30], we present and study a simple model to describe motility initiation in crawling cells. It consists of a non-linear and non-local Fokker-Planck equation, with a coupling involving the trace value on the boundary. In the one-dimensional case we characterize the following behaviours: solutions are global if the mass is below the critical mass, and they can blow-up in finite time above the critical mass. In addition, we prove a quantitative convergence result using relative entropy techniques.

Figures

Equations

$$\partial_{t}c(t,x)=\mathrm{div}\;\Big{(}\nabla p(t,x)c(t,x)+\nabla c(t,x)\Big{)}\,,\ t>0\,,\ x\in\Omega\,,$$

$$\{\begin{array}{ll}-\Delta p=-1,& \text{ on }\Omega ,\\ p=1-c,& \text{ on }\partial \Omega .\end{array}$$

$$v(t)=-{\int }_{\partial \Omega }c(t,y)\mathbf{n}\mathrm{d}\sigma ,t>0,$$

$$\partial_{t}c(t,x)=\partial_{xx}c(t,x)+\partial_{x}\Big{(}\big{(}x+c(t,-1)-c(t,1)\big{)}c(t,x)\Big{)}\,,\quad t>0\,,\ x\in\Omega\,,$$

$$\begin{cases}\big{(}c(t,-1)-c(t,1)-1\big{)}c(t,-1)+\partial_{x}c(t,-1)=0\,,\\ \big{(}c(t,-1)-c(t,1)+1\big{)}c(t,1)+\partial_{x}c(t,1)=0\,.\end{cases}$$

$$\underset{t\in (0,T)}{\sup }{\int }_{-1}^{1}c(t,x)\log c(t,x)\mathrm{d}x$$

$${\int }_0^T{\int }_{-1}^{1}c(t,x){\left({\partial }_x\log c(t,x)\right)}^2\mathrm{d}x\mathrm{d}t$$

$$\{\begin{array}{ll}\mathrm{div}(\tilde{U})=-\frac{1}{\xi }\Delta \tilde{P}=-k_d& \text{mbox}in\Omega (t),\\ \tilde{P}=k_p& \text{mbox}on\partial \Omega (t).\end{array}$$

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

$$\sigma =\mu \left(\nabla \tilde{U}+{}^t\nabla \tilde{U}\right)-\tilde{P}\text{ Id},$$

$$f=-\xi \tilde{U},$$

$$\nabla \tilde{P}=-\xi \tilde{U}\text{ in }\Omega (t).$$

$${\partial }_t\rho +\mathrm{div}(\rho \tilde{U})=-k_d\rho ,$$

$$\rho ={\rho }_0+\epsilon k_p\text{ on }\partial \Omega (t),$$

$$\tilde{P}=\frac{1}{\epsilon }(\rho -{\rho }_0),$$

$${\partial }_t\rho -\frac{1}{\xi }\mathrm{div}\left(\rho \nabla \left(\frac{\rho }{\epsilon }\right)\right)$$

$$\rho$$

$$V(t)=-\gamma {\int }_{\Omega (t)}\tilde{U}(t,X)\mathrm{d}X.$$

$${\partial }_t\tilde{C}(t,X)+\mathrm{div}\left(\tilde{U}(t,X)\tilde{C}(t,X)-D\nabla \tilde{C}(t,X)\right)=0\text{ for }X\in \Omega (t),t>0,$$

$$\left(D\nabla \tilde{C}(t,X)-\tilde{C}(t,X)\tilde{U}(t,X)\right)\cdot \mathbf{n}=0\text{ for }X\in \partial \Omega (t),t>0,$$

$$\{\begin{array}{ll}\nabla \tilde{P}(t,X)=-\xi \tilde{U}(t,X)& \text{ for }X\in \Omega (t),t>0,\\ \mathrm{div}\tilde{U}(t,X)=-K_d& \text{ for }X\in \Omega (t),t>0,\\ \tilde{P}(t,X)=f\left(\tilde{C}(t,X)\right):=\alpha -\beta \tilde{C}(t,X)& \text{ for }X\in \partial \Omega (t),t>0,\end{array}$$

$$\tilde{P}(t,X)=\Big{(}\alpha-\beta\tilde{C}(t,X)\Big{)}_{+}\text{ for }X\in\partial\Omega(t),\,t>0\,,$$

$$V_2(t)$$

$$\begin{array}[]{lrcl}L(.,t):&\Omega_{0}&\longrightarrow&\Omega(t)\\ &x&\mapsto&X=x+\int_{0}^{t}V(s)\,{\mathrm{d}}s=L(t,x),\end{array}$$

$${\partial }_{t}C(t,x)={\partial }_t\tilde{C}(t,L(t,x))+\underset{=V(t)}{\underbrace{{\partial }_{t}L(t,x)}}{\nabla }_X\tilde{C}(t,X)=\frac{\mathrm{D}\tilde{C}}{\mathrm{D}t},$$

$$\left\{\begin{array}[]{rclr}\displaystyle\frac{{\,\rm D}}{{\,\rm D}t}\tilde{C}(t,X)+\mathrm{div}\;\left((\tilde{U}(t,X)-V(t))\tilde{C}(t,X)-D\nabla\tilde{C}(t,X)\right)&=&0&\,X\in\Omega(t),\,t>0,\\ \displaystyle\left(D\nabla\tilde{C}(t,X)-\tilde{C}(t,X)(\tilde{U}(t,X)-V(t))\right)\cdot\mathbf{n}&=&0&\,X\in\partial\Omega(t),\,t>0.\end{array}\right.$$

$$\left\{\begin{array}[]{rcll}\displaystyle\partial_{t}C(t,x)+\textrm{div }\left[(U(t,x)-V(t))C(t,x)-D\nabla C(t,x)\right]&=&0&\,\text{ for }x\in\Omega_{0},\,t>0,\\ \displaystyle\left(D\nabla C(t,x)-C(t,x)(U(t,x)-V(t))\right)\cdot\mathbf{n}&=&0&\,\text{ for }x\in\partial\Omega_{0},\,t>0.\end{array}\right.$$

$$\{\begin{array}{ll}\nabla p=-u& \text{ on }{\Omega }_0,\\ \mathrm{div}u=-\frac{\xi 4R^2}{\alpha }{K}_d=:-k_d& \text{ on }{\Omega }_0,\\ p=1-\delta c& \text{ on }\partial {\Omega }_0,\end{array}$$

$$v(t)=-\gamma {\int }_{{\Omega }_0}u\mathrm{d}x.$$

$$\{\begin{array}{ll}{\partial }_{t}c+\mathrm{div}\left((u-v)c-D^{\prime }\nabla c\right)=0& \text{ on }{\Omega }_0,\\ \left(D^{\prime }\nabla c-(u-v)c\right)\cdot \mathbf{n}=0& \text{ on }\partial {\Omega }_0,\end{array}$$

$${\int }_{{\Omega }_0}c\mathrm{d}x=M.$$

Full text

1 Introduction

Cell migration is a fundamental biological process involved in morphogenesis, tumor spreading, and wound healing [39, 23, 18]. One of its most spectacular instance is cell crawling, which is crucial to immune cells in order to reach an inflammation spot, but is also observed e.g. in tumor cells during metastasis formation [4]. Identifying key mechanisms involved in cell migration is then a major issue both for our fundamental understanding and for clinical research.

We focus here on cells crawling on a flat substrate, without any external cue. This type of movement occurs in several steps: protrusion of the cell at the ”leading edge” (led by actin polymerisation), adhesion to the substrate, contraction and translocation of the cell body. However, if the notion of ”leading edge” is clear for cells undergoing persistent motion (as in the chemotaxis phenomenon), the relation between protrusive activity and directionnality in more general cases is still unresolved. Indeed, migrating cells might transiently polarize - a front and a back appear. If such ability is weak, the resulting motion is close to random. This polarity is reflected at the molecular level by a restriction of certain molecules to particular regions of the inner cell membrane. For example, the activated forms of Rac and Cdc42 molecules are found at the front of the cell, whereas RhoA molecules are found toward the rear (see e.g [30]). It is known that this asymmetry results from strong feedbacks exerted between these signaling pathways and mechanical elements [14]. Self-polarisation is then a multiscale phenomenon for which modeling efforts can bring new insights [33].

Integrated models has been succesfully developped for addressing the issue of motile cell shape, in particular for the keratocyte [40]. Computational models allow for in silico experiments and direct comparison with experimental measures (see [24, 41], and more generally [17]). In particular, it has been highlighted in [41] that redundant minimal mechanisms exist to ensure a motile stable shape for the keratocyte.

Self-polarisation has led to computational models investigating the role of contraction, adhesion and actin polymerisation for keratocytes or epithelial cells, providing useful predictions and allowing confrontation with experimental data [37, 1, 29]. Other mechanochemical models were built to perform in silico experiments ([12, 32]), but they remain unable to provide a minimal framework to explain the persistence issue in cell motion.

Recently, in [31], it was shown experimentally that cell persistence is correlated to the flow of actin filaments from the front to the back of the motile part of the cell, as protrusions with faster actin flows are more stable in time. It was also shown that faster actin flows generate steeper gradients of actin-binding proteins.

In this work, we study a conceptually minimal model for polarity initiation and maintenance based on an active gel description of the actin cytoskeleton, and its interaction with molecular signaling pathways. More precisely, we use a rigid version of the model first proposed in [2] that we enrich, in the spirit of [31], with a feedback loop between actin polymerisation locations and polarity markers [30, 26], that are advected by a dynamic actin cytoskeleton undergoing flows at the cellular scale. This minimal model allows us to obtain qualitative results on cell persistence.

We analyse the long time asymptotics of the model. The principal goal of our analysis is to identify regimes in which non homogeneous stationary states, that will be interpreted as polarized states, emerge.

The main ingredients of the model are as follows. We assume that there exists a mesoscopic length scale, small compared to the cell but large compared to individual molecules at which the properties of the cytoskeleton and of the solvent, which constitutes the cell, can be described by continuous fields [20]. Following [21, 22, 2], we describe the actin cytoskeleton as a 2D Darcy flow, bounded by a membrane with a given shape. We assume that the relevant dynamics is sufficiently slow to neglect elastic effects. Actin is assumed to polymerize at the membrane and to depolymerize uniformly at a constant rate in the interior of the cell. Moreover, polymerisation is assumed to depend on the concentration of a marker that is advected by the actin flow itself. Finally, since we focus here on cell crawling on a flat substrate, the main external force arising is friction with the substrate.

The markers, whose density is denoted $c(t,x)$, are assumed to diffuse in the cytoplasm and to be actively transported along the cytoskeleton. Consider a viscous active fluid with pressure $p$ filling a two-dimensional bounded domain of fixed shape figuring the cell to describe the cytoskeleton. When all coefficients are set to 1, and after a global change of variable so that the domain of definition $\Omega$ is fixed, the resulting motion is a biased diffusion equation with advection field in the cell frame $u(t,x)=-\nabla p(t,x)$:

[TABLE]

together with zero-flux boundary condition on $\partial\Omega$ in order to conserve the molecular content. The active character of the cytoskeleton is contained in the pressure problem. Actin polymerization at the boundary leads to a Dirichlet condition depending on $c$, while depolymerization arises as a sink term inside the domain. More precisely, the pressure is assumed to satisfy, for all $t>0$:

[TABLE]

Furthermore, considering a global friction coefficient set to $1$, the cell velocity $v$, which arises from the inside flow rubbing on the substrate, is given by

[TABLE]

where $\mathbf{n}$ denotes the outward unit normal and ${\,\rm d}\sigma$ the surface measure on $\partial\Omega$.

In the one-dimensional case, assuming that $\Omega=(-1,1)$ the previously described model writes as a non-linear and non-local Fokker-Planck equation:

[TABLE]

together with zero-flux boundary conditions at $x=-1$ and $x=1$:

[TABLE]

The main aim of this paper is to provide some results on the long time asymptotics of the solution to (4) - (5) and to give estimates on the convergence rates in cases of convergence. To this end, we first look for stationary states. Let us denote $M=\int_{\Omega}c(t,x)\,{\mathrm{d}}x$ the (conserved) mass of molecules inside the cell. In the case $M\leq 1$, there exists a unique stationary state $G_{M}(x):=Mexp(-x^{2}/2)/\int_{-1}^{1}exp(-y^{2}/2){\,\rm d}y$ towards which the solution converges and the asymptotic result is obtained through the convergence to zero of a suitable Lyapunov functional $\mathbf{L}$ defined in (29) in the case $M<1$, and through the convergence to zero of the relative entropy $\mathbf{H}$ defined in (28) in the case $M=1$. Moreover, in both cases, using the logarithmic Sobolev inequality with a suitable function, we obtain an exponential decay to equilibrium.

Theorem 1**.**

Assume that the initial datum $c_{0}$ satisfies both $c_{0}\in L^{1}(-1,1)$ and $\int_{-1}^{1}c_{0}(x)\log c_{0}(x){\,\rm d}x<+\infty$. Assume in addition that $M\leq 1$, then there exists a global weak solution of (4) - (5) that satisfies the following estimates for all $T>0$,

[TABLE]

*Moreover, the solution strongly converges in $L^{1}$ towards the unique stationary state $G_{M}(x)$ and the rate is exponential. *

Solutions of (4) may become unbounded in finite time (so-called blow-up). This occurs if the mass M is above the critical mass: $M>1$, and for an asymmetric enough initial profile, measured with the first momentum shifted in $x=-1$ of $c$, $\tilde{\mathbf{J}}(t)=\int_{-1}^{1}(x+1)c(t,x){\,\rm d}x$.

Theorem 2**.**

Assume $M>1$ and the initial first moment shifted of 1 is small: $\tilde{\mathbf{J}}(0)<(M-1)/2$. Assume in addition that $c_{0}$ satisfies $c_{0}(-1)-c_{0}(1)>1$. Then the solution to (4) - (5) with initial data $c(0,x)=c_{0}(x)$ blows-up in finite time.

In Figure 1, we illustrate Theorems 1 and 2.

Although in the present biological context, blow-up of solutions is interpreted as cell polarisation, such a blow-up in finite time might be questionable. Indeed a strong instability drives the system towards an inhomogeneous state and blow-up corresponds to the case where the drift becomes infinite. The boundary condition (5) which is responsible for infinite drift turns out to be unrealistic from a biophysical viewpoint. On the way towards a more realistic model, we distinguish between cytoplasmic content $c(t,x)$ and the concentration of trapped molecules on the boundary at $x=\pm 1$ : $\mu_{\pm}(t)$. Then the exchange of molecules at the boundary is described by very simple kinetics.

The plan of this work is as follows. In Section 2, we introduce the model. In Section 3, we analyse with full details the one-dimensional case. In Section 4, we briefly study a model with exchange of markers at the boundary in the one-dimensional case.

2 The model

Our purpose in this section is to derive the model given by equations (1) - (2) - (3) which describes the behaviour of an active viscous fluid featuring the cytoskeleton.

2.1 Main constitutive assumptions

We consider a two-dimensional layer of viscous fluid, representing the cell cytoskeleton, surrounded by a rigid membrane. Following [21, 22, 2], to model actin polymerisation and depolymerization, we add active properties to the viscous fluid. The first active property we consider is the out-of-equilibrium polymerization of the fluid. In a cell, actin monomers are added to actin filaments by the consumption of the biological fuel ATP. Other proteins regulate the nucleation and polymerization of actin filaments. It is commonly observed that actin polymerization activators such as WASP proteins preferentially locate along the cell membrane [38]. For this reason we suppose that the fluid is polymerised at the membrane. The main novelty of our work relies on the coupling between actin polymerization and a biological marker which is transported by actin flows. Its aggregation in a part of the membrane characterizes the rear of the cell, hence its polarisation. This marker could be an antagonist to polymerization-inducing molecules (Rac1, Cdc42), such as RhoA, Arpin, or even myosin II (see [10, 25]). Furthermore polymerization is balanced by depolymerization, which we assume to occur uniformly at a constant rate in the cell body, to ensure the renewal of ressources for polymerisation. polymerization and depolymerization induce an inward flow which rubs on the substrate. This friction is responsible for the cell displacement.

More precisely, we consider an active Darcy flow, which models the actin cytoskeleton, inside a (moving) domain $\Omega(t)\subset\mathbb{R}^{2}$, where $\Omega(0)$ is a disc of radius $R>0$. The domain moves by translation with a velocity $V(t)\in\mathbb{R}^{2}$. We introduce the fluid pressure $\tilde{P}(t,X)$, and $\tilde{U}(t,X)$ the actin filaments velocity. Finally, we denote $\tilde{C}(t,X)$ the concentration of molecular marker. All these quantities are defined for $X\in\Omega(t)$ and $t>0$. We assume that the fluid density $\rho(t,X)\equiv\rho$ is constant in time and space. Then, the fluid problem writes

[TABLE]

The boundary condition accounts for polymerization at the edge of the cell, while $-k_{d}$ describes depolymerization in the cell body.

Remark 1**.**

Note that equation (6) is similar but different from the model first introduced in [2]. Indeed in our case the polymerization is modeled through a pressure term and not a velocity one.

In the limit of low Reynolds number, viscous forces dominate over inertial forces and the Navier-Stokes equation simplifies to the force balance principle:

[TABLE]

where $\sigma$, the stress tensor, is given by

[TABLE]

with $\mu$ being the viscosity. Since the layer is placed on a substrate, we consider a friction force

[TABLE]

where $\xi$ is an effective friction coefficient.

Following [5], we neglect viscosity arising from the polymer-polymer and polymer-solvent friction forces and consider the limit $\mu\to 0$. This reduces to

[TABLE]

Notice that the pressure, hence the friction force, acts for the integrity of the fluid (no phase separation).

Remark 2**.**

If we were considering a deformable cell, see [2], it would have been natural to assume that the pressure $\tilde{P}$ satisfies $\tilde{P}=-\gamma\kappa$ in $\partial\Omega(t)$, where $\kappa$ is the curvature of $\partial\Omega(t)$ and $\gamma\geq 0$ is the superficial tension of the cell membrane. Here we consider a caricatural situation since we assume that the cell domain is non-deformable and we do not impose any additional condition on $\tilde{P}$ on the boundary of $\Omega(t)$.

Remark 3**.**

About the polymer density: this model can be derived from a more realistic one, where polymerization and depolymerization locally modify the polymer density. More precisely, write

[TABLE]

where $\Omega(t)=\{\rho(t,\cdot)>0\}$ denotes the region occupied by the cytoskeleton. To account for polymerization that takes place at the edge of the cell and which consists in a local increase in actin concentration, we impose on the cell membrane a jump on the actin concentration:

[TABLE]

where $\varepsilon$ is a small parameter and $\rho_{0}$ a constant. One can also assume that the polymerization rate is a continuous non-zero polar function in an annulus initially defined as $\overline{\Omega}(0)\setminus B(0,R-\lambda)$, of width $\lambda>0$, see e.g [27].

As it is classical, see [6], in addition, we assume that $\tilde{P}$ is a function of $\rho$:

[TABLE]

we then get the following problem:

[TABLE]

Formally the limit $\varepsilon\to 0$ of the previous model leads to the Poisson problem (6). This amounts to saying that the osmotic pressure $\tilde{P}$ ensures at all time that the polymer density stays constant. The rigorous justification of this limit is not our purpose here, and will be the object of a future work.

Finally, we have to prescribe the domain velocity arising from friction forces. Since the gel layer is at mechanical equilibrium, the cell moves as a consequence of the inside flow rubbing on the substrate, hence we will consider a moving domain. Friction forces occur at the microscopic scale, and mesoscopic tension forces also are at play. However, for the sake of simplicity, we will neglect this heterogeneity to consider a global friction coefficient $\gamma$. We write for all $t>0$

[TABLE]

The main novelty of our work is to consider that $k_{p}$ is a function of $c$, the concentration of a biological marker. This marker is assumed to diffuse and to be transported by the actin filaments, modeled by the previously described fluid. At the boundary, we prescribe a zero flux condition to ensure mass conservation. The corresponding problem writes

[TABLE]

with the zero-flux boundary condition

[TABLE]

where $\mathbf{n}$ is the outward unit normal to the boundary. The zero flux boundary condition ensures that the total mass is preserved: $\frac{\,{\mathrm{d}}}{\,{\mathrm{d}}t}\int_{\Omega(t)}\tilde{C}(t,X)\,{\mathrm{d}}X=0$.

Remark 4**.**

The non-deformability of $\Omega(t)$ is an important downside of the model, that conceals fundamental modelization issues. Indeed, the cell velocity is defined globally, preventing the description of local effects: experimentally, it is observed that the activity of the cytoskeleton deforms the membrane, leading to a displacement while exerting a geometric feedback on actin flows.

Notice that our choice for the friction force amounts to consider friction arising from both the retrograde flow and the displacement, but their effects on the equations are the same in 1D. Hence, the role of motion in the initiation and maintenance of polarisation cannot be investigated properly, but this will be studied in future works in 2D and for a free-boundary model.

In the spirit of [31], we will assume that polymerization occurs more likely at the boundary points where the marker concentration is the lower (the marker we consider is a rear marker).

Recalling what was explained previously, the fluid problem writes

[TABLE]

Remark 5**.**

A more realistic boundary condition on $\tilde{P}$ would have been to impose that

[TABLE]

where $(\cdot)_{+}$ denotes the positive part. Indeed the case where $\tilde{P}$ takes negative values is not realistic from the modelling viewpoint.

Remark 6**.**

Another option for the cell velocity would be to choose

[TABLE]

However this latter choice is less relevant from a biological viewpoint, since the global cell displacement is a consequence of the actin flow rubbing on the substrate.

Problem formulation on a fixed domain

We now define the problem on a fixed domain $\Omega_{0}:=\Omega(0)$. In such a case the domain motion appears in the problem formulation. As $(\tilde{P},\tilde{U},\tilde{C})$ are functions on $\Omega(t)$, we denote $(P,U,C)$ their analogous on $\Omega_{0}$. Moreover we define the following map:

[TABLE]

which gives $C(t,x)=\tilde{C}(t,L(t,x))$ (for example for $C$).

Let us now rewrite the fluid problem on $\Omega_{0}$. To do so we observe that for all $x\in\Omega_{0}$ and for all $t>0$:

[TABLE]

where $\frac{{\,\rm D}}{{\,\rm D}t}$ denotes the total derivative. The convection-diffusion equation rewrites as

[TABLE]

Moreover, the Jacobian matrix related to $L$ is the identity matrix (as $V$ is space-independent), leading to $\nabla_{X}\tilde{C}(t,X)=\nabla_{x}C(t,x)$. This, combined with $\tilde{U}(t,X)=U(t,x)$, yields that

[TABLE]

We can check again mass conservation: $\frac{\,{\mathrm{d}}}{\,{\mathrm{d}}t}\int_{\Omega_{0}}C\,{\mathrm{d}}x=0$.

Nondimensionalization

Let us introduce the following typical quantities for nondimensionalization: $\overline{L}=2R$, where $R$ is the radius of $\Omega_{0}$, $\overline{P}=\alpha=\frac{\overline{F}\,\overline{L}^{2}}{\overline{L}}$, $\overline{F}=\frac{\alpha}{\overline{L}}=\frac{\alpha}{2R}=\xi\overline{V}$, $\overline{V}=\frac{\alpha}{2R\xi}$, $\overline{T}=\frac{\overline{L}}{\overline{V}}=\frac{4R^{2}\xi}{\alpha}$ and $\overline{C}=\overline{L}^{-2}$.

Remark 7**.**

It has to be noticed that the force we consider is actually a force per unit of area. This explains the expression used for the typical pressure $\overline{P}$, corresponding to a force divided by a length (in 2D).

Remark 8**.**

The typical pressure corresponds to a maximal incoming pressure of the fluid at the boundary of the domain.

Now, write $(p,u,c,v)$ the nondimensionalized analogous to $(P,U,C,V)$. The nondimensionalized problem writes:

[TABLE]

where $\delta:=\frac{\beta\overline{C}}{\alpha}$ and the domain velocity is:

[TABLE]

The convection-diffusion problem is

[TABLE]

where we have introduced $D^{\prime}=\frac{D}{\overline{L}\,\overline{V}}=\frac{D\xi}{\alpha}$. Moreover, the global mass is prescribed:

[TABLE]

2.2 The one-dimensional case

The first equation of (14) rewrites

[TABLE]

and using (15), it leads to

[TABLE]

Moreover, differentiating (18) and using that $\partial_{x}u=-k_{d}$, we find $\partial_{xx}p(t,x)=k_{d}$ and we can write

[TABLE]

where

[TABLE]

Now, replacing $p$ in the expression of $u$ this gives

[TABLE]

Finally recalling the boundary conditions in (14), we obtain the following non-linear and non-local convection-diffusion equation

[TABLE]

with a velocity $v-u$ given by: for all $x\in(a,b)$ and for all $t>0$

[TABLE]

and the domain velocity:

[TABLE]

with $\delta=\frac{\beta\overline{C}}{\alpha}>0$ and we recall that $D^{\prime}=\frac{D\xi}{\alpha}\geq 0$.

Remark 9**.**

Note that up to a constant we recover the cell velocity $v_{2}$, see Remark 6.

From now on, in order to simplify the computations we assume that $\Omega=(-1,1)$. In such a case the velocity simply rewrites as

[TABLE]

Finally, taking $k_{d}=D^{\prime}=1$ and $\delta\left(\gamma+\frac{1}{2}\right)=1$, we obtain the problem (4)-(5).

We find a model that presents some similarities with the model studied in [16, 7, 8, 34, 28] to describe yeast cell polarisation. The main difference here is the presence of an additional drift towards the cell center.

3 The boundary non-linear Fokker-Planck equation in dimension 1

In this part we study the non-local and non-linear Fokker-Planck equation (4) - (5) that we recall now:

[TABLE]

and we prove Theorems 1 and 2.

We are looking for a solution to (24). As it is classical we start by giving a proper definition of weak solutions, adapted to our context:

Definition 3**.**

We say that $c(t,x)$ is a weak solution of (24) on $(0,T)$ if it satisfies:

[TABLE]

and $c(t,x)$ is a solution of (24) in the sense of distribution in $\mathcal{D}^{\prime}(-1,1)$.

Since the flux $\partial_{x}c(t,x)+\left(c(t,-1)-c(t,1)+x\right)c(t,x)$ belongs to $L^{1}\left((0,T)\times(-1,1)\right)$, the solution is well-defined in the distributional sense under assumption (25). In fact we can write

[TABLE]

Moreover, one has

[TABLE]

Note also that to prove existence of weak solutions in the sense of Definition 3, one should perform a regularization procedure. In this work we focus on long-time dynamics and we will not detail such a regularization procedure which is classical, we refer to [8] for more details.

Let us now observe that the non-negativity of a solution is preserved. Indeed if $c$ is solution in $L^{1}_{x}$, since $\mbox{sgn}(c)\partial_{xx}^{2}c\leq\partial_{xx}^{2}|c|$, then

[TABLE]

This proves that if $|c_{0}|=c_{0}$ almost everywhere (initial data non-negative) then $|c|=c$ almost everywhere for later times.

Moreover, weak solutions in the sense of Definition 3 are mass-preserving: $M=\int_{-1}^{1}c_{0}(x){\,\rm d}x=\int_{-1}^{1}c(t,x){\,\rm d}x$. A simple computation on the first momentum defined by $\mathbf{J}(t)=\int_{-1}^{1}xc(t,x){\,\rm d}x$, leads to

[TABLE]

Note that $\mathbf{J}(t)\geq-M$ since $\int_{-1}^{1}(x+1)c(t,x){\,\rm d}x\geq 0$. Here, similarly to [7, 8], we will prove that the following dichotomy occurs:

• •

when $M\leq 1$, the linear Fokker-Planck part drives the equation and the solution converges toward the unique stationary state, see sections 3.1 and 3.2,

• •

when $M>1$, the equation admits three stationary states, two of them being non-homogeneous. Moreover, for an asymmetric enough initial condition, the solution blows up in finite time, see section 3.3.

3.1 Global existence and asymptotic analysis for sub-critical mass $M<1$

In this section we prove Theorem 1 in the case $M<1$. The proof is decomposed in three steps. First we compute the stationary state of (24), then we build a Lyapunov functional, this allows establishing a-priori estimates. Finally we investigate long-time behaviour of solutions in the case $M<1$ using entropy methods. We stress out that the method for proving global existence in this case strongly relies on the existence of a Lyapunov functional, Lemma 5. This is why we analyse the global existence and the long time behaviour simultaneously. Note that the proof is very close to the one given in [8] and in [28]. The differences concern the expression of the stationary state and the domain geometry. In [8, 28], the stationary state is the family $\{\alpha\exp\left(-\alpha x-x^{2}/2\right)\}_{\alpha}$, and the domain is the half-line.

Let us start with the existence of a unique stationary state.

Lemma 4**.**

If $M<1$ equation (24) admits a unique stationary state given by $G_{M}=M\exp(-x^{2}/2)/\int_{-1}^{1}\exp(-y^{2}/2){\,\rm d}y$.

Proof.

An easy computation shows that any stationary state for (24) is either $G_{M}$ (which is symmetric) or of the form $G_{\alpha}$ for $\alpha=G_{\alpha}(-1)-G_{\alpha}(1)\neq 0$ to be found. For $\alpha>0$, it writes

[TABLE]

It remains to find $\alpha$ such that the mass constraint $\int_{-1}^{1}G_{\alpha}(x){\,\rm d}x=M$ is satisfied. This rewrites $P(\alpha)=M$, $P$ being the function defined by:

[TABLE]

We observe that $P(\alpha)>\int_{0}^{2\alpha}\frac{1}{1-\exp(-2\alpha)}\exp(-y){\,\rm d}y=1$, hence there is no stationary state of the form $G_{\alpha}$ with $\alpha>0$.

The case $\alpha<0$ is done similarly. ∎

Lyapunov functional

As it is classical we note the relative entropy

[TABLE]

and the Fisher information

[TABLE]

Note that $\mathbf{H}(c|G_{M})\geq M\log M$ for all $t>0$ by Jensen’s inequality. We introduce a Lyapunov functional for equation (24):

[TABLE]

where $\mathbf{J}$ is given by 26. Let $\Gamma_{c}$ be defined by

[TABLE]

with

[TABLE]

Lemma 5**.**

For the problem (24), if $M<1$, then the Lyapunov functional $\mathbf{L}$ is non-increasing:

[TABLE]

where the dissipation is

[TABLE]

Proof.

We compute the evolution of the entropy:

[TABLE]

where $\int_{-1}^{1}\partial_{t}c(t,x)(1-\log(M/\int_{-1}^{1}e^{-\frac{y^{2}}{2}}{\,\rm d}y)){\,\rm d}x=0$ by mass conservation. Hence,

[TABLE]

We can eliminate $c(t,-1)-c(t,1)$ from (34) in the two following steps:

[TABLE]

leading to

[TABLE]

Combining (34) – (35) – (36), we obtain

[TABLE]

and the proof of Lemma 5 is complete. ∎

A priori estimates

We now derive a priori bounds for solutions to (24) in the classical sense.

Proposition 6** (Main a priori estimate).**

Assume that $\int_{-1}^{1}c_{0}\log c_{0}{\,\rm d}x<+\infty$. Let $c$ be a classical solution to (24). If $M<1$, then the following global estimates hold true for all $T>0$:

[TABLE]

Proof.

The proof follows from Lemma 5. Indeed, integrating (32) in time, it yields that

[TABLE]

Since $\mathbf{H}(c|G_{M})\geq M\log M$, it remains to prove that $c(t,-1)-c(t,1)$ belongs to $L^{2}$ locally in time. We first derive the following trace-type inequality:

[TABLE]

Furthermore, recalling that

[TABLE]

and

[TABLE]

we deduce that

[TABLE]

which together with inequality (38) yields that

[TABLE]

The quantity $M+M^{-1}-2$ is positive since $M<1$. Hence, using (37) and (40) we can prove that $c(t,-1)-c(t,1)$ belongs to $L^{2}$ locally in time. Next, using again (40) togther with (39), we see that

[TABLE]

Hence, $\int_{-1}^{1}c(t,x)\left(\partial_{x}\log c(t,x)\right)^{2}{\,\rm d}x$ belongs to $L^{1}$ locally in time. ∎

Long-time behaviour

To prove convergence of $c(t,\cdot)$ towards $G_{M}$ we develop the following strategy. We use the previous a priori estimates which enable to pass to the limit after extraction of a converging subsequence. The main argument (apart from passing to the limit) consists in identifying the possible configurations $c_{\infty}$ for which the dissipation $\mathbf{D}$ vanishes. This occurs if and only if both positive terms in (33) are zero. Thanks to (26), this means that $\mathbf{J}_{\infty}=(1-M)\left(c_{\infty}(-1)-c_{\infty}(1)\right)$ on the one hand, and on the other hand,

[TABLE]

which yields that $c_{\infty}\equiv G_{M}$.

Rate of convergence

As it is classical when the equilibrium state is a gaussian function, the natural tool is a logarithmic Sobolev inequality established by Gross in [15] that we first recall. Although we are dealing here with a non linear problem this method will be fruitful.

Lemma 7** (Logarithmic Sobolev inequality).**

Let $\nu(x)dx=\exp(-V(x))dx$ be a measure with smooth density on $[-1,1]$. Assume that $V^{\prime\prime}(x)\geq 1$ then, for $u\geq 0$ satisfying $\int_{-1}^{1}u(x){\,\rm d}x=\int_{-1}^{1}\nu(x){\,\rm d}x$, we have

[TABLE]

First, recalling (32) and (33), we deduce that

[TABLE]

Our aim is to apply a logarithmic Sobolev inequality to $\mathbf{H}(c|\Gamma_{c})$. We first observe that

[TABLE]

and we decompose the relative entropy as follows

[TABLE]

Recalling the definition (30) of $\Gamma_{c}$ we deduce that

[TABLE]

thanks to the definition (31) of $A_{c}$. Therefore,

[TABLE]

Finally, applying a logarithmic Sobolev inequality (Lemma 7) to the measure $\Gamma_{c}(x){\,\rm d}x$ and using that $-(1-M)\left(c(t,-1)-c(t,1)\right)^{2}\leq 0$, we obtain

[TABLE]

thanks to Jensen’s inequality. Hence

[TABLE]

leading to

[TABLE]

A rate of convergence for the $L^{1}$ norm is obtained by using the Csiszár-Kullback inequality, [9], that we recall now.

Proposition 8** (Csiszár-Kullback inequality).**

For any non-negative functions $f,g\in L^{1}(-1,1)$ such that $\int_{-1}^{1}f(x){\,\rm d}x=\int_{1}^{1}g(x){\,\rm d}x=M$, we have that

[TABLE]

Hence, the following decay estimate holds:

[TABLE]

3.2 Critical case $M=1$

Similarly, we can prove the existence of a stationary solution.

Lemma 9**.**

For $M=1$, equation (24) admits a unique stationary solution given by $G_{1}$.

Proof.

The proof is similar to the one of lemma 4. ∎

In this part we cannot follow the strategy developped in Section 3.1 since we crucially used $M<1$. In the case $M=1$, from (34) it follows that

[TABLE]

and

[TABLE]

Then, together with the logarithmic Sobolev inequality $2\mathbf{H}(c|\Gamma_{c})\leq\mathbf{I}(c|\Gamma_{c})$, we deduce that

[TABLE]

Consequently, using again Jensen’s inequality, it follows that

[TABLE]

Consequently we deduce that $0\leq\mathbf{H}(c|G_{1})(t)\leq\mathbf{H}(c|G_{1})(0)$, hence

[TABLE]

A priori bound

To obtain a control on the dissipation of entropy, we follow the strategy developed in [8]. Consider the even function $\Lambda:\mathbb{R}\rightarrow\mathbb{R}_{+}$ such that $\Lambda(0)=0$, and $\Lambda^{\prime}(u)=(\log(u))_{+}^{1/2}$ for $u>0$. Then, it is non-increasing on $(-\infty,0)$, non-decreasing on $(0,+\infty)$, convex and superlinear, and for all $D$, there exists $A\in\mathbb{R}_{+}$ such that for any $u\in(-\infty,-A)\cup(A,+\infty)$, $\Lambda(u)^{2}\geq(1+D)C_{0}u^{2}$. Using again a trace-type inequality, we get

[TABLE]

Moreover, recall that we have

[TABLE]

with $\mathbf{K}(t)=\int_{-1}^{1}x^{2}c(t,x){\,\rm d}x\geq 0$. Combining (45) and (46), it leads to

[TABLE]

We can assume $A>\sqrt{2}$, so that $\alpha(t)^{2}>2$. In the case $|\alpha(t)|>A$, as $\mathbf{J}(t)=\mathbf{J}(0)e^{-t}$, we have $-\alpha(t)\mathbf{J}(t)\leq\alpha(t)^{2}|\mathbf{J}(0)|$, and

[TABLE]

Take $D=2+|\mathbf{J}(0)|$. Then,

[TABLE]

Now, on the set $E=\{t:|\alpha(t)|>A\}$, we have

[TABLE]

giving a $L^{2}$ control on $\alpha(t)$.

Finally, using (3.2) and (48), we prove that $\int_{0}^{t}\int_{-1}^{1}c(s,x)\left(\partial_{x}\log(c(s,x))\right)^{2}{\,\rm d}x{\,\rm d}t$ is bounded for all $t\in(0,T)$. Then, proposition 6 is still verified, and we can similarly prove the global existence of weak solutions in the critical case.

Long-time behaviour

The convergence of $c(t,\cdot)$ towards $G_{M}$ is proved as in the subcritical case.

Rate of convergence

By (45), we know that

[TABLE]

Hence, the Csiszár-Kullback inequality leads to the following estimation:

[TABLE]

3.3 The Super-critical case $M>1$

We prove now that depending on the mass, the problem (24) admits one or three stationary solutions. Then, we prove that if the initial condition is asymmetric enough, a blow-up occurs in finite time, and we are able to give a quantitative criterion for its appearance. The proof is based on the blow-up analysis for a modified problem, followed by the use of a concentration-comparison principle to deduce the blow-up in our case. Numerical simulations displayed in Figure 1 show the appearance of a blow-up.

3.3.1 Stationary solutions

We establish now the following lemma.

Lemma 10**.**

Denote $M_{0}=\frac{1}{2}\int_{-1}^{1}exp\left(-\frac{x^{2}-1}{2}\right){\,\rm d}x>1$. Then,

• $\bullet$

for $M\in(1,M_{0})$, equation (24) admits exactly three stationary states: the symmetric solution $G_{M}$, and two asymmetric solutions $G_{\pm\alpha}$, with

[TABLE]

for some $\alpha>0$ defined by the mass constraint $\int_{-1}^{1}G_{\alpha}(x){\,\rm d}x=M$.

• $\bullet$

for $M\geq M_{0}$, $G_{M}$ is the only stationary solution.

Proof.

We have seen in the proof of lemma 4 that the only possible stationary solutions are of the form $G_{M}$ and $G_{\alpha}$. It is straightforward that for all $M>1$, $G_{M}$ is solution and satisfies the mass constraint. For $G_{\alpha}$, let us denote $M_{\alpha}=\int_{-1}^{1}G_{\alpha}(x){\,\rm d}x$. We know that $M_{\alpha}>1$. It remains to characterize the set $I=\{M_{\alpha},\,\alpha>0\}$ of attainable mass values. As we were not able to characterize $I$ explicitly, we performed a numerical simulation of $M_{\alpha}$ as a function of $\alpha$ (see figure 3, a)). Notice that $M_{0}$ is defined such that $\lim_{\alpha\rightarrow 0}G_{\alpha}=G_{M_{0}}$, leading to the result. ∎

3.3.2 Blow-up for a modified equation

In this section we consider the following modified equation of (24):

[TABLE]

where $c(t,x)$ is defined for $t>0$ and $x\in(-1,1)$, together with zero-flux boundary conditions:

[TABLE]

and an initial condition: $c(t=0,x)=c_{0}(x)$. The total mass will also be denoted by $M$, i.e. $M=\int_{-1}^{1}c(t,x){\,\rm d}x$.

In this part we will prove that solutions to (49) - (50) blow-up in finite time when mass is super-critical $M>1$ and $c_{0}$ is decreasing and satisfies $c_{0}(-1)-c_{0}(1)>1$. To do so we prove that the first momentum shifted in $x=-1$ of $c$, $\tilde{\mathbf{J}}(t)=\int_{-1}^{1}(x+1)c(t,x){\,\rm d}x$ cannot remain positive for all time. This technique was first used by Nagai [35], then by many authors in various contexts. In a first step, in Lemma 11, we establish that the assumption that $c_{0}$ is a decreasing function such that $c_{0}(-1)-c_{0}(1)>1$ and that the $\tilde{\mathbf{J}}(0)$ is sufficiently small guarantee that $c(t,\textperiodcentered)$ is also decreasing and satisfies $c(t,-1)-c(t,1)>1$ for any existence time $t>0$. In a second step, in Proposition 12, we prove the blow-up character.

Lemma 11**.**

Assume that $M>1$. Assume in addition that the initial condition $c_{0}$ of (49) is decreasing, satisfies $c_{0}(-1)-c_{0}(1)>1$ and $\int_{-1}^{1}(x+1)c_{0}(x){\,\rm d}x<(M-1)/2$. Then any solution $c$ to (49), if it exists, is non-increasing and satisfies $c(t,-1)-c(t,1)>1$.

Proof.

Since we supposed that $c_{0}(-1)-c_{0}(1)>1$, it remains true at least until a time $t_{0}\in(0,T)$, where $T$ is the existence time. We choose the maximal $t_{0}$ possible. For all $t\in[0,t_{0}]$, $c(t,\cdot)$ is decreasing. In fact the derivative $u(t,x)=\partial_{x}c(t,x)$, satisfies the following parabolic type equation without any source term:

[TABLE]

The solution is initially non-positive, and also initially non-positive on the boundary due to (50) and the assumption $c_{0}(-1)-c_{0}(1)>1$. From classical strong maximum principle, we deduce that $c(t,\cdot)$ is decreasing for all $t\in[0,t_{0}]$.

Therefore, for all $t\in[0,t_{0}]$, $-\partial_{x}c(t,x)/\left(c(t,-1)-c(t,1)\right)$ is a probability density. From Jensen’s inequality, we deduce the following interpolation estimate:

[TABLE]

hence, using that $c(t,-1)>c(t,1)$ for any time $t\in[0,t_{0}]$, it follows that

[TABLE]

Consequently, for all $t$ in $[0,t_{0}]$, the first momentum shifted in $x=-1$, $\tilde{\mathbf{J}}(t)=\int_{-1}^{1}(x+1)c(t,x){\,\rm d}x$, which is non-negative and such that $2\tilde{\mathbf{J}}(t)\geq 4c(t,1)$, for all $t$ in $[0,t_{0}]$, thanks to (52), satisfies:

[TABLE]

as $M>1$. Using again that $2\tilde{\mathbf{J}}(t)\geq 4c(t,1)$, we deduce that

[TABLE]

Since $\tilde{\mathbf{J}}(0)<(M-1)/2$ and $\tilde{\mathbf{J}}(t)\geq 0$ we deduce that for all $t\in[0,t_{0}]$: $\frac{{\,\rm d}}{{\,\rm d}t}\tilde{\mathbf{J}}(t)\leq 0$, hence $\tilde{\mathbf{J}}(t)\leq\tilde{\mathbf{J}}(0)<\frac{M-1}{2}$. Consequently, using (52), it follows that for all $t$ in $[0,t_{0}]$,

[TABLE]

hence $t_{0}=T$ is the existence time of $c$. ∎

We can now prove a blow-up result for (49) - (50).

Proposition 12**.**

Assume $M>1$ and that the first moment shifted in 1 is initially small: $\tilde{\mathbf{J}}(0)<(M-1)/2$. Assume in addition that $c_{0}$ is decreasing and satisfies $c_{0}(-1)-c_{0}(1)>1$. Then the solution to (49) - (50) with initial data $c(0,x)=c_{0}(x)$ blows-up in finite time.

Proof.

From Lemma (11) it follows that

[TABLE]

Hence, since $\tilde{\mathbf{J}}(0)<(M-1)/2$, we deduce the inequality

[TABLE]

Now, $\tilde{\mathbf{J}}$ being non-increasing on $[0,t_{0}]$, we deduce that

[TABLE]

and with $M>1$, we deduce that the maximal time of existence $t_{0}=T^{*}$ of a solution is necessarily finite. Finally, following [19], it can be proved that

[TABLE]

showing that the solution becomes singular in $T^{*}$. Otherwise a truncation method enables to prove local existence by replacing $c$ with $(c-K)_{+}$ for $K$ sufficiently large. ∎

3.3.3 Blow-up in the Super-critical case

In this part, we prove Theorem 2. To do so, following [28], we state a so-called concentration-comparison principle on the equation obtained from (24) after space integration. This principle together with the use of a subsolution that blows-up allow proving the blow-up character of solutions to (24) above the critical mass.

Lemma 13**.**

Let $C,\bar{C},\underline{C}$ be non-decreasing (in space) functions in $C^{1}(0,T;C^{2}([-1,1]))$ satisfying

[TABLE]

Assume that $\underline{C}(0,\cdot)\leq C(0,\cdot)\leq\bar{C}(0,\cdot)$ and that $\partial_{x}\bar{C}(0,-1)-\partial_{x}\bar{C}(0,1)>\partial_{x}\underline{C}(0,-1)-\partial_{x}\underline{C}(0,1)$, then the following inequality holds true for all $0<t\leq T$:

[TABLE]

Proof.

From (57), we deduce that $\delta C=\bar{C}-\underline{C}$ satisfies the parabolic inequation

[TABLE]

with $\bar{C}(t,-1)=\underline{C}(t,-1)=0$. Since we supposed that $\partial_{x}\bar{C}(0,-1)-\partial_{x}\bar{C}(0,1)>\partial_{x}\underline{C}(0,-1)-\partial_{x}\underline{C}(0,1)$, it remains true at least until a time $T^{\prime}\in(0,T)$, and we choose the maximal $T^{\prime}$ possible. Since $\partial_{x}\underline{C}\geq 0$, on the time interval $[0,T^{\prime}]$, we have

[TABLE]

Hence, by strong maximum principle [13], $\bar{C}>\underline{C}$ on $(0,T^{\prime})\times(-1,1)$. Furthermore, by Hopf Lemma (see [13]), we also have

[TABLE]

As $T^{\prime}$ is maximal we immediately conclude that $T^{\prime}=T$. ∎

Proof.

Proof of Theorem 2. Let $\underline{c}_{0}$ be a decreasing function such that $\underline{c}_{0}(x)\leq c_{0}(x)$ for all $x\in[-1,1]$ and which satisfies $\int_{-1}^{1}\underline{c}_{0}(x){\,\rm d}x>1$. Assume in addition that $c_{0}(-1)-c_{0}(1)>\underline{c}_{0}(-1)-\underline{c}_{0}(1)>1$. Finally assume that $\int_{-1}^{1}(x+1)\underline{c}_{0}(x){\,\rm d}x<(M-1)/2$, where we recall that $M=\int_{-1}^{1}c(t,x){\,\rm d}x>1$. First according to Proposition 12, we know that the solution $\underline{c}$ to (49) - (50) with initial data $\underline{c}(0,x)=\underline{c}_{0}(x)$ blows-up in finite time.

We will now prove that the distribution function $\underline{C}(t,x)=\int_{-1}^{x}\underline{c}(t,y){\,\rm d}y$ satisfies

[TABLE]

Indeed integrating (49) - (50) in space, one obtains

[TABLE]

Hence,

[TABLE]

where the last inequality follows from the non-decreasing character of $\underline{C}(t,\cdot)$ together with $x\geq-1$. Hence, from Lemma 13, it follows that $\underline{C}(t,.)<\bar{C}(t,.)$ for all time below the existence time $T$, where $\bar{C}(t,.)$ is any function defined in (57). In particular, for $c$ solution to (24) with initial data $c(0,x)=c_{0}(x)$, we can set $\bar{C}(t,x)=\int_{-1}^{x}c(t,y){\,\rm d}y$, which yields the blow-up character of $c$. ∎

4 The model with dynamical exchange of markers at the boundary: prevention of blow-up and asymptotic behaviour

In Section 3.3, we proved that finite time blow-up occurs in the basic model (24) when mass is super-critical $M>1$. Such a behaviour is not realistic from a biological viewpoint. Here we modifiy the model model (24) by considering markers that are sticked to the boundary and thus create the attracting drift. More precisely, in this part we will study the stationary states associated with the following model:

[TABLE]

together with the flux condition at the boundary:

[TABLE]

The quantities $\mu_{\pm}$ represent the concentrations of markers which are sticked to the boundary and thus create the attracting drift $\mu_{-}(t)-\mu_{+}(t)$. The dynamics of $\mu_{\pm}$ is driven by simple attachment/detachment kinetics. The mass of molecular markers is shared between the free particles $c(t,x)$ and the particles on the boundary $\mu_{\pm}(t)$. The boundary condition (59) guarantees conservation of the total mass:

[TABLE]

From (60), we easily deduce that finite time blow-up cannot occur since $|\mu_{-}(t)-\mu_{+}(t)|$ is bounded by $M$. As a consequence, the problem under study is a conservative convection-diffusion equation with bounded advection term. We denote by $m(t)$ the mass of free particles:

[TABLE]

The conservation of mass reads

[TABLE]

Stationary solutions $(c,\mu_{-},\mu_{+})$ are characterized by:

[TABLE]

where we have denoted $\alpha=c(-1)-c(1)$.

We establish the following lemma.

Lemma 14**.**

Denote $M_{0}=1+\frac{1}{2}\int_{-1}^{1}exp\left(-\frac{x^{2}-1}{2}\right){\,\rm d}x>1$. Then,

• $\bullet$

for $M\leq M_{0}$, the problem (58) - (59) admits a unique symmetric stationary solution $G_{M}$ defined by:

[TABLE]

• $\bullet$

for $M>M_{0}$, there are two other asymmetric solutions $G_{\pm\alpha}$, with

[TABLE]

for some $\alpha>0$ defined by the mass constraint

[TABLE]

Proof.

A simple computation yields that any stationary solution is either of the form $G_{M}$ or of the form $G_{\alpha}$. It is straightforward that for all $M>0$, $G_{M}$ is solution and satisfies the mass constraint. For $G_{\alpha}$, we need to characterize the set $I=\{M_{\alpha},\,\alpha>0\}$. It can be proved that for all $\alpha>0$, $M_{\alpha}=M_{-\alpha}$, and $M_{\alpha}\sim_{+\infty}\alpha$, and that $\lim_{\alpha\rightarrow 0}M_{\alpha}=M_{0}$. As before, we performed a numerical simulation of $M_{\alpha}$ as a function of $\alpha$ (see figure 3, b)), leading to $I=(M_{0},+\infty)$, hence the result. ∎

Numerical simulations for this model are presented in figure 4.

5 Perspectives

In this work, we presented a model of 1D cell migration based on the diffusion and advection of a molecular marker inside the cell, that itself exerts a feedback on this dynamics. The marker’s spatial repartition is supposed to be characteristic of the polarisation state of the cell. The resulting equation, a non linear and non local Fokker-Planck equation, was shown to admit a dichotomy behaviour at equilibrium. Below and at the critical mass, polarisation is not possible. Between the critical mass and a limiting mass $M_{0}$, the system admits two polarised equilibria. Moreover, we have proved that for initial concentrations ”steep enough”, a blow up occurs in finite time, corresponding to a polarised state, but mainly showing that the model lacks regularity and that there is room for improvement. We next presented a more realistic model, where attachment-detachment kinetics of the marker at the membrane is taken into account. In this case, we were able to prove again a dichotomy behaviour for the existence of asymmetric steady states, and that the additional dynamics at the boundary is enough to prevent the blow up appearance.

A few questions remain today unresolved, and will be the object of a future work. The stability of either symmetric or asymmetric stationary states ought to be discussed, as prospective numerical simulations suggest non trivial results. Also, in the super critical case, it is of interest to search for a quantitative criterion to discriminate between the three stationary solutions. Finally, the controllability of the cell velocity can be studied.

This model was designed to investigate the initiation of polarisation during motion. In order to describe a whole displacement, a stochastic instability could be included in the marker’s dynamics. Then, the dichotomy’s result could be investigated at long time scales. The question of an external signal perturbing the molecular dynamics is of similar nature.

Moreover, a natural continuation consists in studying the corresponding 2D model, and the equivalent problem on a free boundary domain, in order to describe geometrical feedbacks of the cell on its motion.

*Acknowledgement: the authors are very grateful to V. Calvez, B. Maury, A. Mogilner and M. Piel for very helpful discussions. *