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

In this work, we present a numerical scheme for the approximate solutions of a 2D crawling cell migration problem. The model, defined on a non-deformable discoidal domain, consists in a Darcy fluid problem coupled with a Poisson problem and a reaction-advection-diffusion problem. Moreover, the advection velocity depends on boundary values, making the problem nonlinear and non local. \parFor a discoidal domain, numerical solutions can be obtained using the finite volume method on the polar formulation of the model. Simulations show that different migration behaviours can be captured.


Introduction
Cell migration ensures fundamental functions in the body (embryogenesis, immune system), but is also involved in the development of pathologies such as tumor metastasis, arising large research efforts (Bravo-Cordero et al., 2012).However, the responsible intracellular mechanisms involve multiscale interaction in time and space, so that modelling cell migration is challenging and produces interesting problems to study.
For 2D cells crawling on a surface, the motion is friction-based and results from the activity of the so-called actin cytoskeleton, that is a dynamics mesh of actin filaments.They are polar: they polymerize at one end and depolymerize at the other end, under the molecular regulation of signalling loops.Overall, the actin mesh can be approximated by an active fluid in our setting (Kruse et al., 2005;Joanny and Prost, 2009).
The mesh grows preferentially at the cell membrane, and shrinks inside the cell body, generating inward actin flows from the membrane.Its mechanical connection to the adhesive substrate generates the friction forces responsible for the cell's displacement.
Modelling this process is a difficult task, because of the large time and space scales, and also because of the large number of effectors of the dynamics.Following key physical ideas of Blanch-Mercader and Casademunt (2013); Maiuri et al. (2015), we proposed in Etchegaray (2016); Etchegaray et al. (2017b,a) a minimal multiscale model for 2D crawling migration, that we recall now.
The cell domain is a non deformable disc Ω, and the problem is formulated in the cell's frame of reference (the domain does not move).First, the actin cytoskeleton is approximated by a Darcy fluid (Blanch-Mercader and Casademunt, 2013), and its dynamics in a crawling situation is modelled by a Poisson problem on the fluid pressure (Etchegaray et al., 2017b).More precisely, for u : R + × Ω → R 2 the fluid velocity, p : R + × Ω → R its pressure, we have for k d the actin depolymerization rate, and k p : R + × Ω → R the polymerization function rate at the boundary.The fluid velocity writes with for γ ∈ R + the domain velocity for n x the unit normal vector to ∂Ω at point x ∈ ∂Ω 0 .These equations show that the actin polymerization at the boundary and depolymerization inside the domain may drive a pressure gradient leading to the cell motion.Note that the dynamics therefore arises from the activity at the boundary, and that the cell velocity is nonlocal.Now, the interaction with the molecular scale consists in studying the dynamics of a molecular inhibitor to polymerization.The molecules diffuse freely inside the cell in a inactive form.They may bind actin filaments and be carryied by their flow.Moreover, if activated at the cell membrane, they become able to locally inhibit actin polymerization.Write c : R + × Ω → R the concentration in an inactive form, and µ : R + × ∂Ω → R the activated concentration at the boundary.Then, → the inactive molecules follow an advection-diffusion dynamics with advection velocity u, → there is an exchange dynamics on ∂Ω between active and inactive forms, → k p is a decreasing function of µ.
For D the diffusion coefficient, and k on/off the activation/desactivation rates, the corresponding problem writes Note that the boundary condition ensures mass conservation: Some remarks can be made to highlight the difficulties in the analysis of the model.The fluid velocity rewrites We notice that this expression depends on the concentration in activated molecules µ from the pressure boundary term k p , so that the reaction-advection-diffusion problem is nonlinear.Moreover, the integral term makes it also non-local.The corresponding 1D model in a special case without activation at the boundary and for a linear k p writes as a nonlinear non-local advection-diffusion problem, and has been studied in Etchegaray et al. (2017b).The analysis shows the existence of different asymptotic solutions, describing both motile and non motile behaviours.Moreover, for a subcritical mass of molecules, the global weak existence of solutions is established, along with their convergence to a non motile gaussian profile at explicit rate.Finally, under conditions on the initial condition, it is shown that solutions blow up in finite time.
Since the model carries non trivial behaviours, computing numerical solutions is of particular interest.In the following, we develop a finite volume method for the polar formulation of the problem (1)-( 2)-( 4), since the domain Ω is a disc.

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

