MICROTUBULES (MT) A KEY TARGET IN ONCOLOGY: MATHEMATICAL MODELING OF ANTI-MT AGENTS ON CELL MIGRATION

. Microtubules (MTs) are protein ﬁlaments found in all eukaryotic cells which are crucial for many cellular processes including cell movement, cell diﬀerentiation, and cell division, making them a key target for anti-cancer treatment. In particular, it has been shown that at low dose, MT targeted agents (MTAs) may induce an anti-migratory eﬀect on cancer and endothelial cells, leading to new prospects in cancer therapy. In that context, we propose to better understand the role of MT dynamics and thus of MTAs on cell migration using a mathematical cell centered model of cell migration taking into account the action of microtubules in the process. The model use a ﬂuid based approach that describes, through level-set techniques, the deformation of the membrane during cell migration. The ﬂuid part of the model is mainly composed of Stokes equations and the biochemical state of the cell is described using Reaction-Diﬀusion equations. Microtubules act on the biochemical state by activating or inactivating proteins of the Rho-GTPases family. The numerical simulation of the model is performed using Discrete Duality Finite Volume techniques. We describe the diﬀerent schemes used for the simulation, focusing on the adaptation of preexisting methods to our particular case. Numerical simulation are performed, showing a realistic behavior of the simulated cells in term of shape, speed and microtubules dynamics. Diﬀerent strategies for a depolymerizing MTA (Vincristin) mechanisms are investigated and show the robutness of our model.


Introduction
Cell migration is a fundamental biological process involved in organogenesis, tissue repair response and vascularization.In parallel, some pathological processes during cancer development, such as tumor progression and metastasis are also cell motility dependent [11].Therefore, the study and understanding of cell migration represent a real challenge in the development of new therapeutic approaches.
The cell migration is defined as a complex and dynamic mechanism initiated by the protrusion and adhesion of the cell membrane involving cytoskeleton and signaling networks regulated by Rho-family GTPases.This cellular event is polarized along a front-rear axis.At the front of the migrating cell, protrusive structures named, lamellipodia are driven by spatially and temporally-regulated actin polymerization and are generally stabilized by adhering to the extracellular matrix [4,41].On the other hand, at the cell rear, actin filaments associated with myosin exert a strength tension that induces a disassembly of adhesion sites allowing the cell tail detachment.Rho-family GTPases are the main regulators of these events that lead to an effective cell migration.Among Rho-GTPases, in this study we focus on Rac and Rho that mainly localize in front of the cell or the rear, respectively, and antagonize mutually.Rac and Rho are pivots of the polarity in migrating cells.Active Rac mediates actin polymerization in protrusion and then membrane extension, while active Rho increases actomyosin contractility and induces focal adhesion detachment [20].Rho regulates actomyosin formation through mammalian homolog of Diaphanous (mDia) and Rho-associated coiled-coil-containing kinase (ROCK) [37].The mDia protein promotes actin nucleation and polymerization [15] while ROCK activates myosin light chain and cooperates with mDia to form actomyosin bundles and generate contractility [37].In this paper, Rho is assumed to activate ROCK at the rear of the cell to regulate actomyosin contractility.Our model does not take into account the Rho effect on mDia at the front of the cell, which has a similar effect than Rac [37].Significantly, microtubules (MTs), another cystokeleton component, directly regulate Rho-GTPase activity and indirectly act on actin filaments and then cell movements [19,51].Especially, MT polymerization activates Rac to form protrusion, and depolymerizing MT activates Rho for adhesion disassembly and rear retraction.MTs are polarized and dynamic structures characterized by their dynamic instability.MTs undergo continuous cycles of polymerization termed growth with add of GTP-tubulin and depolymerization termed shrinkage.The alternance between MT growth and skrinkeage is mainly regulated by GTP-tubulin hydrolysis and GTP-tubulin concentration (see [5,30,50]).
The interest of mathematical modeling for cell migration really began in the 1980's with the works of Abercrombie [1].Since then, different models were proposed to study different mechanisms of cell migration.The field is currently very active, with the recent works of Berlyand and Mizuhara [36], on moving cells on patterned substrate, Meunier and Etchegaray [18] on persistence of the migration or Poignard [23] for models of invadopodia.
To our knowledge, the action of MT on cell migration has never been studied from a mathematical point of view.In this study, we design a first attempt of mathematical modeling to describe MTs impact on cell migration.We start from a polarization model proposed by Edelstein-Keshet and collaborators [47] that describes the effect and the distribution of Rac inside the cell during migration.We extend this model, introducing MT as additional material inside the cell, responsible of the activation and inactivation of the two Rho-GTPases Rac and Rho.The complete system consists in a Stokes equation coupled with transport-reaction-diffusion equations and differential equations posed in a moving domain.We adopt a level set approach to represent the moving domain.One of the difficulties encountered in the discretization of such problem is linked to the influence of the artificial boundary conditions for the Stokes problem on the solution.To reduce this influence, we classically have to work on a domain large in comparison with the cell size, leading to huge computational cost.To reduce this cost, we propose to work on grids that are locally refined in the neighborhood of the cell.To discretize the system on such a nonconformal grid, we adopt a Discrete Duality Finite Volume (DDFV for short) approach, known to be efficient in a large class of elliptic, parabolic problems [3,9,17,33].The DDFV discretization of our Stokes problem is an adaptation of [33].We have developed in [45,46] a DDFV-WENO scheme for the transport of the level set.We show in this paper how the DDFV setting enables us to take into account the possible large deformation of the cell in the discretization of the reaction-transport-diffusion equations in the moving domain.The numerical model that we propose, is then a robust tool to study the influence of microtubules (MTs) and the microtubule targetting agents (MTAs) on cell migration.We will see how the model is able to reproduce a good behaviour in terms of shape of the cell, in terms of migration speed, in terms of microtubule dynamic instability.
The outline of the paper is the following.In Section 2, we present the model used to describe the impact of microtubules on cell migration.In Section 3, we detail the DDFV discretization of the whole system and insist on the DDFV discretization of reaction-transport-diffusion problems inside the moving cell.Then, in Section 4 we describe the calibration of the model and the influence of parameters linked to a MT depolymerizing drug effect.We exhibit the good behaviour of the model.

