A Hybrid Mass Transport Finite Element Method for Keller–Segel Type Systems

We propose a new splitting scheme for general reaction–taxis–diffusion systems in one spatial dimension capable to deal with simultaneous concentrated and diffusive regions as well as travelling waves and merging phenomena. The splitting scheme is based on a mass transport strategy for the cell density coupled with classical finite element approximations for the rest of the system. The built-in mass adaption of the scheme allows for an excellent performance even with respect to dedicated mesh-adapted AMR schemes in original variables.


Introduction
The aim of the present work is to design a numerical scheme capable to deal with concentrations and diffusion phenomena typically arising in one-dimensional taxis-diffusion systems of the form with Lipschitz continuous source terms R ρ , R c that satisfy R ρ (0), R c (ρ, 0) ≥ 0.Here ρ denotes the cell density and c the concentration of a chemo-attractant.These systems constitute adaptations of the classical cell migration model by Patlak, Keller and Segel [29,20].They have been widely used in the modeling of biological processes such as cell organization in tissue, immune system dynamics and cancer growth [13,2,31].The dynamics of their solutions are quite rich; apart from traveling waves [23] the aggregation phenomenon studied in [16,5] that leads to blowup in finite time is of specific interest.One has moreover observed the occurrence of high concentrations that can emerge in a smooth solution, split, and merge with each other [27].Nonlinear diffusions or saturated responses in the chemotactic sensitivity are natural ways to include volume filling effects into the models, see [26,6].They usually avoid blow-up in a biologically meaningful way and lead to interesting phenomena and asymptotic stabilization.Finally, these models are basic building bricks for a variety of cancer invasion models in the literature [11,31,30,32,17] in which the coupling with extracellular matrix, enzymatic activators and other substances are taken into account.One of the common features in all of these models is the simultaneous occurrence of regions of high concentrated densities with diffuse profiles leading to numerical difficulties in choosing well-adapted meshes.The numerical approximation of all of these simultaneous phenomena is particularly challenging.
In [3] a mass transport steepest descent scheme has been proposed to resolve a modified 1D Keller-Segel system for the log interaction kernel proposed in [7].The method satisfies a discrete free energy dissipation principle by design being based on the variational schemes for Fokker-Planck type equations introduced in [18,21] and applied to Keller-Segel type models in [3,4].By considering the problem in transformed variables the method can resolve areas of high concentrations accurately without any mesh refinement.This approach has been extended to several dimensions for nonlinear aggregation-diffusion equations and with different approaches in the discretization in [8,24,9,19] and the references therein.
The aim of this work is to extend the mass transport approach to the general class of systems (1.1).We will test different scenarios that feature in particular the splitting, traveling and emerging of concentrations.
For the adjustment of the scheme we propose a splitting method, where we employ the technique from [3] to the Keller-Segel part of the system (i.e. the first equation of (1.1) with R c = 0).The remaining system of an ODE and a diffusion reaction equation will be decoupled and solved by a suitable finite element method.The advantage of the mass transport approach for the cell densities equations is that the mesh adapts naturally to the mass distribution, and then coarse meshes in the mass variable can still lead to good numerical approximations as we will discuss below.
In more details, we split (1.1) into two subsystems.The solution of the full system (1.1) can then be approximated by appropriately combining short time solution of the subsystems.We introduce at first the diffusion-advection system given by (I) This system makes the assumption of a steady chemo-attractant density c and mass conservation in the cell density ρ.Second, we consider the reaction-diffusion system that contains the remaining terms of the system.Following the mass transport algorithm [3] we transform the system (I) into new variables.With this aim, we consider the pseudo inverse cumulative distribution of the cell density ρ, V (t, w) = inf y : which is defined by The system (I) can now be rewritten, following e.g.[10], as where m denotes a given mass during this splitting step.The advantage of the proposed splitting is that the mass of cell densities does not change over the first step and the cell density is fixed over the second step.
The details of the full discretization of the proposed splitting scheme will be given in Section 2. In Section 3 we discuss the choice of the constraints in the time, spatial and mass steppings due to the choice of the full discrete schemes.Section 4 is devoted to study in detail the performance of this splitting scheme in many complex situations ranging from the simpler Keller-Segel type systems and their small variations to quite more biologically relevant systems in tumor invasion as discussed above.We will analyze the experimental convergence and the computational cost of this discretization with respect to previous schemes with meshrefinement algorithms in original spatial variables.Finally we conclude in Section 5.

