Global existence of classical solutions and numerical simulations of a cancer invasion model

In this paper, we study a cancer invasion model both theoretically and numerically. The model is a nonstationary, nonlinear system of three coupled partial differential equations modeling the motion of cancer cells, degradation of the extracellular matrix, and certain enzymes. We first establish existence of global classical solutions in both two- and three-dimensional bounded domains, despite the lack of diffusion of the matrix-degrading enzymes and corresponding regularizing effects in the analytical treatment. Next, we give a weak formulation and apply finite differences in time and a Galerkin finite element scheme for spatial discretization. The overall algorithm is based on a fixed-point iteration scheme. In order to substantiate our theory and numerical framework, several numerical simulations are carried out in two and three spatial dimensions.


Introduction
The model One of the defining characteristics of a malignant tumour is its capability to invade adjacent tissues [21]; accordingly the mathematical literature directed at understanding underlying mechanisms is vast (see e.g. the surveys [29,36]).
In this paper we focus on the following variant of a cancer invasion model developed by Perumpanani et al. [34] for the malignant invasion of tumours and investigate in Ω × (0, ∞), in Ω × (0, ∞), on Ω × (0, ∞), ( , , )(·, 0) = ( 0 , 0 , 0 ) in Ω. (1.1) We aim for a rigorous existence proof for global solutions and the development of a reliable numerical scheme with an implementation in a modern open-source finite element library.
In (1.1), the motion of cancer cells (density denoted by ) mainly takes place by means of haptotaxis, i.e. directed motion toward higher concentrations of extracellular matrix (density ), of strength ≥ 0. Motivated by experiments of Aznavoorian et al. [7], who reported only "a minor chemokinetic component" of the cell motion, the original model of [34] does not include a term for random (chemokinetic) cell motility at all. Acknowledging that "minor" does not mean "none at all", we deviate from [34] in this aspect and incorporate this motility term in (1.1) ( ∈ (0, ∞), with the formal limit → ∞ corresponding to the model of [34]). Additional growth of the population of tumour cells is described by a logistic term (with being a positive parameter). The extracellular matrix is degraded upon contact with certain enzymes (proteases, concentration ), which, in turn, are produced where cancer cells and matrix meet and decay over time. The reaction speed of these protein dynamics can be adjusted via the parameter > 0. As many proteases remain bound to the cellular membrane -or are only activated when on the cell surface (cf. the model derivation in [34]), no diffusion for is incorporated in the model. This last point is in contrast to the otherwise similar popular models in the tradition of [4,33] or [9], the latter of which additionally included a chemotactic component of the motion of cancer cells.

Global solvability
In order to construct global classical solutions of (1.1), it is necessary to control the haptotaxis term − ∇ · ( ∇ ) in the first equation and thus in particular to gain information on the spatial derivative of the second solution component. For relatives of (1.1) including a diffusion term ∆ in the third equation, this has already been achieved in [40] and [28] by applying parabolic regularity theory to the equation for , first yielding estimates for the spatial derivative of and then also on ; the results of [28] even cover the long-term asymptotics of solutions. Moreover, the presence of diffusion for the produced quantities has also been made use of to obtain global existence results for different cancer invasion models, see for instance [42].
However, the absence of any spatial regularization in both the second and third equation makes the corresponding analysis much more challenging. Up to now, global classical solutions have only been constructed for a rather limited set of initial data: already in [34], where (1.1) has been proposed for = ∞, it has been shown that the model formally obtained by taking the limit ↘ 0 admits a family of travelling wave solutions. Corresponding results for positive have then been achieved in [31]. Moreover, if = 0, travelling wave solutions may contain shocks [30] and solutions of related systems without a logistic source may even blow up in finite time [35]. In general, the destabilizing effect of taxis terms such as − ∇ · ( ∇ ) may not only make it challenging but even impossible to obtain global existence results for certain problems. We refer to the survey [25] for further discussion regarding the consequences of low regularity in chemotaxis systems.
Despite these challenges, in the first part of the present paper we are able to give an affirmative answer to the question whether (1.1) also possesses global classical solutions for widely arbitrary initial data in the two and three-dimensional setting. Our analytical main result is the following which, moreover, is nonnegative.