Mathematical model of cell migration
The model that we have developed to describe the action of microtubule dynamics instability on cell migration is mostly based on the ideas proposed in [47].

Global description of the model
In the model we consider a two-dimensional (2D) representation of the cell from the top-down point of view.This representation is very suitable to describe in vitro situations where the cell adopts a very flat profile.It is commonly accepted in the literature to consider the nucleus as a passive material and to neglect its presence in the model [47].This assumption is justified by the fact that cell fragments devoid of nuclei can still migrate [48].
We adopt a fluid based model that considers the cell and its environment as a viscous fluid.Then the cell membrane is seen as an elastic material moving through this fluid.In order to describe the interaction of the cell with the extracellular matrix at the membrane, we use the Immersed Boundary Method.The main idea of this method is to model the different interactions at the membrane as forces exerted on the cell membrane.
The forces responsible for the displacement of the membrane are mostly protrusive forces due to the polymerization of actin and contractile forces due to the actomyosin contraction.Here, we don't model explicitly the actin network, but, like in [47] we consider the family of the Rho-GTPases as a regulatory signal for the activity of actin.In our approach, we introduce two kinds of proteins, the Rac protein, responsible for protrusive forces and the Rho protein, responsible for contractile forces.The concentration of those proteins is driven by diffusion processes inside the cytoplasm.
The main goal of this model is to understand the action on cell migration of another element of the cytoskeleton, the microtubules.They are polymers of tubulin that exhibit a very dynamic behavior.They quickly alter between phases of growth and shrinkage.In this model, the activation of Rac and Rho is regulated by the dynamic state of microtubules.Growing microtubules activate Rac whereas shortening microtubules activate Rho (See [32,49,51] for more details).

Determination of the cell speed
As mentioned before, the model is a fluid based model using the immersed boundary method.We assume that the cell and its environment can be seen as a fluid and we use fluid mechanics tools to describe its motion.
In the context of cell migration, we can estimate the characteristic parameters of the fluid motion in particular its Reynolds number: The low Reynolds number supports the use of Stokes equation to describe the motion of the fluid.We denote by u, the velocity of the fluid and p its pressure.
We consider three types of forces localized at the membrane of the cell, protrusive forces F pro , contractile forces F cont and elastic forces F el .Protrusive forces model the action of the actin polymerization on the membrane whereas contractile forces correspond to the strength tension of the actomyosin system.They are induced by biochemical reaction and their modeling will be described in Section 2.3.Elastic forces model the elastic properties of the cell membrane, their expression depends strongly on the representation of the membrane.
The Stokes equations are defined on R 2 and the velocity u is assumed to vanish at infinity.

The level-set method
An important part of the immersed boundary method consists in the evolution of the boundary, the cell membrane in our case.Different techniques have been used to handle the displacement of the boundary.The most commonly used are, the Lagrangian Marker Points technique (LMP), like in [47], the Volume Of Fluid and the Level-Set method.In this work, we have preferred the use of the level-set methods that offers a full eulerian formulation for the mechanical model which is very convenient for the numerical schemes.Moreover, the framework of level-set allows to handle easily large deformations of the membrane.
The method consists in representing implicitly the cell membrane as the zero level set of a certain function φ.The evolution of the function φ is then computed through a transport equation: Usually, the level-set function is initialized as the signed distance to the membrane.It takes positive values outside of the cell and negative values inside the cell.This choice, besides offering an improved numerical precision, gives to the level-set function an additional meaning.It is notably used in order to locate the forces at the membrane.The forces are regularized in the numerical scheme in a neighborhood of the membrane.That can actually be computed with the level-set function itself.More details about the Level-set technique can be found in [24,39].Elastic force: Elastic force used in this model depends on the local variation in space of the local stretching of the membrane.The expression of the local stretching of the interface depends strongly on its representation.In [47], they used an LMP model which makes very easy the computation of both the local stretching and its variations.Therefore, their works on elastic forces cannot be immediately transposed in the level-set context.
We use instead the expression proposed in [14].They have shown that the local stretching of the interface can be computed with the quantity |∇φ|.But despite this expression, the local variation of the stretching is hard to compute due to the implicit representation of the interface.To overcome this difficulty, they propose in [14] to derive the elastic force from an elastic energy for the membrane.We present here the regularized form of this energy, that is the one used in the numerical schemes: where λ el is the magnitude of the force and the term ) is an approximation of the Dirac function.Renormalization: A key component of the level-set techniques is the redistanciation step.It is well known that the level-set function does not remain a signed distance function at all time.It is necessary, to maintain the information of distance and localize the forces to regularly reset the level-set function into a signed distance function.The idea is to replace a given level-set function by a signed distance function that has the same zero level set.To do that, there exist many techniques, one can cite for example the use of the redistanciation equations introduced in [44].The main problem in our case is that the level-set function does not encode only the distance to the membrane but also keeps track of the stretching of the membrane.Hence we can't reset the level-set function without loosing all the information about the stretching.As an alternative, we choose to use the method renormalization technique proposed in [12].It consists in using the quantity φ |∇φ| as a first order approximation, near the interface, of the signed distance.This method is quite easy to implement and has been proved to be efficient in [13,14].

