A congestion model for cell migration

This paper deals with a class of macroscopic models for cell migration in a saturated medium for two-species mixtures. Those species tend to achieve some motion according to a desired velocity, and congestion forces them to adapt their velocity. This adaptation is modelled by a correction velocity which is chosen minimal in a least-square sense. We are especially interested in two situations: a single active species moves in a passive matrix (cell migration) with a given desired velocity, and a closed-loop Keller-Segel type model, where the desired velocity is the gradient of a self-emitted chemoattractant. We propose a theoretical framework for the open-loop model (desired velocities are defined as gradients of given functions) based on a formulation in the form of a gradient flow in the Wasserstein space. We propose a numerical strategy to discretize the model, and illustrate its behaviour in the case of a prescribed velocity, and for the saturated Keller-Segel model.

1. Introduction, modelling aspects. We propose a model to handle interactions between organisms such as unicellular organisms e.g. bacteria or amoebia. The behavior of an entity is based on a will to move with its desired velocity, regardless of others, but the fulfillment of individual wills is made impossible because of congestion. We consider here the case of two organisms (see Section 4 for a straightforward extension to several populations). We describe the densities by measures ρ 1 and ρ 2 . We assume global saturation of the domain, i.e. ρ 1 + ρ 2 = 1. In order to take into account the individual tendencies together with congestion constraints, we consider that individuals have two velocities, namely a desired velocity that will be denoted by U 1 (resp. U 2 ) and a common correction velocity that will be denoted by w. Densities ρ 1 and ρ 2 satisfy: We consider the correction velocity which minimizes L 2 norm among all those velocity fields which ensure preservation of the saturation constraint ρ 1 + ρ 2 = 1, which can be expressed in a dual, Darcy-like, form: This basic model corresponds to a competition between two species which tend to achieve a prescribed motion, and realize a sort of compromise. We shall be especially interested in the following situations: (i) Competition between species which tend to minimize a given function: U i is defined as −∇D i , with ρ 1 U 1 + ρ 2 U 2 not feasible (i.e. leading to violation of the saturation constraint) in general. This setting, which extends the model proposed in [11] to a two-population situation with distinct tendencies, shall be the core of the theoretical analysis proposed in section 2. (ii) Desired velocity for species 2 is zero, and U 1 is prescribed as the gradient of a given function S (concentration of a chemoattractant). Note that this case, which corresponds to the migration of cells in a passive biological matrix, as encountered in the developpment of atherosclerosis, is a particular case of (i). (iii) Keller-Segel model with congestion (closed-loop version of the basic model): U 2 is again 0, and U 1 is defined as ∇S, where S is a chemotactic agent created by the species 1 itself. Assuming S diffuses in the mixture domain, we shall consider ∂ t S − ∆S = ρ 1 , or, assuming diffusion is instantaneous, −∆S = ρ 1 . The question of boundary conditions rises delicate issues. Setting equation in the whole space amounts to deal with infinite quantities of 1 and/or 2 (the sum of densities is equal to 1), which rules out standard tools for theoretical analysis and numerical simulation. One may consider that both populations occupy a bounded, moving domain Ω(t), on the boundary of which pressure is zero. It is then natural to consider that the boundary of Ω moves with a normal velocity which identifies to the normal velocity of the mixture (ρ 1 U 1 + ρ 2 U 2 + w) · n. This approach is surely relevant in some situations, yet we shall not consider it in this paper. Let us give an idea of one of the difficulties which are likely to occur: consider the monodimensional situation with initial condition and desired velocities: Translation of those intervals at constant speed 1/2 is obviously a solution to our problem (with piecewise affine pressure, 0 at both ends, 1/2 at the interface). Now consider the scenario where 1 runs away with velocity 1, and 2 stays where it is (both species are separated, and the domain Ω(t) is no longer convex). One can even imagine that half of 1 (initially located in 1 (1/2,1) ) takes off with velocity 1, while the other half of 1 trails 2 at speed 1/3. In this spirit, one can create an infinite number of scenarios which all seem to be, in a reasonable sense, solutions to our problem. Note that those problems are likely to occur in the case of expansion fields, as the one considered previously, because of the loss of monotonicity of the underlying evolution equation.
Another approach consists in considering a fixed domain Ω, delimited by rigid walls. This assumption is reasonable if one considers migration phenomena in living organisms. The boundary condition for Darcy equations is of Neumann type: noflux condition at walls leads to (ρ 1 U 1 + ρ 2 U 2 + w) · n = 0, i.e.
This condition ensures no-flux for the mixture, but not for both species independently; it may allow some out-or in-flux of either 1 or 2, which calls for further prescriptions on the boundary. The situations we shall consider in next sections suggest that, at least in the case U 2 = 0, the model tends to create pure zones (either 1 only, or no 1 at all) in the neighbourhood of boundaries. As a consequence, global no flux conditions turn into no flux conditions for both species. Yet, in 2 or 3 dimensions and for non zero desired velocities, one may have to consider more general situations. To overcome those difficulties, we shall consider in the theoretical part (Section 2) the fully periodic setting. In a different context, such saturated two-component mixtures have been studied in [13], with a motion driven by chemical potentials. In the biological context, a Keller-Segel model with logistic sensitivity (KSLS) has been proposed recently [6,4]. This KSLS model reads where the "desired" velocity U is of the chemotactic type: First of all, in one dimension (say Ω = (0, 1)) with no-flux boundary conditions at x = 1, our model reads with (no flux condition) ρU − ∂ x p = 0 at 1, so that ∂ x p = ρU in the whole interval, and therefore it identifies exactly with KSLS. Note that in case the velocity is given and constant, we recover an inviscid Burgers equation.
Both models are different for d ≥ 2. Yet, they present some common structure which can be described as follows: consider (situation (iii) above) the case of zero velocity for species 2, and velocity for 1 given as the gradient of a density S of chemoattractant created by 1 itself. We consider the biperiodic setting, and define ∆ −1 as the operator which maps a function g with zero mean value over the biperiodic domain onto the zero mean value solution to ∆u = g.
The pressure in our model can then be written p = ∆ −1 ∇ · (ρU), so that transport equation for species 1 becomes which is formally equivalent to KSLS model (3) with identity operator replaced by Helmoltz projection ∇∆ −1 ∇· (projection onto the space of irrotational fields). Nonlocality of the latter operator differentiate both models. Yet both PDE's systems are quite similar from the spectral point of view.
As for modelling aspects, they rely on different assumptions: KSLS model can be seen as standard KS model with a reduced velocity U = (1 − ρ)∇S, which can express for example the fact that when entities reach some critical density (set to 1 here), they disturb each other and are no longer able to sense the underlying gradient. In our approach, we consider that entities continue to exert some action which would lead them in the right direction if they were alone, but they are (mechanically) prevented from fulfilling their purpose, because of the presence of other entities (possibly with different "strategies"). The model is in some way less constrained as motion of saturated zones is possible, but also more constrained because the other species (which replaces empty space in KSLS) has to be swept away, at some price, by species 1.
2. Theoretical framework. We consider in this section the case of 2 species in the flat torus Ω = S 1 × S 1 , with desired velocities prescribed as follows where D 1 and D 2 are given smooth functions over Ω. The system corresponding to this situation writes Although concentrations for both species automatically admit a density which is in L ∞ (thanks to the congestion constraint), we shall keep the general setting of measures, which is usually adopted in optimal transport, and which allows for generalisations (e.g. relaxing the congestion constraint in some zones). Besides, we disregard normalization to obtain probability measures: we shall consider that ρ 1 ⊗ ρ 2 belongs to P(Ω× Ω), although the total mass is not 1. Similarly, we consider that ρ i ∈ P(Ω).
In what follows, regularity of functions defined over Ω accounts for periodic conditions. In particular H 1 (Ω) is defined as the set of biperiodic functions with H 1 regularity (closure for the H 1 norm of biperiodic regular functions).
Definition 2.1 (Weak solutions). We say that (ρ 1 , ρ 2 ) is a weak solution of (4) with initial condition (ρ 0 1 , ρ 0 2 ) if ρ 1 + ρ 2 = 1 for a.e. t, and if there exists p ∈ H 1 (Ω) such that for all ϕ 1 , ϕ 2 ∈ C ∞ c ([0, T [×Ω), for a.e. t ∈ (0, T ), and for all q ∈ H 1 (Ω), we have Let ρ = ρ 1 ⊗ρ 2 for a.e. t, U = (U 1 , U 2 ). We can rewrite the previous formulation as follows: The proof of this theorem relies on the notion of gradient flow in the Wasserstein space (see e.g. [1]). Theoretical analysis of a similar problem in the context of crowd motions was proposed in [11]. The proof proposed therein is quite general, and the same approach could be carried out in the present context. Yet, to avoid some technicalities and to give a clearer view of the underlying gradient flow structure, we propose here an alternative approach which is directly based on the main existence and characterization theorems in [1]. The rest of this section describes the main outlines of this approach, which is based on theories of optimal transport and gradient flows, which we will use to prove theorem 2.2. For more details, see the books of Villani [17], [18], and Ambrosio et al. [1], [2]. See also the recent application of this framework [7] to recover entropic solutions to the Burgers' equation (which identifies in some way to the proposed model for d = 1, as pointed out in the introduction).
We define a discrete scheme for (4): . Let τ > 0 be given, and ρ 0 = ρ 0 1 ⊗ ρ 0 2 ∈ P(Ω × Ω) an initial density. We define ρ 1 , . . . , ρ n (we drop the dependence upon τ to alleviate notations) recursively according to where J is given by and K is the set of admissible densities The notation I K stands for the indicatrix function of K, i.e.
This is a well-known scheme in gradient flow theory (see [5], [8], [3], [1], [2]), which has been widely used to prove existence theorems. Taking a minimizing sequence of the right-hand side of (5), we can verify that every of these minimizing problems admit at least a solution.
Let us underline that here the JKO scheme is only a tool to prove the existence of a solution to our problem, and has not been used for numerical purposes. As a matter of fact, the numerical resolution at each time step of the minimisation problem leads to many difficulties, and we have chosen a totally different approach for the numerical tests in section 3.
Thanks to the previous proposition, the JKO scheme can be rewritten as follows Proposition 2 (Convergence of the JKO scheme). Let ρ τ be the piecewise constant interpolation (in time) of the sequence (ρ n ) n defined by the JKO scheme (5) for the time step τ . Under the assumptions of theorem 2.2, there exists a subsequence of (ρ τ ) τ which converges weakly to a weak solution of the transport equation Proof. We have to verify that the functionnal Φ := J + I K satisfies the assumptions of the theory developed by Ambrosio et al. in [1]. Φ is clearly proper and lower semicontinuous, and its sublevels [Φ(ρ) ≤ c] are compact. Moreover, we can prove that Φ is regular according to definition 10.1.4 in [1], i.e. for all sequence ρ n = (ρ n 1 ⊗ ρ n 2 ) ∈ K, for all u n ∈ ∂Φ(ρ n ) such that (ρ n ) narrowly converges towards ρ in P(Ω × Ω), and sup n ||u n || L 2 (µn) < +∞, u n ⇀ u weakly then u ∈ ∂Φ(ρ). All assumptions of Prop. 2.2.3, Th. 2.3.1, and Th. 11.1.3 of [1] are satisfied, therefore there exists a subsequence of (ρ τ ) τ which converges to a weak solution of ∂ t ρ + ∇ · (ρu) = 0, with u ∈ −∂Φ(ρ) for a.e. t.
Proposition 3. The limit of (ρ τ ) τ in proposition 2 is actually a weak solution of equation (4).
We now use the fact that −u is a strong subdifferential, i.e. that for all transport map r, we have Let us underline that if r # ρ ∈ K, the previous inequality does not give any information, as Φ(r # ρ) = +∞.
We define the set of admissible velocities as follows We just proved that u ∈ C ρ . Let v ∈ C ρ , ε = 0, and r ε = id + εv. In most cases, we have r ε # ρ ∈ K. However, it is possible to find a transport t ε such that Using estimates between (t ε • r ε ) # ρ and r ε # ρ, we get Letting ε go to 0 from positive and negative values, we obtain Ω×Ω −∇D − u, v dρ = 0, which means that −∇D − u ∈ C ⊥ ρ . As u ∈ C ρ , u is indeed the projection of −∇D = (U 1 , U 2 ) onto C ρ , i.e. This minimizing problem can be written as a saddle-point problem for which we can prove existence and uniqueness of a solution (u, p) which satisfies Moreover, since ∇ · (ρ 1 u 1 + ρ 2 u 2 ) = 0, the equation on p is given by If we replace the expression of u in equations (7), we finally get i.e. ρ is a solution of (4). Under reasonable assumptions on D 1 and D 2 , it is to be expected that the JKO process leads to a unique solution (see [11] for remarks regarding uniqueness in a similar setting). Yet, as we pointed out in the introduction, the system is under some conditions equivalent to the standard inviscid Burgers' equation, which rules out uniqueness in general. It is natural to wonder whether the JKO scheme selects a particular solution, namely the entropic one. This delicate question is still widely open, but a recent work on Burgers equation ( [7]) suggests a positive answer.
3. Numerical methods and results. We consider here the transport model with congestion, with one active species only (U 1 = U, U 2 = 0) expressed in terms of the active species only: −∆p = −∇ · (ρU ) .
3.1. Discretization and maximum principle. First, a time discretization of the system (8)-(10) is given. We set t n = n δt, and Equation (8) is approached by using a backwards Euler finite difference method. Semi-discretized version of System (8)-(10) hence reads: For the sake of simplicity, we only give the 1D spatial discretization of (11)-(13), the extension to 2D or 3D on cartesian grids being straightforward. We introduce: Since the equations of the model are written in a conservative form, the natural framework to be used for the spatial discretization of (11)-(13) is the finite volume framework. We hence introduce the control volume defined by: We denote by (ρ n i ) i and (p n i ) i the piecewise constant approximations of densities and pressures at time t n (i.e. ρ n i stands for the value of ρ at cell center x i , at time t n ). As for flux variables, we define approximate values at interfaces: w n i (resp. U n i ) stands for w (resp. U) at interface x i− 1 2 , By integrating (11) over C i and by using the divergence theorem, we obtain We approximate the integral terms in the following way: An upwind discrete flux is then introduced to approximate the remaining flux term in (14): where the numerical flux A up is defined by We finally obtain ; It is important to notice that fluxes ρU and ρw were treated separately in (15). As we will see, this plays an essential role in preserving the maximum principle on ρ for the numerical solution. This advection scheme is stable under the Courant-Friedrichs-Lewy condition : Eq. (13) is discretized in space with Finally, by using an Euler finite difference scheme in (12), correction velocity w is approximated with For the sake of simplicity, only periodic boundary conditions are considered in the following proof: ρ n 1 = ρ n Nx , p n 1 = p n Nx . Note that, in practical, the proposition remains true if the no-flux boundary condition (mentioned in Section 1) is taken over p. In this case, appropriate boundary conditions have to be taken upon ρ in the upwind scheme: First the positivity of the numerical scheme is given by the well known positivity of the upwind scheme under the C.F.L. condition (17). For the upper bound ρ n i ≤ 1, the proof relies on the positivity of µ n i := 1 − ρ n i , that is given by a similar numerical scheme. Indeed, equation (16) implies By replacing the values of w n i with the discrete gradient of p as shown in (19), one gets We recognize in the first terms of the right hand side of (20) the exact discrete expression of (18). Hence the numerical scheme satisfied by µ is which is an upwind discretization of the advection problem: ∂ t µ + ∇ · (µ w) = 0. Thanks to the positivity of the upwind scheme under the C.F.L. condition (17) we deduce the discrete maximum principle for ρ.

3.2.
Numerical results. In this section we present several numerical simulations performed for the migration model (8)-(10) with the numerical scheme (16)-(19) introduced above. First, in order to validate the numerical method, two 1D test cases for which the exact solution is known are presented. Then 2D simulations are shown for two types of situations: first the case where the desired velocity U is known, and then the case where U is given as a function of ρ.

1D simulations.
In this section the domain is the interval (0, 1), the desired velocity is U = 1, and we consider the following boundary conditions for p: Two different initial conditions are considered: Fig. 1a) , and Fig. 2a).  expected to diffuse fronts: it is due to the behavior of the model we consider, which tends to sweep ρ toward the discontinuity, inducing permanent correction of the numerical diffusion.