Numerical modeling
In the second part of our paper we then analyze the behavior of these solutions numerically with an implementation in the modern open-source finite element library deal.II [5,6]. Related numerical studies in various software libraries using different numerical schemes are briefly described in the following. The traditional method of lines has been widely used for simulations of the cancer invasion process [15,19]. In addition, finite difference methods [23] have been considered and in [10,22], the authors proposed a nonstandard finite difference method which satisfies the positivity-preservation of the solution, that is an important property in the stability of the model. Moreover, the finite volume method [11], spectral element methods [41], algebraically stabilized finite element method [37], the discontinuous Galerkin method [16], combinations of level-set/adaptive finite elements [1,44], and a hybrid finite volume/finite element method [2,8] have also been proposed in the literature for some cancer invasion models and chemotaxis. Finally, we mention that in [38] the authors illustrate their theoretical results for a related cancer model employing discontinuous Galerkin finite elements implemented as well in deal.II.
The main objective in the numerical part is the design of reliable algorithms for (1.1) and their corresponding implementation in deal.II. First, we discretize in time using a -method, which allows for implicit -stable time discretizations. Then, a Galerkin finite element scheme is employed for spatial discretization. The nonlinear discrete system of equations is decoupled by designing a fixed-point algorithm. This algorithm is newly designed and then implemented and debugged in deal.II.
These developments then allow to link our theoretical part and the numerical sections in order to carry out various numerical simulations to complement Theorem 1.1. Specifically, several parameter variations of the proliferation coefficient and the haptotactic coefficient will be studied in two-and three spatial dimensions. These studies are non-trivial due to the nonlinearities and the high sensitivity of (1.1) with respect to such parameter variations.

Plan of the paper
The outline of this paper is as follows. In Section 2, we study the global existence of classical solutions. Next, in Section 3, we introduce the discretization in time and space using finite differences in time and a Galerkin finite element scheme in space. We also describe the solution algorithm. In Section 4, we carry out several numerical simulations demonstrating the properties of our model and the corresponding theoretical results. Therein, we specifically study parameter variations. Finally, our work is summarized in Section 5.

Global existence of classical solutions
As a first step in the proof of Theorem 1.1, we apply two transformations in Subsection 2.1; the first one allows us to get rid of some parameters in (1.1), the second one changes the first equation to a more convenient form. We then employ a fixed point argument to obtain a local existence result for the transformed system in Lemma 2.5.
The proof that these solutions are global in time consists of two key parts, both relying on the fact that the second and third equation in (1.1) at least regularize in time (which allows us to prove Lemma 2.7 and Lemma 2.10). First, in order to prove boundedness in ∞ , the comparison principle allows us to conclude boundedness in small time intervals (cf. Lemma 2.8). We then iteratively apply this bound to obtain the result also for larger times (cf. Lemma 2.9). As to bounds for the spatial derivatives, we secondly apply a testing procedure to derive estimates valid on small time intervals (cf. Lemma 2.11), which then is again complemented by an iteration procedure (cf. Lemma 2.12). Finally, we are able to make use of parabolic regularity theory (inter alia in the form of maximal Sobolev regularity) to conclude in Lemma 2.14 that the solutions exist globally.

Two transformations
We first note that with regards to Theorem 1.1 we may without loss of generality assume = 1 and = 1. for˜∈Ω := 1 √ Ω. By Theorem 1.1, there then exists a global classical solution of (1.1) (with all parameters and initial data replaced by their pendants with tildes) (˜,˜,˜). Then Then ∇ = e − ∇ − e − ∇ , so that and (1.1) (with = = 1) is equivalent to in Ω (2.4) for 0 := 0 e − 0 . Expanding ∇ · ( ∇ ) to ∇ · ∇ + ∆ shows that this transformation allows us to get rid of a term involving ∆ in the first equation at the price of adding several zeroth order terms. In particular as the second equation does not regularize in space, (2.4) turns out to be a more convenient form for the following analysis.

Local existence
In this subsection, we construct maximal classical solutions of (2.4) in Ω × [0, max ) for some max ∈ [0, ∞) by means of a fixed point argument. Moreover, we provide a criterion for when these solutions are global in time (that is, when max = ∞ holds), which then will finally be seen to hold true in Lemma 2.14.
As a preparation, we first collect results on (Hölder) continuous dependency of solutions to ODEs on the data.