Biochemical model
The biochemical model is the part of the model that describes the polarization process of the cell.The Rac and Rho concentrations are the main components of this model and determine the protrusive and contractile forces that allow the cell to move.The microtubules are the key element of the Rac and Rho evolution, they regulate the activation rates of the proteins.

Rac and Rho equations
We consider the protein Rac and Rho on both active and inactive forms.The concentrations of the active forms will be referred as Rac and Rho respectively whereas the inactive forms will be named Rac and Rho.
For the evolution of the proteins, we use reaction-diffusion-advection equations like in [47]: with g P = −g P .The equations are posed on Ω(t) the interior of the cell which is determined with the level-set function φ.This domain is moving over the time at the speed u given by the mechanical model.The advection term in the equations takes into account this motion.These equations are endowed with Neumann boundary conditions on ∂Ω(t).
The term g P models the activation or inactivation of the proteins.Unlike [47], this term will be determined by the dynamical state of the microtubules that we are going to describe in the next section.The coefficients D Rac , D Rac , D Rho , D Rho are the diffusion coefficients of the proteins.

Microtubules modeling
The microtubules are, as mentioned before, very dynamic elements of the cytoskeleton.They are polymers of tubulin that alternate very quickly between phases of growth and phases of shrinkage.This continuous alternation is termed microtubule dynamic instability and characterizes their roles inside the cell.When a microtubule is polymerizing, it activates the Rac protein whereas when it depolymerizes, it activates the Rho protein.The specific process that regulates the dynamic instability is quite complex.More details about this phenomenon can be found in [19,32,49,51].
We introduce new unknowns to the model.First we consider the concentration of Tubulin, Tub, inside the cell.It will be used as material for the polymerization of microtubules.The evolution of the tubulin is modeled, like for Rac and Rho, by a reaction-diffusion-advection equation: The reaction term g Tub models a consumption term and will be detailed later.
The role of the microtubule in cell migration is completely determined through its dynamical state.Hence, we only need to locate the (+)-end dynamic extremity of the microtubule.Let N MT be the number of microtubules present in the cell.We will denote by: MT i , i = 1, . . ., N MT the position of the microtubule (+)-end.We assume that microtubule polymerization speed γ pol only depends on the tubulin concentration.This assumption has been demonstrated for a long time and is well documented in [35] for example.We assume here that the dependence is linear: with α pol the magnitude of the speed and c c the critical concentration at which polymerization speed and hydrolysis rate equilibrate [5,30,50].
The experimental observations of microtubules in a moving cell suggest that the microtubules are globally aligned in the direction of polarization of the cell, with many microtubules oriented to the front and fewer oriented to the rear.We define the direction of polarization, v pol of the cell as a combination of the velocity of the cell u and the gradient of the protein Rac inside the cell: From a numerical point of view, the schemes seem to be more stable when we rely more on u than ∇Rac so in practice we always use α > β.From a biological poitn of view, polarization is mostly due to the concentrations of proteins like CdC42 that are not considered in this model.We distinguish two types of microtubules, the ones oriented to the front of the cell and the ones oriented to the rear.In the following, we will call them MT + i and MT − i respectively (see Fig. 1).The direction v pol describes the main direction for the ones going to the front whereas −v pol is the main direction for the ones going to the rear.We use as a secondary direction v sec for the microtubules the gradient of free tubulin in the cell: This choice means that the microtubules go in direction that maximizes their polymerization speed.The direction of polymerization v ± dir of the microtubules is a combination of the main and secondary direction: with ε a regularizing parameter.The notation v + dir (respectively v − dir ) stands for MTs going to the front (respectively to the rear).The direction ηv sec ± v pol is used as the direction of polymerization as well as the direction of depolymerization.The actual direction of depolymerization is the direction of the microtubule itself.Because we don't keep track of the whole microtubule but just the (+)-end, we assume that the direction of depolymerization is the opposite of the direction of polymerization.This assumption needs that both the quantities u, ∇Rac and ∇Tub have a slow dynamic with respect to the growth lifetime of microtubules, which is the case in practice.Then the equation for the motion of the (+)-ends of microtubules is: where the term u reflects the fact that the cell is moving and the term v ± dir models the microtubule dynamics.The question that remains at this point of the modeling process concerns the behaviors of the microtubules near the cell membrane.In equation (2.11), the microtubule is not aware of the presence of the membrane and a numerical discretization of this equation leads to microtubules able to cross the membrane.In a moving cell, different behavior for the microtubules at the membrane can be observed.Some microtubules stay to the membrane for a while, others keep polymerizing along the membrane but most of the microtubules seem to depolymerize after reaching the membrane.The observations suggest that those microtubule then depolymerize up to the beginning of the lamellipodium.A biological hypothesis to explain this behavior could be that when reaching the membrane, the microtubules are destabilized and then enter in a shrinkage phase.
To reproduce this behavior, we define two populations of microtubules, the "stable" one, referred as MT ± i [Stab], and the "unstable", referred as MT ± i [Pert] that has experienced a contact with the membrane.We assume that the microtubule in an unstable state is depolymerizing and so its direction is the opposite of the direction of polymerization: ) A microtubule can pass from one state to another when reaching a given zone inside the cell.The transition from Stab to Pert is done when reaching the membrane whereas the transition from Pert to Stab is done at a distance d Stab from the membrane of the size of the lamellipodium (see Fig. 1): In the following we will keep using the notation MT i when talking of a microtubule without stating if it is a stable or unstable one.