2D simulations for a constant desired velocity U.
In this section we perform some numerical simulations for the migration model (8)- (10), in the case where the desired velocity U is given and constant in time. In all the cases shown here, the domain Ω is bounded and surrounded by walls, so that p verifies Neumann boundary conditions: The first situation studied here is a direct 2D extension of the 1D test case described previously, see Fig. 2. The domain is the unit square Ω = (0 , 1) × (0 , 1). We also prescribe the following desired velocity: U = 1. Figure 3 shows the evolution of ρ from the initial condition (Fig. 3a) to the steady state (Fig. 3f) given by a numerical simulation performed in the situation described above. The same way as in the 1D situation (Fig. 2b), the evolution of ρ begins with the spreading of a mixing zone along the x axis (Figs. 3b-c). Along with the evolution of the mixing zone, the  solution also tends to spread along the y axis. Hence, a real difference is shown here between our model and Burgers-type models. When the mixing zone reaches the wall, a congested zone instantly appears and develops backwards (Fig. 3d-e). The steady state is reached when all the right side of the domain is saturated with the active species (Fig. 3f).
The second situation studied in this section involves a more complex geometry for the domain, which is defined by which corresponds to two rectangular rooms joined by a thin corridor (see Fig. 4). The desired velocity is set to U = −∇D, where D(x, y) represents the geodesic distance from (x, y) to the right wall. In practice D is computed thanks to the toolbox provided by [15], which is based on the Fast-Marching method (see [10]). Figure 5 represents the simulated evolution of the concentration ρ with the setting described above. The initial condition is shown in figure 5a. As we can see a congested zone rapidly forms near the entrance of the corridor (figure 5b), which size decreases in time as the matter runs through it (figures 5c-e). We notice that a quasi-homogeneous flow regime tends to be established in the corridor (Fig. 5d-e),  Figure 3. Evolution of the concentration ρ in the case of a constant desired velocity U = 1, with the "wall" boundary condition (24). Mesh resolution: N x × N y = 300 × 300.
where a constant concentration ρ ∼ 0.5 is observed. It illustrates conservation of species 2 (in white in the figures): as the domain is bounded by walls, flow of 1 to the right has to be balanced by an opposite flow of 2 to the left. As we define the correction velocity in the least-square sense, actual speeds for 1 and 2 are close, so that the mixture is necessarily balanced for global mass conservation reasons. After some time, as observed in the previous test case, all the active species is concentrated on the right side of the domain (Fig. 5f).
The numerical treatment we give to equations (26)-(27) is similar to the numerical treatment of equations (9)-(10) described previsouly: For all the following numerical simulations the initial data is set randomly according to Bernoulli's law where q ∈ (0, 1) conditions the initial mass of active cells. The domain is the unit square, the boundary conditions set here are periodic boundary conditions for ρ, p and S.  Figure 5. Evolution of the concentration ρ in the case of a constant desired velocity U = ∇D (see Fig. 4), with the "wall" boundary condition (24). Mesh resolution: N x × N y = 300 × 300. and 7a). In each case we see aggregation of bigger and bigger structures as time goes (Figs. 6b-e and 7b-e). Two different steady states are shown, both involving a steady, fully congested zone. In the case q = 0.1 the steady state congested zone is a disc, whereas in the case q = 0.5 a band has formed because of the periodic boundary conditions we chose to consider. of the model as a gradient flow in a product space of densities, and proposed a discretization strategy which enjoys reasonable stability and accuracy properties. In terms of modelling, as mentioned in the introduction, any number N of species can be handled, as far as equations are considered (possible issues concerning boundary conditions are disregarded here). Saturation writes simply ρ 1 + · · · + ρ N = 1, we have an advection equation for each species ∂ t ρ i + ∇ · (ρ i (U i + w)) = 0 , i = 1, . . . , N, and the common correction velocity verifies It could be of particular importance to include also proliferation phenomena. Denoting by β i the local rate of creation (possibly depending explicity upon ρ i or other densities), conservation equations become ∂ t ρ i + ∇ · (ρ i (U i + w)) = β i , i = 1, . . . , N, and the constraint on w accounts for mass creation: Note that, in this situation, global non conservation rules out the use of no-flux conditions (or periodic setting). Considering a bounded domain Ω one can consider its boundary as a free outlet (interaction pressure is set at 0), so that a non conservative flux through boundary can compensate the unbalance of mass in the domain. In case of mass creation (β i > 0), one can expect outflow through the boundary, so that no additional condition is needed for the advection equations. In case some individual velocities U i +w may point inward the domain, the system has to be complemented with appropriate conditions (e.g. prescribed value of ingoing densities).
As for theoretical issues, non conservation rules out the standard framework of gradient flow in the Wasserstein space, which is dedicated to measures with constant total mass. Yet, as suggested in [11], generalizations of the JKO scheme, in the spirit of prediction-correction algorithms (or catching-up algorithms for sweeping processes, see for example [12]) might be expected to provide well-posedness results.
One may also wonder whether it could be possible to recover evolution equations for a single species, as in the case of crowd motion models ( [11]). In this situation, species 2 is replaced by empty space, which can be moved at no cost. Our proposed model could be seen as an attempt to replace a unilateral constraint ρ 1 ≤ 1 (which is quite delicate to handle, see again [11]), by an equality ρ 1 + ρ 2 = 1 with ρ 2 ≥ 0. It can be done formally by lowering the importance of species 2, more precisely by defining the correction velocity as the one which minimizes a weighted L 2 norm: where the minimum is taken among all those fields which satisfy ∇·w = −∇·(ρ 1 U 1 ) (in the case U 2 = 0), and k ε (ρ 1 , ρ 2 ) = ρ 1 + ερ 2 . One recovers formally an evolution equation for ρ 1 with unilateral constraint ρ 1 ≤ 1, yet with a significant difference: if one considers a saturated zone of 2 surrounded by a saturated zone of 1, the asymptotic (ε → 0) model sees the inclusion as globally incompressible, which is not the case if one imposes simply ρ 1 ≤ 1.