for every
∈ Ω the existence and uniqueness of a solution By an ODE comparison argument, nonnegativity of follows from that of 0 , and nonnegativity of from that of 0 , and . Therefore, 0 ≤ ( , ) ≤ 0 ( ) for all ∈ Ω and ∈ (0, max ( )) due to the sign of = − . Consequently, }︂ for all ∈ Ω, ∈ [0, max ( )), which also shows that max ( ) = for all ∈ Ω. The remainder of part (a) follows from an application of Lemma 2.1 with 1 = 2 = 0.
(b) We apply Lemma 2.1 to = (∇ , ∇ ), the solution of noting that sufficient regularity is given by the assumption on and part (a). (c) Follows from Lemma 2.2.
Both as an ingredient to the proof of Lemma 2.5 below and also for its own interest, we note that classical solutions of (2.4) are unique.
Making use of Schauder's fixed point theorem and applying Lemmas 2.1-2.4, we now obtain a local existence result for the system (2.4).
By fixing initial data as in (2.8), we henceforth implicitly also fix the unique classical solution ( , , ) of (2.4) given by Lemma 2.5 and denote its maximal existence time by max .

∞ bounds
The results in the previous subsection show that Theorem 1.1 follows once we have shown that for the solutions constructed in Lemma 2.5 the second alternative in (2.9) cannot hold. That is, we need to derive sufficiently strong a priori estimates. In this subsection, we begin with bounds in ∞ .
For the second solution component, such a bound directly follows from the comparison principle.
As a preparation for obtaining ∞ estimates also for the other two solution components, we note that the time regularization in the third equation in (2.4) implies that we can bound by a quantity including an arbitrarily small contribution of the ∞ norm of -at least if we are willing to shrink the time interval on which this estimates holds accordingly.
We now turn our attention to ∞ estimates of . The fact that the transformed quantity fulfills an equation whose first-and second-order terms reduce to e − ∇ · 3)), an expression without any explicit ∇ , opens the door for certain testing procedures. In related works, these have been used to first derive boundedness in and then, after an iteration argument, also in ∞ (see for instance [42,Prop. 5.1], [40,Lemma 3.5] and [39,Lemma 3.10]). However, here we are able to employ a slightly faster method: another advantage of the transformation = e − is that sufficiently large constant functions are supersolutions of the equation for in (2.4), at least as long both and are bounded. Therefore, the previous two lemmata allow us to prove boundedness for on small timescales.
We also emphasize that the following proof crucially relies on the presence of a logistic source in the first equation, that is, on positivity of . (The same would be true for testing procedures similar to those performed in the works referenced above.) In fact, this is essentially the only place where we directly make use of the assumption > 0. Proof. We fix ∈ (0, ⋆ ] ∩ (0, max ). By Lemmas 2.6 and 2.7, we may estimate in Ω × (0, ). Therefore, the comparison principle, applied with the constant supersolution and thus ‖ ‖ ∞ (Ω×(0, )) ≤ 2 max Since ‖ ‖ ∞ (Ω×(0, )) ≤ + ‖ ‖ ∞ (Ω×(0, )) by Lemma 2.7, this implies the statement for : Next, iteratively applying Lemma 2.8 allows us to derive boundedness for all solution components on all bounded time intervals. We note that a prerequisite for such an iteration procedure is that the time ⋆ given by Lemma 2.7 (and to a lesser extent also the constant given by Lemma 2.8) only depends on the data in a manageable way. This justifies why we have kept track of the dependencies of the constants in the previous lemmata. Lemma 2.9. Suppose (2.7) and let > 0. Then there are 1 , 2 > 0 with the following property: for all > 0 and all ( 0 , 0 , 0 ) satisfying (2.8) and for all ∈ (0, max ). (2.14) Proof. We set (·, ) = 0 and (·, ) = 0 for < 0, and let ⋆ > 0 and > 1 be as given by Lemmas 2.7 and 2.8, respectively. Moreover, setting for ∈ N 0 and applying Lemma 2.7 (which is applicable for the same by Lemma 2.6) and Lemma 2.8 to initial data ( , , )(·, ( − 1) ⋆ ), we can estimate for all ∈ N with ( − 1) ⋆ < max . However, the same estimates hold trivially also in the case of ( − 1) ⋆ ≥ max , as then = ∅. Thus, A straightforward induction then yields ≤ (2 + 1) 0 ≤ e ln(2 +1) (2 + 1) for all ∈ N. If ∈ for some ∈ N 0 and thus ( − 1) ⋆ < , that is < ⋆ + 1, then , this implies (2.14) for 1 = ln(2 +1) ⋆ and 2 = 2(2 + 1).