Reaction terms
We define now the reaction terms in the equations (2.5) and (2.6) that depend on the dynamic state of the microtubules.
The tubulin is attached to polymerizing microtubules and released by depolymerizing microtubules.Then we set for the reaction term: with d the density of polymerized tubulin and δ 0 the Dirac function.
The reaction terms for the Rac and Rho protein are done as in [47], with a positive feedback on the activation.This positive feedback is here to maintain the polarization of the cell over time.The main change in the model is that Rac protein, activation takes place only near a polymerizing microtubule.The inactivation of the protein is modeled as a competition between the Rac and Rho proteins [42].Thus, in our model when a microtubule is depolymerizing, the Rho protein is activated leading to a Rac inactivation.The reaction terms then read: with, τ Rac→Rac the intrinsic activation rate, τ Rac→Rac the inactivation rate, γ Rac , K Rac parameters for the positive feedback.The function k 0 is used to know if a microtubule is polymerizing or not: The function 1 B(MTi,dMT) is the indicator function of the ball of radius d MT centered in MT i .The distance d MT is the radius of influence of the microtubules on the proteins.The function G models the influence of microtubule dynamic instability on Rac activation.We assume that the activation occurs only in the neighbourhood of a MT in polymerization ( vdir = 0) at a rate that increases with the polymerization speed γ pol : The reaction term for Rho is similar, exchanging the role of polymerization and depolymerization: with the parameters defined as for the Rac protein and G * (v ± dir ) = ||v ± dir || independent of γ pol as mentioned before.

Protrusive and contractile forces
The protrusive and contractile forces are defined by the action of the Rac and Rho proteins.We follow the idea of [47] to define the protrusive force, only taking into account the level-set definition of the membrane: The term ∇φ |∇φ| is the outward pointing normal to he membrane.The function h Rac define the Rac dependency of the force.We use a force with inhibition like in [47]: The constant H Rac is the magnitude of the force, Rac c the minimal concentration needed to observe a force and Rac + the concentration that generate an optimal force.The contractile force is defined in a similar way: with h Rho defined as h Rac .

Numerical model of cell migration in DDFV settings
The discretization of the different equations of the model is made with Finite Volume techniques.Because the cell moves during the simulation and to avoid boundary effect, we must take a computing domain large compared to the size of a cell, which can be very costly.
Locally refined meshes are commonly used in level-set problems.The main idea is to refine the mesh near the interface and to use a coarser mesh away from the interface.As the interface moves, we can use dynamical refinement techniques to adapt the mesh to the moving boundary.Here, we will use a less complicated approach by refining the mesh on a zone where the cell will stay during the whole simulation (see Fig. 2).The zone to refine is, for the purpose of our simulations, predicted depending on the initial polarization of the cell.We highlight the fact that we cannot compute the solution only on this refined zone because of the imposed boundary conditions on the velocity.In particular, a small computing zone will lead to large errors coming from the Stokes equations.We choose to use DDFV schemes that have been successfully used to solve Stokes equations on very general meshes like locally refined meshes.
The discretization of the complete model is made according to the following strategy.At each time step, we solve: -The Stokes equations discretized using a DDFV scheme as in [33].
-Reaction-diffusion-advection equations discretized using a splitting technique, between a diffusive part and an advective part.The method is based on the one proposed in [47].We propose an improvement of the method based on DDFV that allows a sharper discretization, in particular in the context of high deformations of the domain.-Transport equations (both for the previous splitting scheme and for the Level-Set equation) discretized using an original DDFV formulation of a WENO scheme ( [31]) that we have developed for the use of this model.

The DDFV discretization
Discrete Duality Finite Volume (DDFV) are Finite Volume methods introduced first in [17,28].They consist in a decomposition of the computing domain Ω in a set of polygons.Those polygons form the primal mesh and one unknown is associated to each polygon.Then additional unknowns are added on the vertices of the  polygons.Those vertices are therefore seen as centers of other polygons that define a dual mesh as in Figure 3. Introducing new unknowns allows to compute a discrete gradient in any directions, in duality with the classical finite volume divergence.
Let us recall notation for DDFV scheme introduced in [3].We denote by K a polygon of the primal mesh M and K * a polygon of the dual mesh M * .In the following, we will refer to an element of the primal or dual mesh by C ∈ M ∪ M * .We will also denote by ∂M the set of the edges contained in the boundary of the domain ∂Ω, which can be seen as degenerated cells and by ∂M * the cells of the dual mesh associated to vertices x K * ∈ ∂Ω.We name T = M ∪ M * the global mesh.A third mesh, the diamond mesh D, can be associated to the DDFV structure.We define the diamond cells D σ,σ * being the quadrangles whose diagonals are a primal edge σ = K|L = (x K , x L ) and a dual edge σ * = K * |L * = (x K * , x L * ) (see Fig. 4).
Hence, one diamond can be associated to each edge of the primal and dual mesh to construct the diamond mesh like in Figure 3.This set of 3 meshes can seem to be complex but in practice only the diamond mesh is used for the implementation of the mesh and needs to be constructed.
DDFV methods have been used to discretize a large class of differential operators, like Stokes or non linear diffusion.The complete discretization of the gradient in this method presents the advantage to not imply geometrical dependencies to the mesh that we choose to use.Consequently, in order to deal with locally refined meshes and to release us from the orthogonality constraint imposed by the classical Finite Volume methods (see [21]), we choose to use a DDFV strategy.
We briefly recall here the complete definition of the gradient and divergence operator that we use in the following.For a discrete scalar field u T ∈ R T , we set the discrete gradient defined on each diamond D by: where the notation |.| stands for the measure of the corresponding element (edge, cell, diamond,...).The discrete divergence is defined for a vector field ξ T ∈ (R 2 ) D on each cell of the mesh by: with: where we define D C = {D ∈ D, D ∩ C = ∅} for a given cell C and D σ,σ * as the diamond defined with the two diagonals σ and σ * .
We can also define a discrete divergence div D on the diamond mesh:

Discretization of the Stokes equations
The Stokes equations of the mechanical model are discretized using the DDFV scheme developed in [33].
The method consists in a staggered scheme, where the unknowns for the components of the velocity are located on the nodes of the both primal and dual meshes, whereas unknowns for the pressure are located on the diamonds.The equation (2.1) is discretized on the primal and dual meshes and the equation (2.2) is discretized on the diamond mesh D.
In order to have a well-posed and stable scheme, we use a stabilization of the scheme by a Brezzi-Pitkäranta technique.The stabilized DDFV scheme for the Stokes equations can be written: with strong Dirichlet conditions: ∀C ∈ ∂M ∪ ∂M * , u C = 0.More details about the discrete operators can be found in [33].

The force terms in the Stokes equations
The originality of this work stands here on the discretization of the force terms.
Protrusive and contractile forces: The protrusive and contractile forces are discretized in the same way.We integrate the force on each cell C ∈ M ∪ M * , then, we decompose the integral into integrals on the semidiamonds, which give for example for the protrusive force: The unknowns Rac and φ are defined at the nodes of primal and dual meshes.We can then compute a discrete gradient of φ on each diamond.We approximate the gradient of the semi-diamond by the discrete gradient on the diamond.The function h Rac is evaluated at the center of the cell C. Then the complete discretization is given by: The discretization is equivalent for the contractile force.
Elastic force: The divergence expression of the elastic force makse it easy to discretize in a Finite Volume context.We use here the discrete divergence div T on the primal and dual cells as well as the discrete gradient ∇ D defined on the diamonds.The only difficulty remaining concerns the localization term which involves a value of φ on the diamond.Here, in order to maintain the conservative form of the force, we choose to use a weighted mean value of φ on the diamond: with d C,σ the distance between x C and x σ .The discretization of the elastic force is then given by:

Discretization of the diffusion equations
Consider the reaction-diffusion advection equation: with Ω(t) a moving domain with a divergence free speed u.
The main difficulty in the discretization of this equation lies in the time dependence of the domain Ω(t).We follow and adapt the idea proposed in [47] to use a splitting method by separating the transport process and the diffusion process.They propose to first solve a diffusion equation and then to transport the cell with a Lagrangian method.Here we will present an improvement of this method that consists in reducing the problem to a non-linear diffusion equation on a fixed domain.More details about this approach can be found in [45].

Lagrangian change of variables
The first step to the splitting technique is to solve the equations on a fixed domain.The most natural change of variable consists in using the characteristic curves of the flow.
Let us introduce Φ(.; t, x) solution to: The function Φ(s; t, x) represents the position at time s of the particule that was at the position x at time t.
Then we set the change of variables: The function F defines a bijection from Ω(t) on Ω(0).We recall that the change of variable in the divergence operator is given by the following formula.
Theorem 3.1.Let f be a regular vector field: Here, as div(u) = 0, the flow preserves the volume and then: Equation on a fixed domain: That provides us a formulation in the X coordinates for the reaction-diffusionadvection equation: with ã defined by ã(t, X) = a(t, Φ(t; 0, X)).
We can note that the advection term is no longer present in the equation (3.16) but we have introduced an anisotropic diffusion matrix ∇ X Φ(t; 0, X) −1 ∇ X Φ(t; 0, X) −T .Note that the hypothesis of ∇.u = 0 is important to obtain an equation in a divergence form.Boundary conditions: In the model we consider Neumann boundary conditions.It can actually be shown that we find back Neumann boundary conditions after the change of variables: with n and n 0 the outward pointing normal to respectively ∂Ω(t) and ∂Ω(0).Anisotropy: The anisotropic diffusion matrix introduced in the equation (3.16) keeps tracks of the deformations of the domain.One can easily ensure that in the case of non-deformating move, like a translation or a rotation, there is no anisotropy, meaning ∇ X Φ(t; 0, X) −1 ∇ X Φ(t; 0, X) −T = Id.In the general case however, when the domain is actually deformed, the diffusion is anisotropic.

Splitting method
The previous change of variables shows that it is equivalent to solve the equations (3.10) or (3.16).The equation (3.16) is defined on a fixed domain and so is easier to discretize and solve numerically.But because we are interested in the function a and not ã, we must find a way to go back to the x-coordinates after solving the equation (3.16).One difficulty is that we don't have access to the function Φ explicitly.But we can use the definition of the characteristic curves to observe that knowing the function ã we can solve a transport equation to find back a.
The second part of this splitting technique is solved using the same scheme as for the level-set equation.This scheme is detailed in Section 3.5.

Numerical scheme for the anisotropic diffusion
The DDFV method is known to be an efficient way to solve anisotropic diffusion in a finite volume framework.The discretization of the equation using classical finite volume method implies an orthogonality condition on the mesh which involves the anisotropic operator and which can be hard to satisfy in practice.The DDFV method allows to discretize those schemes in a large class of meshes.
In the following we set: the anisotropic operator.
The time discretization of equation (3.16) is done with an Euler scheme.We use an implicit discretization for the diffusion and an explicit discretization for the reaction.
Let h rd be the time step of the scheme, we note ãn ) D∈D an approximation of M (t, X) on each diamond D at time t n .Description of this discretization is given in Section 3.4.4.
The implicit DDFV scheme for 3.16 is given by: The Neumann boundary conditions are no flux conditions:

Discretisation of the anisotropic operator
In the general case of our model, we don't know explicitly the anisotropic diffusion operator and so we can't discretize it easily.However, we know an equation for the gradient of the characteristic curves that we write thereafter in an integral form: Using a left rectangle quadrature rule on (3.22), we obtain: To discretize M , we need first to discretize the inverse of ∇Φ(t; 0, X).Assuming nh rd to be small enough, we can use a geometrical series troncated at order p: We do the same for the transpose.Taking p = 1 and neglecting the order 2 term in h rd we set: This provides us a discretization of the anisotropic diffusion term.One can think to use either a more complex integration rule for equation (3.22) or using a forward Euler method in the differential equation defining ∇Φ(t; 0, X).The main difficulty are the term ∇Φ(t; 0, X) in the right side of the equation and the term Φ(t; 0, X) inside the gradient of u.Because, in our model, u is not known explicitly, those solutions seem quite difficult to use in practice.
Remark 3.3.Due to CFL conditions, the transport equation imposes the strongest constraint on the time step.Then in practice we will usually choose h rd = T k+1 − T k , meaning we only do one step of diffusion for multiple steps of transport.
Remark 3.4.The scheme proposed in [47] can be derived formally in the same way using a semi-implicit version of equation (3.20)where the anisotropic diffusion operator is explicit: Then, using this semi-implicit scheme and the choice of time step proposed in Remark 3.4 we find back the scheme proposed in [47].The complete derivation of the scheme presented in [47] can be found in [45].

Submesh for the diffusion equation
In our simulation we use a fixed mesh of a domain Ω such as Ω(t) ⊂ Ω∀t.Hence, at each time t k , we must localize the sub-domain Ω(t k ) and solve the reaction-diffusion equation only on this sub-domain.This can be done using the level-set function.Let us define: We define then the dual mesh associated with this submesh: The submesh used to solve the diffusion equation is then M k ∪ M * k .

Discretization of the transport equations
Transport equations are present in the model both for the level-set technique as well as in the solving of the diffusion equations.To solve the equation, we use a high-order WENO scheme.WENO schemes are commonly used in the level-set community.High-order schemes are needed for the transport of the level-set function in order to conserve the volume of the cell.Other high-order schemes could be used for transport equations.One can cite the MUSCL schemes [40] or the MOOD schemes [10].
First introduced by Harten, Osher and others [25][26][27]43], WENO schemes are known to be efficient on convection problems.The WENO strategy consists in computing a non-oscillating high degree polynomial interpolation of the solution on the edge of the mesh.
A complete numerical study of this scheme can be found in [45,46].

Time discretization
The transport equation on a bounded open set Ω ⊂ R 2 , with a divergence-free velocity u, can be written as: For the time discretization, we follow [22,43] and use a TVD Runge-Kutta of order 2. Let ∆t be the time step of the method, we will denote by φ n the approximation of function φ at time t n = n∆t.The RK2 scheme is then given by the following steps:

Discretization of operator div(φu)
Let φ τ = (φ C ) C∈M∪M * , the vector of the approximations φ C of the mean values φC = 1 |C| C φ of function φ on the cells C ∈ M ∪ M * that we want to compute.
Following the Finite Volume strategy, we integrate the operator L on each cell: where n is the outer unit normal to the boundary ∂C of C.
The curve integral can be approximated using a p point Gaussian quadrature.In practice, we use p = 2 in our numerical simulations.
The WENO scheme consists in approximating for each cell C and each edge σ, the value φ(x σ ) by a convex combination of the value in x σ of several polynoms whose mean values coincide with the mean values of φ on a set of selected cells.This set of cells is called the stencil of the method.Let us assume that we dispose of such an approximation φ C,σ .Then we define the spatial discretization L C of operator L using an upwind flux as: where C and C share the edge σ and φ b prescribed through the boundary condition.Let define φ n,τ = (φ n C ) C∈M∪M * the vector of the approximation φ n C of the mean value of φ on the cells C at time t n .The full discretization is then given by: This discretization is done in the same way on both primal and dual meshes.The coupling between the two meshes is ensured by the reconstruction process.

Approximation of φ C,σ
Given a cell C and an edge σ, we want to reconstruct an approximation of φ(x σ ) though we only know the mean values ( φC ) of φ on M ∪ M * .Following the WENO strategy, the approximation φ C,σ is computed as a convex combination of several polynomial interpolations of φ.
To find those polynomial interpolations, we fix a subset S ⊂ M ∪ M * , depending on C and σ, and we choose the polynom P S [φ] among the polynoms of degree 2 as the solution of the following problem: The weights in the convex combination are chosen in order to favor non-oscillating polynoms: The coefficients a S are choosen in order to avoid oscillating polynoms.We choose to use the oscillating criterium for the polynoms provided by Abgrall (c 0 ) in [2].This criterium is defined for a polynom P = |α|≤m p α X α as: Other criteria and weights can be found in [22] and [29].Then, we set the weights a S as in [22,31]: The stencils are associated to a couple cell-diamond (C, D).They are chosen randomly among the unknowns in a neighbouring of the diamond.The complete process for the stencil selection is presented in [46].