Numerical method
In what follows, we describe a numerical treatment for both systems (I') and (II).The inverse distribution V is given on the time evolving mass space (0, m h (t)), whereas the chemo-attractant c is given in the Eulerian coordinates in (a, b).This leads to two meshes that the proposed numerical method employs.
First, we discretize the normalized mass domain (0, 1), on which the pseudo inverse distribution V resides by the mesh with length M ∈ N and width h w = 1/M that corresponds to the width ∆w(t) = m h (t)h w in the time evolving mass domain (0, m h (t)).We denote the point values of V by V j (t) = V (m h (t) w j , t) for j = 0, . . ., M and introduce the linear spline in w connecting the discrete values that we denote by V h (t, m h (t) w).Here we have used the discrete mass of the cells where ρ h is a discrete representation of the cell density to be defined later on.
A second mesh partitions the physical space (a, b) for the chemo-attractant density c into The chemo-attractant mesh is thus of length N and width ∆x = (b − a)/N .We employ a linear finite element representation for the chemo-attractant density c.Therefore let {φ k , k = 1, . . ., N − 2} be the basis of piecewise linear hat functions on the grid (2.1) satisfying the boundary conditions.In particular, we have in the center of the domain and otherwise near the boundary.By using the basis functions we can define the approximate chemo-attractant density as For the construction of the splitting method we define solution operators for both systems (I') and (II).
To this end we design T to be a numerical solution operator of system (I') in the following sense: if ) is a numerical solution of system (I') at time t = t + ∆t.In the same manner, we define also a solution operator S for system (II).

2.1.
The solution operator T for system (I').For a discretization of the system (I') we need to evaluate the derivative of the chemo-attractant concentration in the state variable V .With this aim we consider an interpolation by cubic splines of the discrete chemo-attractant concentration.Let (V h (t), c h (t), m h (t)) be given initial data.By ĉh we denote the cubic spline over the data points (x k , c h (t, x k )) for k = 1, . . ., N that satisfies the boundary conditions ∂ x ĉh (a) = ∂ x ĉh (b) = 0. We use this spline for the approximation of the advection term.Concerning the time integration we split the taxis and diffusion terms and treat the stiff diffusion terms implicitly.In this way we allow for both large time steps and stability of the scheme.We apply in particular the two stage implicit-explicit midpoint scheme (see e.g.[28]) that reads in our case both for j = 0, . . ., M .We have approximated the diffusion terms above by a central difference formula as in [3].At the boundary we impose Neumann boundary conditions, i.e. 1 The intermediate stage Ṽj (t) is given by a nonlinear implicit equation (2.2a) and we use the Newton's method for its computation.For the computation of the taxis terms in (2.2a) and (2.2b) we evaluate the afore determined spline ĉh .
The chemo-attractant density as well as the mass of the cells are not affected by system (I'), hence we define the numerical operator accordingly by Note that if instead of linear diffusion, i.e.D ρ constant, we have a power-law nonlinear diffusion D ρ (ρ) = D ρ ρ γ−1 , γ > 1, modelling cell volume size effects as in [26,6], we obtain a similar approximation with D(t) = D ρ γ −1 ∆w(t) γ−1 , j = 0, . . ., M and similar boundary conditions as above.Remember that the continuous function T ∆t V h (t) is built as the linear interpolant of the values T ∆t V j (t) for j = 0, . . ., M , and thus we can define a reconstructed density T ∆t ρ h (t) by its own definition as long as the sequence V j (t) is strictly increasing.
2.2.The solution operator S for system (II).In the splitting method that we propose we will apply the reaction-diffusion operator S starting with the data (T ∆t V h (t), T ∆t c h (t), T ∆t m h (t)) obtained from a previous evaluation of the operator T .For simplicity we will describe the numerical operator S for general initial data (V h (t), c h (t), m h (t)).
System (II) is formulated for physical concentrations of cells.To provide adequate initial data using the given approximations (V h (t), c h (t), m h (t)) we transform the discrete pseudo inverse distribution V h (t) on (0, m h (t)) to a finite volume representation of ρ(t, •) on (a, b).Since the approximate density ρ h satisfies Vj (t) for all j = 1, . . ., M by construction (2.4), we can introduce the cell averages and the piecewise constant function ρ h in the following way This approximation of the cell density resides on physical space (a, b).Note though that the cell averages are given on a non-uniform grid which differs from the grid for the chemo-attractant density c given in (2.1).Now, we are in the position to write down the scheme for system (II).Again we split diffusion from reaction and apply the implicit-explicit midpoint scheme and obtain As usual, we employ precomputed integrals of the basis functions For their computation we use suitable quadratures together with an indicator function to identify the position of a particular point x ∈ [a, b] on the grid corresponding to the cell density ρ h .The reaction update in the cell density c h alters the mass of the cells over the interval Ω.Thus we update m h (t) by To be able to apply the advection-diffusion operator after the reaction-diffusion update we transform S ∆t ρ h (t) to its inverse distribution representation S ∆t V j (t).Therefore, we use the formula S ∆t ρ j (t) χ (Vj−1(t),Vj (t)) (x) dx = S ∆t m h (t)h w , j = 1, . . ., M. (2.6) As long as S ∆t V j (t) is monotonically increasing in j, identity (2.6) allows for an efficient update of the inverse distribution V h .

2.3.
The splitting method.To approximate the full system (1.1) we propose the classical Strang splitting method [33] employing both numerical operators defined above.For given non-negative and sufficiently smooth initial conditions ρ 0 and c 0 of system (1.1) we deduce discrete initial data (V h (0), c h (0), m h (0)).To compute a discrete representation V h (0) of the normalized concentration ρ 0 /m h (0) we integrate as in (2.6).
Then we define the fully discrete Strang splitting scheme for system (1.1) iteratively by where 0 = t 0 < t n = n i=1 ∆t i is a discretization of the time axis.In this way we alternate between applying the diffusion-taxis and the diffusion-reaction operator.The symmetrical structure leads to the second order splitting error.
To optimize the efficiency we adapt the time increment ∆t in each time step.Since the discretization of system (I) is more sensitive to instabilities that are caused by large time increments ∆t than the discretization of the diffusion-reaction system, we start the method in each time step with the numerical operator T in which we determine ∆t n .We will elaborate on the stability of the scheme in Section 3.
The scheme (2.7) is not limited to the case of a single pair of a cell and an chemo-attractant.An extension to multiple attractants (i.e. a replacement of χρ∇c by a sum The case of multiple cell densities coupled through the taxis terms, such as in the model discussed in [31], can be treated as well.Note though that each cell species brings along another non-uniform mesh on the domain (a, b) which requires further projections in the numerical operator S.

Monotonicity preservation of the diffusion-taxis operator
As demonstrated in [12] unphysical negative values that arise in the numerical solutions of the Keller-Segel type systems can cause instabilities and wrong behavior of the scheme.Hence, the so called positivity preserving finite volumes schemes for these kind of models have been developed, e.g. in [12].A non-negative density ρ implies a monotonously increasing pseudo inverse distribution V by its definition (1.2).If a method operates on inverse distributions it should in turn preserve the discrete monotonicity of V .This monotonicity preserving property of such schemes was studied in the case of filtration and convolution-diffusion equations in [15,14].In more details, We call a method monotonicity preserving if from V j (t) − V j−1 (t) > 0 for all 0 < j ≤ M follows that also In the rest of this section we focus on a simplified problem that motivates a way to adapt the time increment ∆t in such a way, that non-monotone solutions and thus possible related instabilities are avoided.We consider in particular the system (I') for the case of a steady chemo-attractant c ∈ C 1 (a, b).For the numerical resolution we consider a forward Euler scheme of the form for a discrete inverse distribution as defined in Section 2. This scheme can be understood as an explicit firstorder version of the advection-diffusion operator introduced in the previous section.In this setting we can follow the lines of [15,14] and derive a bound on ∆t that makes the scheme (3.1) monotonicity preserving: Lemma 3.1.The scheme (3.1) is monotonicity preserving, if for a fixed θ ∈ (0, 1) both CFL conditions are satisfied.
Proof.We consider a single time step in the scheme (3.1) and drop the superscript n.For brevity we will use the notation ∆V j+1/2 = V j+1 − V j .We assume the monotonicity of the discrete inverse distribution at the time instance t and compute for an arbitrary 0 ≤ j < M the difference By applying the mean value theorem to the function f (x) = x γ /γ we find two function evaluations of its derivative, κ j and κ j+1 , such that we obtain .
Note that by the non-negativity of ∆V j+1/2 both κ j and κ j+1 are non-negative.In the next step, we define Finally we estimate by the monotonicity at time instance t By using the conditions (3.2a) and (3.2b), the non-negativity of the right hand side in (3.3) follows.This implies the monotonicity-preserving property of the scheme (3.1).
For our splitting method (2.7) we assume that we avoid time step restrictions due to the diffusion terms by our implicit treatment and take a closer look on the condition (3.2b) (θ = 0).The point values of the inverse distribution V j (t) for 0 ≤ j ≤ M coincide with the mesh cell interfaces of the non-uniform mesh corresponding to the cell densities ρ h (t).Thus the quantity L j+1/2 in the proof of Lemma 3.1 can be understood as a finite difference formula for the second derivative of the chemo-attractant density c.In effect, the above CFL condition (3.2b) motivates to choose the time increment ∆t n according to For our numerical experiments with the more complex scheme (2.7) we have accordingly computed the time increments by for constants CF L, K > 0. The additional bound proportional to ∆w balances the temporal and the spatial errors; large values of K can be used in practice.We have chosen CF L = 0.49 and K = 100 in our numerical experiments.Using this condition we have not observed any non-monotone numerical solutions in our experiments and no instabilities have occurred.

Numerical experiments
In this section we apply our newly developed mass transport method to several models arising in biomedical applications that bring along merging, emerging, and traveling concentrations phenomena.In particular, we consider the classical Keller-Segel model both elliptic and parabolic.We study also a simple as well as a detailed cancer invasion model.The latter takes the role of the serine protease urokinase-type plasminogen activator into account.The numerical study of such systems constitutes a challenge due to the complex behavior that the solutions exhibit.Numerical experiments presented below demonstrate the robustness and reliability of our newly developed mass transport finite element method.

4.1.
A parabolic-elliptic Keller-Segel model with logistic growth.In the first test case we consider the modified KS system from [7] with added logistic growth which reads Note that the adaptation of system (4.8) to 2D with µ = 0 is equivalent to the simplified Keller-Segel system from [16], where the chemo-attractant c is determined by a Poisson equation.The logistic term accounts for additional cell growth that is locally limited by resources and space.Global existence of solutions to the parabolic-parabolic model with logistic growth in 2D was shown in [25].Except for the logistic source term this model has been numerically investigated by the mass transport scheme employing only inverse distributions in [3].
Since the chemo-attractant density c is given by a convolution term, we do not need to use a finite element approximation.Instead we proceed as in [3] and expand the diffusion taxis operator by During the computation of (4.2a) we control the convergence of the Newton method by comparing subsequent iterates.If the iteration fails to converge, we abort the computation assuming blowup of the numerical solution.The time increments in these experiments have been adapted such, that the Jacobian of (4.2a) that occurs in the Newton iteration is strictly diagonally dominant.The second operator S in this setting accounts only for the logistic growth term.For the numerical simulations we use a grid with only M = 50 points.
We consider a numerical experiment with the parameters D ρ = 1, χ = 2.5π and the initial datum given by This experiment has been studied in the case µ = 0 in [3], where blowup in final time around t = 0.33 has been obtained numerically.We confirm the same phenomenon using the splitting method, see Figure 1.The blowup was indicated by the method during the computation.
When conducting the experiment with altered µ = 0.2, no blowup occurs, as can be seen in Figure 2. The aggregation stops and reverses since the logistic term attracts the cell concentration to a lower density.The total mass of the cells decreases after the aggregation stops and increases again after around t = 1.5.No blowup could be observed even for later times, instead the numerical solution seems to converge to a stationary state.The CFL condition given by (3.4) has caused an increase of the time increment over the computation time.
4.2.Nonlinear diffusion and chemotaxis models.Our method can also resolve models that include generalized diffusion and migration terms as we will demonstrate in this section.To this end we consider at first the model with exponent γ > 1.In the case χ = 0 the first equation in (4.4) is known as the porous media equation modeling the gas flow through a porous interface.We refer to [34] for an introduction to the subject.
Similar as in (4.2), the scheme to resolve (4.4) corresponds to (2.3a)-(2.3b)where the chemo-attractant gradient is computed as We have tested our scheme using again the initial condition (4.3) and the chemo-sensitivity χ = 2.5π. Figure 3 exhibits the results from the numerical simulation for the exponents γ = 2 and γ = 1.5.In both cases the nonlinear diffusion prevents the blowup that would occur for γ = 1 and the numerical solution converges to a stationary state.
Another model that we consider here has been proposed in [26].In this work the authors endowed the Keller-Segel model with a volume filling mechanism.For this purpose they reconsidered the derivation of the model from a random walk and added a function q(ρ) describing the probability that a cell finds sufficient space to jump to a particular position.We adopt here the probability function q(ρ) = 1 − ρ γ that models the volume filling together with enhanced diffusion for γ > 1 and reduced diffusion for γ < 1 [26].Independent of the choice of γ > 1, the model does not allow for cell migration to a position where the maximal density ρ = 1 has already been reached.The corresponding model for the cell density includes nonlinear diffusion and advection terms and reads For the numerical experiments with the volume filling model (4.6) we have adapted the update steps (2.3a) and (2.3b) in the diffusion-advection operator by where we set 2 ∆w Vj+1(t)−Vj−1(t) = 0 for j = 1 and j = M to account for the boundary conditions.In Figure 4 we present simulation results with the parabolic-elliptic model (4.4) where we have replaced the original evolution equation of the cell density by the volume filling approach (4.6).Again we have used the initial condition (4.3) and the chemo-sensitivity parameter χ = 2.5π.We exhibit the numerical results for γ = 2 and γ = 0.5.The computed cell densities do not exceed a density of one in both cases and no blowup occurs.Instead the cells evolve quickly to a bounded spatial profile from which they slowly diffuse afterwards.The parameter choice γ = 2 leads to a larger maximal cell density throughout the computation when compared to the case γ = 0.5.

4.3.
The parabolic-parabolic Keller-Segel model.In this section we apply our scheme to the well known parabolic-parabolic Keller-Segel system which reads As opposed to (4.1) this system features an additional parabolic equation to be treated by the splitting method.In order to exhibit the phenomena that the scheme can resolve, we consider two test cases with distinct initial chemo-attractant densities that control the cell movement.In both tests we adopt the initial datum (4.3) for the inverse distribution V .Figure 5 presents the cell dynamics, showing their movement to the right side of the domain.As the cells produce the chemical with density c, a negative gradient is created that leads to an aggregation of the cells which counteracts the movement.We point out that both the migration and the growth are well resolved by the splitting scheme.to a splitting of the initial concentration into two peaks.The discretization grid for the cells on the density level concentrates its grid points on the locations of both peaks and adapts to the solution over time.
In the setting of the present experiment we study the convergence of the introduced splitting scheme experimentally.For a fixed instance in time and given M , let V h i , i = 1, . . ., M − 1 denote a numerical solution corresponding to the mesh discretization parameter M .Then we define the approximate L 1 finite difference error by where we have used a numerical solution on a finer mesh with 2M points, V h/2 j , j = 1, . . ., 2M − 1, as the reference solution.The experimental order of convergence (EOC) for the discretization error in V h can now be defined by for any even integer M .Similarly, we define the EOC for the cell densities on their non-uniform mesh.To this end let ρ h i , i = 1, . . ., M denote the finite volume representation corresponding to V h i , i = 1, . . ., M − 1 and let E ρ (M ) denote the discrete L 1 error using as reference ρ ref i , i = 1, . . ., 2M the finite volume representation of V h/2 j , j = 1, . . ., 2M − 1.This L 1 error is computed by projecting the reference solution to the coarser non-uniform grid corresponding to the cell densities ρ h i , i = 1, . . ., M .Then we define for even integers M   1. Mesh convergence in the peak splitting experiment up to T = 0.01 with respect to the discretization parameter M .We have adapted the number of points on the Finite Element mesh by N = M .The EOCs approach two in the inverse distributions and in the corresponding densities.Yet, for large M the EOC drops which is probably due to limitations by the Finite Element mesh.
analog to (4.10) In Table 1 we exhibit the errors and the EOCs computed at the final time T = 0.01 and with constant time increment ∆t = 10 −4 when doubling the mesh resolution on the mass space mesh iteratively.We have coupled the resolution of the Finite Element mesh to the number of points for the inverse distribution by using N = M .We can clearly see that the method converges as the mesh size is refined.The EOCs indicate a convergence order of two in both, the inverse distributions and the densities.However, we see that the EOC decreases as the grids becomes very fine.We suppose that this is caused by the finite element mesh that is only uniformly but not locally refined: as M increases the number of mesh cells on the non-uniform finite volume mesh for the cell densities aggregates around the positions of the peaks.Throughout the computation the finite element solution c h must in turn be interpolated in many points in a small physical area which leads to a loss of accuracy as the number max i |{j : 1 demonstrates that the method has provided accurate numerical results using only a few mesh points.

A cancer invasion system.
In this test case we address a model of cancer invasion of the extracellular matrix (ECM), the first step in cancer metastasis.The macroscopic modeling of this process commonly uses an Keller-Segel approach that models the densities of the cancer cells, the concentration of the extracellular fibers on which cancer cells adhere and move, and the density of an enzyme of the matrix metallopreteinases (MMPs) family that is produced by the cancer cells and is responsible for the degradation of the ECM.
There is a wide variety of cancer invasion models in the literature, see e.g.[11,31,30,32,17].In order to test our scheme we employ a simple test case based on the pioneering model [2] augmented with a proliferation term in the cancer cell density equation which reads In this model the cancer cells with the density ρ move using their motility apparatus with a preferred direction towards higher concentrations of the ECM with concentration denoted by v.This is the haptotaxis phenomenon.Being a network in a static equilibrium the ECM does not translocate.The MMPs however, whose density we denote by m, diffuse freely in the extracellular environment.Additionally, the cancer cells proliferate towards a preferred density ρ = 1 and they produce MMPs with a constant rate.The MMPs attach to the ECM which they dissolve upon contact.
For a numerical experiment with this model we consider the computational domain (0, 1) with initial conditions ∆t/∆t ref error E V ∆t EOC t 0.1/ 0.05 2.244e-04 0.05/ 0.025 4.728e-05 2.247 0.025 / 0.0125 1.107e-05 2.094 0.0125 / 0.00625 2.766e-06 2.001 Table 2. Experimental convergence in time (left) and in space (right) in the numerical experiment with system (4.11) at T = 1.In all computations we have set N = M .The EOCs suggest a convergence of second order in time and space.
In our method we discretize both the ECM density v and the MMP concentration m on the same finite element basis.The corresponding approximations are updated in the reaction-diffusion operator of the splitting method.The interpolations are only needed with respect to the ECM density v.We resolve the migration of the cancer cells in transformed variables with the advection-diffusion operator and the cell proliferation in original variables with the reaction-diffusion operator.
The considered numerical experiment simulates the propagation of cancer cells into the ECM on the right side of the computational domain.To account for the corresponding temporal expansion of the support of the cancer cell density c we have adapted the treatment of the right boundary.In more details, we have neglected the discrete cancer cell density entry adjacent to the right boundary in the proliferation update (2.5), i.e. S ∆t ρ M (t) = ρ M (t).Though we have not excluded the corresponding boundary entry in the cumulative function, V M−1 , from the diffusion and haptotaxis updates of the scheme.We present the according numerical results in Figure 7. Apart from the propagation of the cells into the tissue, we observe a build up of cancer cells at the leading front of the tumor.Degradation of the tissue and MMP production are also visible.Throughout the computation the not invaded part of the tissue is resolved by a single grid cell in the cancer cell density.
In this experiment we have also studied the convergence of the scheme experimentally.Along with the errors in space, we have also computed the errors in time by the formula where V h, ∆t i , i = 1, . . ., M − 1 denotes a numerical solution computed on M mesh points with constant time increment ∆t.For the computation of the temporal errors we have considered a fine spatial resolution with M = N = 600 mesh cells.The corresponding EOC is given by ).The spatial errors and EOCs are computed according to (4.9) and (4.10) with constant time increment ∆t = 2 × 10 −4 and coupled N = M .Both, temporal and spatial errors have been computed at the final time T = 1.
In Table 2 we present the computed errors and EOCs in the invasion experiment.We see that the method converges as either the mesh size or the time increment is refined.The EOCs in time and space range around two which confirms our expected second order.As in the "peak splitting" experiment, the EOC decreases slowly as the mesh is refined to very high resolutions.We point out that previous numerical tests which did not employ our proposed boundary treatment have yield only a spatial EOC of one.4.5.The uPA model.In the last series of experiments we apply our scheme to a detailed tumor invasion system derived in [11].This model focuses on the enzymatic urokinase plasminogen activator (uPA) system which is known to play an essential role in the context of cancer progression and metastasis.The uPA is an extracellular serine protease which is responsible for the activation of the protease plasmin.This activation occurs mainly if uPA is bound to its uPA receptors (uPAR) on the cancer cell membrane.The receptor bound uPA enhances the affinity of uPAR to the ECM constituent vitronectin [35] and integrins.Thus, the uPA/uPAR-complex regulates indirectly also the vitronectin-integrin interactions.Both proteases plasmin and uPA catalyze the degradation of vitronectin and other ECM components.Another actor in the system is the plasminogen activator inhibitor type 1 (PAI-1) which is produced by the tumor cells and limits the activation of plasmin to prevent tissue damage and to maintain homeostasis.
The considered model complements the system (4.11) by chemotactic movement of the cells due to uPA and PAI-1, remodeling of the ECM modeled by a logistic term and the dynamics of the uPA system modeled in terms of mass-action kinetics.We refer to [11] for more details.The full model reads where the cancer cell concentration is represented by ρ, the ECM by the density of its constituent vitronectin v, and uPA, PAI-1, and plasmin densities are denoted by u, p, and m.We assume non-negative initial data.
We consider a numerical experiment that we have studied in [22] by a Finite Volume method.It employs the parameter values from [1] given by As done to treat the model (4.11) we use a single finite element basis to discretize the concentrations of the ECM, the uPA, the PAI-1, and the plasmin.The cubic spline in the advection-diffusion operator interpolates the linear combination χ v v + χ u u + χ p p. Similar as in the models (4.8) and (4.11) the scheme approximates the cell proliferation in Eulerian coordinates but diffusion and advection of the cancer cells in transformed variables.We have used the same boundary treatment as in Section 4.4.
In Figure 8 we present the simulation results obtained by our scheme with mesh parameters M = N = 400.The method is capable to approximate accurately the dynamics that we have obtained in [22] including the emergence and movement of multiple steep peeks.The present simulation clearly demonstrates the robustness of the newly developed scheme to simulate complex taxis-diffusion systems arising in cell biology.
To investigate the dynamics of such a cancer invasion system in the case that the cell migration is restricted by the occupied extracellular space we have endowed the model (4.14) with the volume filling approach (4.6).We have used M = N = 400 grid points on both meshes in the numerical simulation.
In more details we have replaced the evolution equation for the tumor cell density in (4.14) by and resolved the same numerical experiment as above.To this end the scheme has been adapted in a similar way as in (4.7).To study how the new method compares in efficiency to more conventional numerical methods we consider again the above experiment without volume filling.For the comparison we consider the Finite Volume/Finite Difference from [22] for both uniform and adaptive meshes.In particular we have chosen a second order method with implicit-explicit Strang operator splitting.For the adaptive mesh refinement (AMR) method we have chosen the gradient monitor function to determine the mesh-cells to be either refined or coarsened * .* In more details we have used the refinement and coarsening threshold values θ ref = 10, θcoars = 2.5, a single refinement and coarsening operation per time step n ref = ncoars = 1 and a maximal refinement level of lmax = 2, cf [22]. .The new MTFE method seems to be most efficient in terms of error per CPU time, its relation between the error and the average number of cells is similar as in the FVFD scheme.
For brevity, we will refer to the adaptive method as AMR and to the uniform method as FVFD.The new mass-transport/finite element method will be denoted by MTFE.
For our comparism we consider the set S = {40, 80, 160, 320, 640, 1280} and run the MTFE method for M ∈ S, the FVFD method for N = 6k for any k ∈ S, and the AMR method for N 0 ∈ S with N 0 denoting the number of cells on the lowest level.We couple the two meshes in the MTFE scheme by setting N = M .We do not consider finer resolutions due to restrictions by the uniform reference solution in the error computations of solutions obtained by the MTFE scheme.For comparison reasons we let N denote the average number of cells in the AMR method.In addition, all three methods employ the same Courant number CF T = 0.49 and all numerical solutions are computed on the domain Ω = (0, 5).
We compute the numerical solutions of the considered experiment at the time instance t = 23 that features two steep peaks in the cancer cell concentration.In this process we measure the CPU time that is needed for the corresponding simulations and compute the error of the approximation at the final time.For the error computation we have used a reference solution that employs a uniform mesh with cell size h = 1.25 × 10 −5 in the relevant part of the domain † .The discrete L 1 error is then computed with respect to the densities using a suitable projection of the reference solution.Note that the following test results are dependent on our (non-reference) implementation of the numerical methods.
We show the results of our comparison in Figure 11.Here we present the relation between the error and the computation time and the relation between the error and the average number of cells for all three methods.We see that for all tested methods the error decreases as either the cell number or the CPU time increases.Figure 11 (left) exhibits an advantage of the new MTFE method over the other schemes in efficiency for most of the conducted simulations.This can be seen as the MTFE method achieves in most cases lower errors than the FVFD or the AMR scheme using the same CPU time.As the runtime increases the MTFE method approaches the efficiency of the AMR method with the new method being at a slight advantage over the mesh refinement method.Clearly, the AMR and the MTFE scheme both outperform the FVFD method for sufficiently large CPU times.
Figure 11 (right) shows that the AMR method achieves the lowest errors when compared with simulations by the FVFD and MTFE scheme employing the same average number of cells.The error of the MTFE scheme has a similar dependence on the number of cells as the error of the FVFD scheme.We conjecture thus that the better efficiency of the MTFE scheme in terms of CPU time seen in Figure 11 (left) is probably caused by the CFL condition in the MTFE scheme allowing for larger time steps compared to the FVFD method.

Conclusion
In this paper we have proposed a new splitting scheme for one-dimensional reaction-taxis-diffusion systems related to the Keller-Segel system.The solutions of these systems are well known for having concentrated and diffusive regions simultaneously.In addition, traveling waves and merging phenomena typically occur.
Our splitting has separated a part of the model which is mass conservative in the cell density from the rest of the system.The latter has been approximated by a classical linear finite element method, whereas the approximation of the conservative part has been based on the mass transport strategy.More precisely, we have first transformed the cell density to the corresponding pseudo-inverse cumulative distribution.Then we have discretized the transformed system by the finite difference method and used a cubic spline to account for the chemo-attractant whose evolution is described in the rest subsystem.The splitting method is described in Section 2. In Lemma 3.1 we have studied the stability of the explicit mass transport method for the conservative part in which we allowed for general nonlinear diffusion.The obtained result has been used to derive a time-step restriction for our scheme.
In Section 4 we have presented a series of numerical experiments demonstrating the robustness and reliability of the scheme.In particular, we have used the new method to resolve the Keller-Segel model in the parabolicelliptic and in the parabolic-parabolic form numerically.We have applied our scheme also to augmentations of these systems by reaction terms, nonlinear diffusion and a volume filling approach.The method has resolved the movement, splitting and aggregation phenomena accurately.We have verified the mesh convergence of the scheme in both time and space in an application a simple tumor invasion system in Section 4.4.The obtained experimental order of convergence has ranged around two spatially and temporally.Moreover, we have applied the scheme to the uPA-tumor invasion model from [11] in Section 4.5.The proposed hybrid mass transport finite element scheme has been capable to resolve its complex dynamics featuring multiple peaks in the cancer cell concentration without using a fine spatial discretization.By the help of our new method we could also study a combination of the uPA model with the volume filling approach from [26].In addition, we have compared the efficiency of the hybrid mass transport finite element method with a finite volume scheme with adaptive mesh refinement from [22].The hybrid mass transport finite element method has not only outperformed the uniform finite volume scheme but it has also delivered slightly better results than the finite volume scheme equipped with adaptive mesh refinement.
the computation of the linear systems (2.5b) and (2.5d).The integrals of the form b

Figure 1 .
Figure 1.Numerical results (cell concentration and inverse cumulative function) for the parabolic elliptic Keller-Segel model, experiment 4.3.1 in [3].The cell concentration blows up.The numerical cell concentration has attained a maximum of approximately 176.

Figure 2 .
Figure 2. Numerical results (cell concentration, inverse cumulative function, and mass) for the parabolic elliptic Keller-Segel model with added logistic growth (4.1).The additional reaction term has prevented blowup.
Peak splitting.In the next test we use the parameters D ρ = D c = 0.1, χ = 5, α = β = 1 and the computational domain (a, b) with boundaries chosen as in the "peak movement" experiment.The initial chemo-attractant density though is replaced by the function c 0 (x) = 1 − e −20 x 2 , x ∈ (a, b).

Figure 6 Figure 4 .
Figure6shows the computational results for the discretization parameters M = 90 on the mass space mesh and N = 450 on the Finite Element mesh.The cells move out of the center of the domain on which the most part of the attracting chemical is already consumed.The symmetrical movement to both sides leads

Figure 5 .
Figure 5. Numerical results (cell concentration, inverse cumulative function, and chemo-attractant density in space and time) for the parabolic-parabolic KS model (4.8) in the "peak movement" experiment.The movement and aggregation is accurately resolved using M = 45 grid points.

Figure 6 .
Figure 6.Numerical results (cell concentration, inverse cumulative function, and chemo-attractant density in space and time) for the parabolic-parabolic KS model (4.8) in the "peak splitting" experiment.At the final time we show the approximated cell averages of the density ρ at their respective position on the grid (top).The grid for the cell density adapts to the two splitting peaks.

Figure 7 .
Figure 7. Numerical results (cell concentration, inverse cumulative function, tissue and MMP density) for the cancer invasion model (4.11) using only N = M = 45 mesh points.At the final time we show the approximated cell averages of the tumor density ρ at their respective position on the grid (top left).A high concentration of tumor cells emerges and invades the tissue.The grid for the tumor density omits the part of the tissue that is not yet invaded.

Figure 8 .
Figure 8. Numerical results (cancer cell concentration, inverse cumulative function, ECM, uPA, PAI-1 and plasmin density in space and time) in an numerical experiment with the model (4.14) computed by the new scheme.The dynamics, particularly the steep peaks in the cancer cell density, are well resolved by the scheme.We have used M = N = 400 grid points on both meshes in the numerical simulation.

Figure 9 .Figure 10 .
Figure 9. Numerical results (cancer cell concentration, inverse cumulative function, ECM, uPA, PAI-1 and plasmin density in space and time) in the model (4.14) with volume filling by (4.15) and exponent γ = 2.We have used M = N = 400 grid points on both meshes in the numerical simulation.

Figure 11 .
Figure 11.Relation between the CPU time and the error (left) and between the (average) number of cells and the error (right) for the FVFD, AMR, and MTFE scheme in log-log scale in a numerical experiment with the uPA model(4.14).The new MTFE method seems to be most efficient in terms of error per CPU time, its relation between the error and the average number of cells is similar as in the FVFD scheme.