Poisson problem on p
For the polymerization function k p , let us choose a simple form, that is with δ > 0 and x ∈ ∂Ω.
The pressure p is solution of the following problem: where the pressure condition on C(0, R min ) is arbitrary, and fix the pressure values.Let us consider these equations in polar coordinates with 1 r p(t, r, θ

Discretization
Let t n = n ∆t be the time discretization, and Let cn (j,k) (resp.μn k ) be the approximated value of the exact solution c(t n , r j , θ k ) (resp.μ(t n , θ k )), and pn (j,k) be the approximated value of the exact solution p(t n , r j , θ k ): cn Moreover, we write u n and v n the corresponding discretized velocity functions at time t n .
The resolution is made as follows: for n ≥ 0, knowing (c n , μn ) allows to compute pn , then (u n , v n ).Finally, (c n+1 , μn+1 ) is computed using u n , and so on.

Problem on p
The pressure problem is time-dependent because of the Dirichlet boundary condition.Therefore, it is solved explicitly in time.Write F for the numerical flux.Then, we have the following scheme for equation (15 The finite volume numerical fluxes are defined by The Dirichlet boundary conditions ( 16)-( 17) are imposed using ghost values p(0,k) and pn (Nr+1,k) .For k ∈ {1, ..., N θ }, r Nr + , the corresponding flux writes Therefore, the term r Nr + will be included in the right hand side of the matricial problem.
(22) or equivalently Let us now write the corresponding matricial problem.We define the column vector P by P(k + (j − 1)N θ ) = p(j,k) with (j, k) ∈ {1, ..., N r } × {1, ..., N θ }: For ∆r = ∆θ the stiffness matrix A p is defined by where the second matrix accounts for the angular diffusion, with A ∈ M N θ (R) the classical diffusion matrix with periodic flux boundary conditions: The right-hand side in (15) and the flux boundary condition ( 16) on C(0, R max ) imposes this right hand side column vector of length N r N θ : .
We use a standard numerical method to invert the symmetric positive definite matrix 1 ∆r 2 A p and then resolve at each time step Equations for u and v The velocities u and v depend on the pressure p, which is obtained through a time explicit scheme.Therefore, they also depend explicitly in time on the concentration.The equation on v in polar coordinates writes We compute numerically the velocity in cartesian coordinates v n cart := (v n x , v n y ) T : Then, a polar change of coordinates leads to v n := (v n r , v n θ ) T .For the fluid velocity, we have that rewrites since u r = u(t, x) • e r and u θ = ru(t, x) • e θ .We define at time
We define now the numerical fluxes.The diffusion part is implicit, so that no CFL condition is needed (Allaire, 2005), while the advection is explicit due to the nonlinearity in the expression of v.
We write the corresponding scheme and group the implicit (resp.explicit) terms on the left-hand-side (resp.right-hand-side).Note also that each equation is written for k Now, for j ∈ 2, ..., N r − 1, we obtain 2 ) , cn (j,k) , cn (j,k+1) .
(34) The membrane activation equation writes in polar coordinates At each time step, the implicit discretization of equation ( 35) for k ∈ {1, ..., N θ } writes For simplicity, we treat both the free and activated concentrations in the same linear problem of unknown We have at each time step where A is the diffusion matrix, and B n the advection matrix.

The diffusion matrix writes
where the second matrix accounts for the angular diffusion, with A ∈ M N θ (R) the classical diffusion matrix with periodic flux boundary conditions: Now, for the advection term A up defined by equation ( 31), we write (u) + = max(u, 0) and (u) − = min(u, 0) so that A up (u, cj,k , cj,k+1 ) = (u) + c j,k + (u) − c j,k+1 .Therefore, we introduce the following diagonal matrices for j ∈ {1, ..., N r }, U + .
Thus we can define the advection matrix: where the discrete advection matrix B n ∈ M N θ (R) with periodic flux condition on the boundary is defined as in Allaire ( 2005) .

Results
The discretization scheme was implemented using MATLAB.We performed some numerical simulations to test the scheme's numerical convergence.Since the problem has a boundary nonlinearity, comparing a numerical solution to an exact one is out of reach in general.We consider the case of a polarised initial condition that writes c(0, r, θ) = cos(θ − π) + 1, μ(0, θ) = 0.5c(0, R, θ).We fix some parameters:

Illustrative example
By changing the value of the activation rate k on , we observe qualitatively different stationary solutions.The figure 1 Note that for k on = 0, it is clear that the stationary concentration in activated molecules is 0 everywhere on the boundary, so that v = 0 and the nonlinearity disappears.In the following, we rather focus on cases where the stationary state can be asymmetric.

Non-zero stationary polarisation
We also performed the same simulation for k on = 0.3 and R = 1, with an angular discretization step ∆θ = 2π/160 3.9 * 10 −2 , and varying values for ∆r and ∆t.The system is considered at numerical steady state when the concentration µ in activated molecules has stayed unchanged for 1 numerical hour.Then, we obtain the time to attain the steady state, as well as the stationary polarisation of the system, quantified up to a constant by the cell velocity vector.
We consider the following parameter values: The figure 2 shows that the time to attain the stationary state is a consequence of the polarisation state of the stationary solution.Indeed, considering that the initial condition is polarised, and the simulation shows the system's depolarisation towards a less polarised steady state, then the more polarised it is, the sooner it is attained.
We also clearly notice that the smaller ∆t gets, the more polarised the stationary state gets.This trend can be distinguished independently of ∆r.The radial step has a smaller effect, but more precise grids are correlated with lower polarised solutions.
Overall, the simulations show that a fine time step is fundamental to catch the relevant stationary type of behaviour.The polarisation module being decreasing with the time step, we can infer that the numerical solution approach a polarised state.However, this also shows that the previous numerical simulations are better understood in a qualitative sense rather than for quantitative purposes.

Low stationary polarisation
Finally, we perform the same numerical test with k on = 0.1 to check how the discretization steps may generate an error between a polarised and a non polarised stationary solution.We obtain the values depicted in figure 3.
We can see that the stationary polarisation levels are very low compared to the case where k on = 0.3.The same trend appears for the effect of ∆t, while the effect of ∆r is less visible since the nonlinear term is small.
To determine if these polarisation levels could be comparable to a true symmetric stationary state, we performed the same simulation for k on = 0 and for the most precise time-space grid.The stationary polarisation module was approximately equal to 1.63 * 10 −4 , so that no clear distinction can be made between these cases.

Conclusion
In this work, we have presented the finite volume discretization of a multiscale model for 2D cell crawling migration consisting in a Darcy fluid dynamics coupled to a Poisson problem and to a nonlinear and non-local reaction-advection-diffusion problem for the concentration in a molecular specie.The simulation of numerical solutions of this type of problems is very useful since the mathematical analysis is necessarily limited, whereas the model show varied behaviours.
The discretization method showed good qualitative numerical result.In particular, the molecular mass is preserved numerically.However, in critical cases, the scheme seems not able to make the distinction between polarised and unpolarised states.A natural continuation will consist in taking into account the interaction with the environment (chemical signal, mechanical obstacle) as a bias for motion.Finally, further studies should be based on an implicit treatment of the nonlinearity Cancès et al. (2017b,a).
shows the time evolution in the molecular concentration in the cell body in two characteristic cases.The left figure shows a typically non motile profile, while the right figure displays a polarisation situation.

Figure 2 :
Figure 2: Left: Time to attain the stationary state and Right: stationary polarisation module for varying ∆r and ∆t, and k on = 0.3.

Figure 3 :
Figure 3: Left: Time to attain the stationary state and Right: stationary polarisation module for varying ∆r and ∆t and k on = 0.1.