Numerical illustrations of a MTas effects on cell migration
We present now numerical experiments to illustrate the behavior of the model and the influence of the different parameters of the model.These simulations exhibit the importance of the microtubules dynamics to maintain the cell polarisation and to drive cell migration.
In particular, some details are presented about the modeling of the action of microtubules targeted agents on cell migration.
During in vitro experiments, several indicators are used to describe the microtubules dynamics as well as the migration of the cells: -growth speed: the mean polymerization speed of microtubules.
-growth lifetime: the mean duration of growth phases of microtubules.
-growth distance: the mean distance covered by microtubules during polymerization phases.
-migration distance: the distance covered by the cell.
-migration speed : average speed of the cell.
In the numerical simulations, we can compute the same quantities in order to compare them with the biological data.That provides us quantitative criteria to estimate the global behavior of the model.

Calibration of the parameters
One of the crucial step of the modeling process is the calibration of the different parameters of the model in order to simulate the migration of a given cell type.To do so, we use some data provided by the literature as well as the data provided by biological experiments done in the Institute of Neurophysiopathology, in Marseille, on U87 glioblastoma cells and designed especially for this purpose.Details about those experiments and the calibration process can be found in [16,45].
We just recall here the different values for the parameters provided by the calibration process.
We detail here the settings of the initial condition for the reference test.In the following we will only discuss the change of settings regarding this test called reference test.
The initial form of the cell is a circle of radius 5µ m.The number of microtubules is settled to 45 at the front of the cell and 15 at the back.At t = 0, the tubuline concentration is uniform inside the cell with a concentration with the parameters used for each initial concentration presented in Table 2.
All simulations are performed on 3 min, with a time step dt = 0.01.The mesh used for the simulations is shown in Figure 5, with space step of dx = dy = 0.46875 µm.
The shapes and positions of the cell over time are presented in Figure 6.The indicators of the migration introduced previously are computed and confronted to data experiments from [16,45] in Table 3.
We can also compare the polarisation of the cell between t = 0 and t = 3 min.Results are shown in Figure 7.We observe that the microtubules maintain and even reinforce the cell polarisation during migration concomitantly with an activation of RAC (maximum value of RAC increasing from 7 µM to 8.5 µM).

Parameters involved in a drug effect
In the sequel, we focus more precisely on the influence of the parameters involved in the modeling of an antimigratory effect by MTAs.Five important mechanisms, that can potentially be modulated by the action of drugs, have been identified in the cell migration process: the number of active microtubules, cell polarisation, MT instantaneous polymerization speed (γ pol ), the Rac activation rate and the tubulin-GTP hydrolysis rate as shown in Figure 8.Those main regulating processes have been chosen based notably on the experiments of E. Denicolai in [16].Those experiments show in particular that the action of a drug like vincristine affecting microtubules dynamics mainly reduces Rac activation and doesn't seem to act on the activation of the Rho protein.Indeed, we choose to consider only mechanisms linked to the Rac protein.
The hydrolysis rate in Figure 8 describes here the processes that allow to define the equilibrium of the MT dynamic instability.The dynamical state of the microtubules depends in particular of the balance between two processes occuring at the (+)-end of the microtubule: the hydrolysis of the GTP-tubulin into GDP-tubulin and the incorporation of new GTP-tubulin.This equlibrium is modeled in a simplified way here by the definition of the critical concentration c c (see [30]).
From a modeling point of view, modulating those mechanisms consists in changing respectively the parameters: N MT , α Rac , the function G (see Def. (2.18)), τ Rac→Rac , c c .We could choose other parameters in order to affect those mechanisms but we select here a minimal set of parameters that seems biologicaly relevant.

Impact of the initial cell polarisation
The cell polarisation strongly impacts the migratory behavior of the cell.This influence can be seen in the model by the influence of the initial polarisation.Even if microtubules have an effect, in this model, on the     biochemical state of the cell, the polarisation process is not modeled here.We check that by simulating the migration of the cell without initial polarisation, i.e. setting the initial concentration of Rac and Rho to uniform concentrations.Then we can observe in Figure 9 the cell does not migrate between t = 0 min and t = 3 min, even if the set concentration is enough to generate forces (Fig. 10).The polarisation of the cell also remains constant during the whole simulation and without polarisation the microtubules stay at rest (Fig. 9).

Influence of the number of microtubules on the model
We are interested in this model by the impact of the microtubules on cell migration.To highlight the role of microtubules in the model, we perform several tests with a different number of microtubules.First, the reference  test shown previously was done with 45 microtubules at the front and 15 at the rear.We compare it to two cases.We remove all the MTs in the first one, and take a small number of MTs in the second one (20 at the front and 5 at the rear).
Figure 11 shows that without microtubules the cell does not move at all and that decreasing the number of MTs slows down the migration.The migration speed is 1.1 µm min −1 compared to 1.7 µm min −1 in the reference test.Note that we also have tested higher number of MTs, which does not lead to strong increase in migration speed but affects the quality of the simulation process.
With these simulations (Figs.11 and 9) we can see that the initial polarisation is very important to create cell motion but also that microtubules dynamics actually control cell polarisation during motion.

Impact of the speed of polymerization
In the activation term used for Rac protein, we use a term G (2.18) that models the influence of the microtubule polymerization speed on the activation of the Rac protein.Here we support that choice by comparing two sets of simulations.In the first set, we use the function G as mentioned in Section 2.4 and in the second one, we use instead G * (γ pol , v ± dir ) = ||v ± dir || where we skipped this influence.We perfomed simulations for both functions G and G * for the reference test proposed in Section 4.1 and for a second one corresponding to α pol = 6 µm min −1 .µM−1 .In this particular case, the reference test can be seen as a vincristin effect compared to the second one, as vincristin is known to suppress MT growth speed.We compare the different indicators of the MT dynamics and cell migration in Table 4.
It seems that in absence of influence of G on microtubule polymerization speed in the activation rate, we are not able to show an effect of the microtubule polymerization speed on the cell migration.This is due to the difference of time scale between microtubules dynamics and cell migration.If the activation only depends on the dynamic state, then at the time scale of the cell, the protrusive force is mainly due to the time spent by polymerizing microtubules near the membrane which is the same for slow and fast microtubules.