Gradient bounds in 4
While the ∞ estimates proven in the previous subsection surely form an important step towards proving global existence, the extensibility criterion (2.9) also requires boundedness of the gradients -which will be the topic of the present and the following subsection.

Hölder estimates for the gradients
Lemmas 2.6, 2.9 and 2.12 provide several bounds for the right-hand side of the first equation in (2.4), which allow us to make use of parabolic regularity theory to iteratively improve our bounds. In particular, we adapt the techniques developed in [39, pp. 791-792], where only planar domains are considered, to the three-dimensional setting.
As it is used multiple times in the proof of Lemma 2.14 below, we first state the following consequence of maximal Sobolev regularity results.
This in turn renders [ , contradicting the extensibility criterion in Lemma 2.5. Thus our assumption that max is finite must be false.

Proof of Theorem 1.1
The proof of Theorem 1.1 has now been reduced to referencing some of the lemmata above.
Proof of Theorem 1.1. Lemma 2.5 asserts the local existence of a unique maximal classical solution of (2.4), which is global in time by Lemma 2.14. Therefore, the statement follows by transforming back to the original variables; that is, first setting ( , ) := ( , )e ( , ) for ∈ Ω and ∈ [0, ∞) and then applying the transformation in (2.1).

Weak formulation, discretization and numerical solution
In this section, we address the numerical realization of (1.1). To this end, we first derive a weak formulation and then apply the Rothe method, namely, first temporal discretization using finite differences, and afterward spatial discretization based on a Galerkin finite element scheme. Due to the highly nonlinear behavior, we then propose and implement a fixed-point algorithm to solve all three equations sequentially. Similar algorithms and implementations are available in the deal.II library [5,6], and we have former experiences in solving highly nonlinear coupled PDE systems, e.g., [43], but the algorithmic design, implementation and code verification of (1.1) in deal.II is novel to the best of our knowledge.

Temporal discretization and fixed point scheme
and and for = 1, 2, . . . , * , where * is the iteration index where some stopping criterion is met, and for = 0, 1, . . . , − 1. For details on the specific steps and stopping tolerances we refer the reader to Section 3.4.
Denoting by 1 (ˆ) the space of polynomials on the reference cellˆ(square in two dimensions and cube in three dimensions) which are linear in each variable, the shape functions from 1 ( ) are obtained using 1 (ˆ) transformations of functions in 1 (ˆ) onto , so-called isoparametric finite elements. We refer the reader to the classical textbook [12] for more details. Moreover, we denote by (·, ·) the scalar product in 2 (Ω). The discrete solutions +1 ℎ , +1 ℎ and +1 ℎ are written as linear combinations of standard basis functions of ℎ : where denotes the number of spatial degrees of freedom, i.e., dim ℎ = . The fully discrete system then reads as follows: and and ∈ ℎ , respectively, analogously as in (3.2). Each linear system is solved with a sparse direct solver.
Remark 3.2. The system of algebraic equations of each equation at each step is solved using the sparse direct solver UMFPACK [14]. Remark 3.3. We notice that the relaxation parameter can also be obtained via a backtracking procedure starting with = 1 and constructing a sequence with → 0 for → ∞ until convergence is achieved.

Numerical simulations
In order to illuminate the evolution of solutions and show their qualitative behavior, beyond the mere existence assertion of Theorem 1.1, in this section, we perform several numerical simulations in two and three spatial dimensions. The main objective are investigations of the influence of variations in the proliferation coefficient and the haptotactic coefficient , whose size did not matter for Theorem 1.1. The specific values are chosen for illustrative purposes, not due to their biological relevance. For a discussion of realistic diffusion and taxis coefficients of tumor cells see, for instance, [3].  Table 1.

Simulations for different proliferation coefficients
First, we study the influence of the cancer cell proliferation coefficient on the cancer invasion for = 10 −10 , 0.5, 1.0 with small haptotactic rate = 0.01. We notice that our theory in Section 2 requires > 0 and for this reason we made the previous choice = 10 −10 . Numerically we are interested in a value being close to zero in order to study the behavior of the cancer invasion model. Proliferation shows the ability of a cancer cell to copy its DNA and divide into 2 cells, therefore an increase in the proliferation rate of tumours causes an accelerated invasion of cancer cells into connective tissues domain. In all the computations we use −1 = 0.1, = 0.2.
The results obtained with the standard Galerkin discretization of the system (1.1) are displayed in Figures 1-6, at time instances 5, 15, 25 and 35. The snapshots of cancer cell invasion, connective tissue and protease are plotted in Figures 1, 3 and 5. We start with = 10 −10 , that is, almost no growth in the cancer cell density. As we can see from Figure 1, there is no growth in the cancer during the initial stage, and despite a small amount of concentration at the initial period, the cancer cell density and also protease (which is produced by cancer cells upon contact with connective tissues) are decreased and spread slowly due to diffusion effect and the invasion does not continue after time = 15. Now, let us consider = 0.5. As we can see from Figure 3, in this case, an increase of the concentration of cancer cells becomes visible, and it continues during the time. The cancer invasion gradually increases and degrades nearly half of the connective tissue by the time = 25. For = 1.0, Figure 5 shows the growth effect. Due to high proliferation rate, cancer cells produce more protease, which helps them to invade the connective tissues area rapidly. In particular, cancer cells complete invasion in three-quarters of the connective tissue domain at = 25 when = 1.0 is used. The snapshots of cancer cell invasion for different values of proliferation rate are given in Figures 2, 4 and 6. As explained, by increasing the value of , cancer cells increase and the invasion happens more rapidly for all the considered time intervals.

Identical proliferation and haptotactic coefficients
In this subsection, we consider the case when the proliferation rate is equal to haptotactic coefficient, i.e., = = 1, and all other parameters are the same as in the previous subsections. As it can be seen from Figures 13 and 14, due to the proliferation rate, the concentration of cancer growths quickly even from the beginning resulting from a high amount of haptotaxis, therefore the tumour migrates rapidly inside the domain and degrades the connective tissue in a much shorter amount of time.

Three dimensional simulations
In this final subsection, we perform numerical simulations in three spatial dimensions to consider some more realistic movement. Here, the experiments are performed on a mesh with 32 768 hexahedral elements covering the domain Ω. Figures 15 and 16 show the snapshots of cancer cells and connective tissues for growth rate = 1 and haptotactic coefficient = 1. Further, we use the parameters −1 = 0.1 and = 0.2. As it can be seen, at = 5 the connective tissue covers the entire domain and only a small amount of cancer cells exists at the corner, by the time cancer cells growth and invade the domain of connective tissue quickly and by = 35 almost all the domain is occupied by cancer cells.

Conclusions
In this paper, we established theoretical proofs, numerical algorithms, implementations and numerical simulations for a cancer invasion model. In our theoretical part, existence of global classical solutions in both twoand three-dimensional bounded domains was established. In the proofs, we employed the fact that the second and third equation in (1.1) at least regularize in time. For showing boundedness in ∞ , the comparison principle allowed us to conclude boundedness in small time intervals, which then was iteratively applied to obtain the result also for larger times. For the spatial derivatives, we secondly applied a testing procedure for deriving estimates valid on small time intervals, again followed by an iteration procedure. Parabolic regularity theory yielded global existence of the solutions.
The numerical stability of the system heavily depends on the haptotactic coefficient . By fixing proliferation rate and varying the one can make either the diffusion or transport of the cells dominant. The later usually gives rise to spurious oscillations or numerical blow up in the system. In order to study such properties, (1.1)  was discretized using finite differences in time and Galerkin finite elements in space. A fixed-point scheme was designed to decouple the three equations, yielding a robust nonlinear procedure. These developments and their implementation allowed us to study numerically variations in and in two and three spatial dimensions and to illustrate our theoretical results.
Compared to other models, the system in this article did not feature any spatial regularizing effects in the third (or second) equation. This was based on the modelling in [34], where it was argued that there should be no diffusion term for the protease equation. Key challenges both in the analytical and numerical part are precisely caused by this biologically motivated choice. Related works treating systems including a diffusion term also for the third equation crucially make use of the corresponding smoothing effects -a direct adaptation of their methods would evidently have been insufficient for the system at hand.
As to future work, we notice that higher parameter variations resulting into convection-dominated regimes, require the design and implementation of stabilization methods such as streamline upwind Petrov-Galerkin stabilizing formulations or algebraic flux corrected transport. This would introduce additional terms in the equation (3.3) of our nonlinear fixed-point scheme. In case of an algebraic stabilization, an additional nonlinearity would be possibly created since it involves limiters that often depend on the unknown discrete solution. Nevertheless, this nonlinearity can be treated in the framework of the considered fixed-point iterations so that it does not increase the computational cost significantly.