Impact of the Rac activation rate
We show here that a modification in the activation of Rac protein can lead to a change of cell migration.We perform a new simulation taking τ Rac→Rac = 0.1 µM min −1 .The results are shown in Figure 12.The cell goes faster with a strong activation of the Rac protein.We observe in Figure 12 a delay of about 2 min in the cell migration in the case of a low τ Rac→Rac = 0.1 µM min −1 .
The activation of the Rac protein is an important parameter that can be modulated in order to simulate a drug effect.

Impact of the GTP-hydrolysis
The microtubule polymerization is regulated by the balance between the hydrolysis of the GTP-tubulin and the incorporation of new GTP-tubulin dimers.This equilibrium is characterized, in our model, by the parameter c c (see [30]).As drugs like vincristine are known to impact this equilibrium by increasing the c c , we propose to investigate the influence of this parameter on the dynamics.We perform a simulation taking c c = 5 µM instead of 4 µM in the reference test.The results of the experiments are shown in Figure 13.We can observe that increasing the parameter c c leads to a total inhibition of cell migration in this case due to a too low Rac activation.We observe at t = 3 min a maximum value of RAC of 6.5 µM against 8.5 µM for the reference test.
The results of this simulation show optimistic results in the modeling of MTAs action on cell migration.A specific study of the effect of vincristine and its modeling is presented in [16].

Conclusion
In this paper, we proposed a novel modeling approach to describe how MT dynamics may impact cell migration and developed numerical tools to efficiently solve the equations of this model.The model was developed from a previous model of the literature [47] describing the effect of Rho-GTPases on cell migration.Microtubules acting as a catalyst for the activation or inactivation of these proteins, we have coupled the previous model to a relatively simple model of the microtubule dynamic instability.It describes on the one hand the mechanical aspects of migration with a submerged boundary method using Stokes equations and a level-set technique for monitoring interface involving transport equations.On the other hand, it describes the biochemical aspects of the cell with mobile domain reaction-diffusion equations for Rac, Rho and tubulin.Microtubules are described through the localization of their (+)-end.They act as catalysts for the activation of Rac and Rho proteins through their dynamic state.The resolution of the different equations of the model has been implemented by DDFV techniques allowing to use locally refined meshes around the cell.The use of a splitting method for the resolution of the reaction-diffusion equations on a mobile domain was motivated by using a suitable change of variables.We show that the resolution is equivalent to solving an anisotropic diffusion equation on a fixed domain followed by the resolution of a transport equation.A DDFV version of a high order WENO scheme for the resolution of transport equations has also been proposed.
Numerical illustrations performed were callibrated according to the experiments described in [16] in such a way that our simulated microtubule dynamics, cell shape and cell migration indicators were biologically relevant.Moreover, we were able to analyse the impact of several parameters identified as a crucial target of MTAs.According to our results, our approach reveals to be robust to investigate MTAs effect on cell migration.This paper opens large perspectives.At the level of modeling, it would be interesting to use a more precise model for the MT dynamic instability, with the implication of GTP-tubulin hydrolysis and/or EB proteins [5,38,50], and/or to include the effect of Cdc42 to regulate the initial polarization of the MT cytosqueleton and protusion and thus the cell polarity [7,8,34,52].
From a numerical point of view, we are still limited by the numerical cost of the scheme.In order to be able to perform numerical simulations on larger times, it will be necessary to improve the WENO procedure and use refinement/derefinement procedures to reduce the global size of the mesh.
Concerning the confrontation of the model with the biological experiments, the calibration of the parameters is a long process described in [16].The next step will be to identify the dose response effects of drugs on the parameters of the model.Another development track concerns also the action of other MTAs such as paclitaxel, a stabilizing agent or new MTAs such as BAL 101553 [6].

Figure 1 .
Figure 1.A cell with the microtubules MT + i and MT − i represented (left).Example of tracks of microtubules with the corresponding change of state (right).

Figure 2 .
Figure 2. Example of a cell and the direction of its move (left).Corresponding mesh refinement used for the simulations (right).

Figure 3 .
Figure 3. DDFV structure.From the left to the right: primal mesh, dual mesh and diamond structure.

Figure 4 .
Figure 4. Example of a diamond with corresponding notations.

Figure 5 .
Figure 5. Localy refined mesh used for the simulations.Cell membrane is shown in black, distances are in µm.

Figure 6 .
Figure 6.Time sequence of cell shapes over 3 min.Positions and shapes are shown at t = 0, t = 1 min and t = 3 min.

Figure 7 .
Figure 7. Concentration of active Rac protein at t = 0 min (left) and t = 3 min (right) for the reference test.

Figure 8 .
Figure 8. Schematic view of the dynamics: links between MT dynamics and Rac activation.

Figure 9 .
Figure 9. Shape, microtubules and Rac concentration at t = 0 min (left) and t = 3 min (right) in the case of no polarisation.

Figure 10 .
Figure 10.Vector field of the contractile force at t = 0 in the case of no polarisation.

Figure 11 .
Figure 11.Shapes of the cell at t = 3 min with no MT (blue), a 25 MT (black) and 60 MT (red).

Table 1 .
Values provided by the calibration process.

Table 2 .
Parameters for the initial linear concentration of proteins.

Table 3 .
Differents indicators of the migration for the control cell evaluated on a 1 min interval.The other concentrations of proteins are linear function of x of the form: