From a cell model with active motion to a Hele–Shaw-like system: a numerical approach

FROM A CELL MODEL WITH ACTIVE MOTION TO A HELE-SHAW-LIKE SYSTEM. A NUMERICAL APPROACH FRANCISCO GUILL´ EN-GONZ´ ALEZ†AND JUAN VICENTE GUTI´ ERREZ-SANTACREU‡ Abstract. In this paper we deal with the numerical solution of a Hele–Shaw-like system via a cell model with active motion. Convergence of approximations is established for well-posed initial data. These data are chosen in such a way the time derivate is positive at the initial time. The numerical method is constructed by means of a finite element procedure together with the use of a closed-nodal integration. This gives rise to an algorithm which preserves positivity whenever a right-angled triangulation is considered. As a result, uniform-in-time a priori estimates are proven which allows us to pass to limit towards a solution to the Hele–Shaw problem. 2010 Mathematics Subject Classification. 92C50, 35B25, 35K55, 35Q92, 35R35, 76D27. Keywords. Finite-element approximation; nonlinear diffusion; free boundary problems; Hele- Shaw flows. Contents 1. Introduction 1 1.1. The models 1 1.2. Notation 3 1.3. Outline 3 2. Spatial discretization 4 2.1. Finite-element approximation 4 2.2. The numerical scheme 7 2.3. Main result 8 3. Proof of Theorem 2.2 8 3.1. A priori energy estimates 8 3.2. Passing to the limit 12 4. An algorithm on unstructured meshes 17 5. Numerical simulation 19 5.1. Temporal integration 19 5.2. Computational experiments 19 References 22 1. Introduction 1.1. The models. Tumour cells are active mechanical systems that are able to produce forces which cause random migration [3,8,14]. This movement is due to rather complicate mechanisms which occur inside cells and give rise to changes in cell shape. Another important mechanism under which cells move is pressure [5,8,13] as a consequence of space competition generated by Date: November 9, 2018. This work was partially supported by Ministerio de Econom´ıa y Competitividad under Spanish grant MTM2015- 69875-P with the participation of FEDER. 1 2 F. GUILL´ EN-GONZ ´ ALEZ AND J. V. GUTI´ ERREZ-SANTACREU cell proliferation itself. In the setting up we take into consideration a very simplified model which incorporates the two spatial effects for describing tumour growth. Let Ω be a connected, open, bounded set of Rd, with d= 2 or 3, and [0, T] a time interval. Consider the cell model with active motion [11] which consists in finding a tumour cell population density n: Ω ×[0, T]→R+satisfying (1) ∂tn−∇·(n∇p(n)) −ν∆n=n G(p(n)) in Ω ×(0, T), subject to the (natural) boundary condition (2) ∇n·n= 0 on ∂Ω×(0, T), with nbeing the outwards unit normal vector on the boundary ∂Ω, and the initial condition (3) n|t=0 =n0in Ω. Here p: [0,+∞)→[0,+∞) is defined by (4) p=p(n) := k k−1nk−1∀n≥0,(k∈N, k ≥2), and G=G(p) is a truncated decreasing function such that there exists Pmax >0 (the homeostatic pressure) with (5) G(0) >0, G(p) = 0 ∀p≥Pmax >0,and G0(p)<0∀p∈(0, Pmax). In the above, Gstands for the decrease in the tumuor cell growth rate when space is limited; the lack of space is governed by the local pressure p, the parameter Pmax is the maximum pressure threshold that tumour cells can exceed before entering a quiescent state, and the parameter ν > 0 represents the effect of including the active (random) motion of cells. It should be noted that the relationship of p(n) given in (4) is invertible for n≥0: (6) n(p) := k−1 kp1/(k−1) ∀p≥0. In this work we assume that {n0 k}k∈Nis a sequence of initial data (3) for (1) such that (7) 0 ≤p(n0 k)≤Pmax in Ω, and that there exists a limit function n0 ∞such that (8) n0 k→n0 ∞in Lp(Ω)-strongly for any p < ∞as k→ ∞. Consequently, defining Nmax(k) := n(Pmax) with n(·) being given in (6), we have (9) 0 ≤n0 k≤Nmax(k) in Ω. from which we infer that there must exist N0>0 such that Nmax(k)≤N0. Under the above assumptions, equation (1) generates a sequence of solutions {nk}k∈Nwhich lead to a solution describing the dynamics of tumour growth as a free-boundary problem. To be more precise, the convergence of the solutions {nk}k∈Nof the active motion cell model problem (1)-(3) towards a weak solution to a Hele–Shaw-like system, as the parameter kgoes to infinity, was proven in [11]. This limit system reads as follows. Find n∞: Ω ×[0, T]→R+and p∞: Ω ×[0, T]→R+such that (10) ∂tn∞−∆p∞−ν∆n∞=n∞G(p∞) in Ω ×(0, T), subject to (11) n∞|t=0 =n0 ∞in Ω, (12) ∇n∞·n= 0 and ∇p∞·n= 0 on ∂Ω×(0, T), APPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 3 jointly to the complementary relation (13) p∞(∆p∞+G(p∞)) = 0 in Ω ×(0, T). The key point in establishing convergence is imposing that ∂tnk(0) ≥0. Moreover, equation (10) is equivalent to solving (14) ∂tn∞−∇·(n∞∇p∞)−ν∆n∞=n∞G(p∞) in Ω ×(0, T). This equivalence will be accomplished due to the equality ∇p∞=n∞∇p∞, which comes from the equalities p∞∇n∞= 0 and p∞n∞=p∞. In this paper, we shall be concerned with the convergence of a finite element scheme, the time variable being continuous, for the active motion cell model problem (1)-(3) towards the Hele-Shaw system (10)-(13) as the space discrete parameter hgoes to zero and kgoes to infinity. 1.2. Notation. We will assume the following notation throughout this paper. Let O ⊂ RM, with M≥1, be a Lebesgue-measurable set and let 1 ≤p≤ ∞. We denote by Lp(O) the space of all Lesbegue-measurable real-valued functions, f:O → R, being pth-summable in Ofor p < ∞ or essentially bounded for p=∞, and by kfkLp(O)its norm. When p= 2, the L2(O) space is a Hilbert space whose inner product is denoted by (·,·). To shorten the notation, the norm k·kL2(Ω) is abbreviated by k·k. Let α= (α1, α2, ..., αM)∈NMbe a multi-index with |α|=α1+α2+... +αM, and let ∂αbe the differential operator such that ∂α=∂ ∂x1α1...∂ ∂xdαM. For m≥0 and 1 ≤p≤ ∞, we define Wm,p(O) to be the Sobolev space of all functions whose mderivatives are in Lp(O), with the norm kfkWm,p(O)= X |α|≤mk∂αfkp Lp(O)  1/p for 1 ≤p < ∞, kfkWm,p(O)= max |α|≤mk∂αfkL∞(Ω),for p=∞, where ∂αis understood in the distributional sense. For p = 2, Wm,2(O) will be denoted by Hm(O). We also consider C∞(O) to be the space of functions continuously differentiable any number of times, and C∞ c(O) to be the subspace of C∞(O) with compact support in O. Spaces of Bochner-measurable functions from a time interval [0, T] to a Banach space Xwill be denoted as Lp(0, T;X) with kfkL2(0,T ;X)=RT 0kf(s)kp Xdsif 1 ≤p < ∞or kfkL∞(0,T,X)= ess sups∈(0,T)kf(s)kX<∞if p=∞. 1.3. Outline. Next we sketch the remaining content of this work. In section 2 we present our finite-element spaces and some preliminary result mainly concerning interpolation operators. Furthermore, we set out our finite element numerical method, where the time variable remains continuous, and the main result of this paper. Next is section 3 which is devoted to demonstrating the main result. Firstly, a discrete maximum principle for finite-element approximations is achieved by assuming a partition of the computational domain being made up of right-angled simplexes, and a priori estimates are also established independent of (h, k) with hbeing the space parameter associated to our finite-element space. As a result, we are able to prove positivity for the time derivative of finite-element approximations. Then better a priori energy estimates lead to obtaining compactness for passing to the limit as (h, k)→(0,+∞). In section 4, we propose a variant of our numerical algorithm for nonobtuse triangulations which keeps with a discrete maximum principle 4 F. GUILL´ EN-GONZ ´ ALEZ AND J. V. GUTI´ ERREZ-SANTACREU and positive for the discrete time but whose convergence is not clear. Finally, in section 4, some numerical experiments are presented for studying the behavior of several parameters. 2. Spatial discretization 2.1. Finite-element approximation. Herein we introduce the hypotheses that will be required along this work. (H1) Let Ω be a bounded domain of Rd(d= 2 or 3) with a polygonal or polyhedral Lipschitzcontinuous boundary. (H2) Let {Th}h>0be a family of shape-regular, quasi-uniform triangulations of Ω made up of right-angled simplexes being triangles in two dimensions and tetrahedra in three dimensions, so that Ω = ∪K∈ThK, where h= maxK∈ThhK, with hKbeing the diameter of K. Further, let Nh={ai}i∈Idenote the set of all the nodes of Th. (H3) Conforming piecewise linear, finite element spaces associated to Thare assumed for approximating H1(Ω). Let P1(K) be the set of linear polynomials on K; the space of continuous, piecewise P1(K) polynomial functions on This then denoted as Nh=nh∈C0(Ω) : nh|K∈ P1(K)∀K∈ Th, whose Lagrange basis is denoted by {ϕa}a∈Nh. We now give some auxiliary results for later use. We begin by an inverse inequality whose proof can be found in [4, Lem. 4.5.3] or [9, Lem. 1.138]. Proposition 2.1. Under hypotheses (H1)–(H3), it follows that, (15) k∇nhkL2(K)≤Cinv h−1 KknhkL2(K)∀K∈ Th,∀nh∈Nh, where Cinv >0is a constant independent of h. Let Ihbe the nodal interpolation operator from C0(Ω) to Nhand consider the discrete inner product (nh, nh)h=ZΩIh(nhnh) = X a∈Nh nh(a)nh(a)ZΩ ϕafo ∀nh, nh∈Nh, which induces the norm knhkh=p(nh, nh)hdefined on Nh. We recall the following local error estimate. See [4, Thm. 4.4.4] or [9, Thm. 1.103] for a proof. Proposition 2.2. Under hypotheses (H1)–(H3), it follows that, (16) kϕ−IhϕkL∞(K)≤Capph2 Kk∇2ϕkL∞(K)∀K∈ Th,∀ϕ∈W2,∞(K), where Capp >0is independent of h. We next state the equivalence between the norms k·khand k·kin Nhand a discrete commuter approximation property for Ih. Proposition 2.3. Under hypotheses (H1)–(H3), it follows that, for all nh, nh∈Nh, (17) knhk≤knhkh≤51/2knhk and (18) knhnh−Ih(nhnh)kL1(Ω) ≤Capphknhkk∇nhk, where Capp >0is independent of h. APPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 5 Proof. We have knhk2=X a∈Nh n2 h(a)ZΩ ϕ2 a+X a6=e a∈Nh nh(a)nh(e a)ZΩ ϕaϕe a and knhk2 h=X a∈Nh n2 h(a)ZΩ ϕa. Since 1 = Pe a∈Nhϕe a, we write knhk2 h=X a,e a∈Nh n2 h(a)ZΩ ϕaϕe a=X a∈Nh n2 h(a)ZΩ ϕ2 a+X a6=e a∈Nh n2 h(a)ZΩ ϕaϕe a. Then knhk2 h−knhk2=X a>e a∈Nh (n2 h(a) + n2 h(e a)−2nh(a)nh(e a)) ZΩ ϕaϕe a =X a>e a∈Nh (nh(a)−nh(e a))2ZΩ ϕaϕe a≥0. From the above equality and Young’s inequality, we have knhk2 h=knhk2+X a>e a∈Nh (nh(a)−nh(e a))2ZΩ ϕaϕe a ≤ knhk2+ 2 X a>e a∈Nh (n2 h(a) + n2 h(e a)) ZΩ ϕaϕe a =knhk2+ 2 X a∈Nh n2 h(a)ZΩ ϕaX e a<a ϕe a+ 2 X e a∈Nh n2 h(e a)ZΩ ϕe aX a>ea ϕa ≤ knhk2+ 4knhk2≤5knhk2. We now prove (18). By using (16), we obtain kIh(nhnh)−nhnhkL1(Ω) =X K∈ThkIh(nhnh)−nhnhkL∞(K)ZK 1 ≤Capp X K∈Th h2 Kk∇2(nhnh)kL∞(K)ZK 1. Since nh, nh∈P1(K) on K∈ Th, we write ∇2(nhnh)=2 d X i,j=1 ∂inh∂jnh. Then, from (15) and on noting that ∇nh,∇nhare piecewise constant on each K∈ Th, we deduce that kIh(nhnh)−nhnhkL1(Ω) ≤Capp X K∈Th h2 Kk∇nhkL∞(K)k∇nhkL∞(K)ZK 1 ≤Capp X K∈Th h2 KZK|∇nh||∇nh| ≤CappCinv X K∈Th hKknhkL2(K)k∇nhkL2(K) ≤CappCinv hknhkk∇nhk, from which we conclude that (18) holds.  6 F. GUILL´ EN-GONZ ´ ALEZ AND J. V. GUTI´ ERREZ-SANTACREU We will need to use an (average) interpolation operator into Nhwith the following properties. In particular we use an extension of the Scott-Zhang interpolation operator to L1(Ω) function. We refer to [15,10] and [2]. Proposition 2.4. Under hypotheses (H1)–(H3), there exists an (average) interpolation operator Qhfrom L1(Ω) to Nhsuch that (19) kQhψkWs,p(Ω) ≤CstakψkWs,p(Ω) for s= 0,1and 1≤p≤ ∞, (20) kQh(ψ)−ψkWs,p(Ω) ≤Capph1+m−skψkWm+1,p(Ω) for 0≤s≤m≤1, and, for all ψ∈C∞(Ω) and nh∈Nh, (21) kQh(nhψ)−nhψkWs,p(Ω) ≤Capph1+m−sknhkWm,p(Ω)kψkWm+1,∞for 0≤s≤m≤1. The key point in proving a discrete maximum principle is the following property which is accomplished for right-angled simplexes assumed in (H2). Proposition 2.5. Under hypotheses (H1)–(H3), it follows that, for any diagonal nonnegative matrix D= diag(di)d i=1 (with di≥0), (22) D∇ϕa·∇ϕe a≤0a.e. in Ω if a6=e awith a,e a∈ Nh. Proof. For every right-angled d-simplex K∈ Thof vertices {ai}i=0,...,d with a0being the vertex supporting the right angle, we denote by Faithe opposite face to aiand by naithe exterior (to the d-simplex K) unit normal vector to the face Fai. Let b Kbe the reference unit d-simplex with vertices b a0=0and b ai=ei,i= 1,··· , d, where {ei}i=1,··· ,d is the canonical basis of Rd. Let FKbe the invertible affine mapping that maps b Konto Kdefined by FKb x=a0+BKb x, where BK∈Rd×d is orthogonal. Let bϕb ai(b x) = ϕai(FKb x). Then we have b ∇bϕb ai=−1 d|b Fb ai| |b K|nb ai. In particular, nb ai=−eiif i6= 0 and nb a0= [1,··· ,1]T. Thus, we obtain b ∇bϕb ai·b ∇bϕb aj=1 d2|b Fb ai||b Fb aj| |b K|2nb ai·nb aj≤0 if i6=j. Therefore, by means of the change of variable x=a0+BKb x, it follows that ∇ϕai=BKb ∇bϕb aiand hence D∇ϕai·∇ϕaj=DBKb ∇bϕb ai·BKb ∇bϕb aj=1 d2|b Fb ai||b Fb aj| |b K|2nT b aiBT KDBKnb aj≤0 if i6=j because, since BKis a orthogonal matrix, the inner products defined by Dand BT KDBKpreserves angles.  Remark 2.1. When D=Idwith Idbeing the d×didentity matrix, property (22)can be proved for nonobtuse triangulations [7]. Then property (22)can be somewhat seen a generalization restricted for right-angled triangulations. APPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 7 Let us now introduce the discrete Laplacian associated to the mass-lumping scalar product (·,·)h. For any Σh∈Nh, let −e ∆hΣh∈Nhsolve (23) −(e ∆hΣh, nh)h= (∇Σh,∇nh)∀nh∈Nh. We end up with a compactness result [1, Lm. 2.4] needed in proving the equivalence between problems (10) and (14). Theorem 2.1. Assume that (H1)-(H3) holds. Let 2d d+2 < ` < ∞. Suppose that {ρh,k}h,k≥0⊂ L2(0, T;L2(Ω)) is such that ρh,k(t, ·)∈Nhfor all t∈[0, T]and satisfies kρh,kkH1(0,T ;L`(Ω)) +kρh,kkL∞(0,T ;L2(Ω))∩L2(0,T;H1(Ω)) +ke ∆hρh,kkL2(0,T ;L2(Ω)) ≤Cdat. Then there exist a subsequence {ρh,k}h,k>0(not relabeled) and a limit function ρ, such that ρh,k →ρin L2(0, T, H1(Ω))-strongly as (h, k)→(0,+∞). Hereafter Cwill denote a generic constant whose value may change at each occurrence. This constant may depend on the data problem and the constants Cinv,Capp,Ccom and Cdat. 2.2. The numerical scheme. In order to avoid dense technical calculations, we assume for simplicity that each element K∈ Thhas its edges lined up with the axes. The numerical scheme relies on a finite-element method combined with a closed-nodal integration applied to the time-derivative and pressure-migration terms. Thus our numerical method which consists in finding nh,k ∈C1([0, T]; Nh) such that (24) (∂tnh,k, nh)h+ (∇Ih((nh,k)k),∇nh) + ν(∇nh,k,∇nh) = (G(p(nh,k))nh,k, nh)h∀n∈Nh nh,k(0) = n0 h,k, with p(nh,k) = k k−1(nh,k)k−1. Equivalently, we may write (24)1as (25) (∂tnh,k, nh)h+ (D(nh,k)∇nh,k,∇nh) + ν(∇nh,k,∇nh)=(G(p(nh,k))nh,k, nh)h, where D(nh,k) is a piecewise constant, d×ddiagonal matrix function with respect to Thdefined as follows. Let K∈ Thwith vertices {ai}i=0,··· ,d where a0corresponds to the right angle. Then (26) [D(nh,k)|K]ii =   (nh,k)k(ai)−(nh,k)k(a0) nh,k(ai)−nh,k(a0)if nh,k(ai)−nh,k(a0)6= 0, 0 if nh,k(ai)−nh,k(a0)=0. By the mean value theorem, one can write (27) [D(nh,k)|K]ii =k(nh,k)k−1(ξi), where ξi=αai+ (1 −α)a0for a certain α∈(0,1). The above choice for the sequence of {n0 h,k}h,k>0is as follows. Let {n0 k}k∈N⊂H1(Ω) ∩L∞(Ω) satisfy (7) and (9). Then we select n0 h,k =Qh(n0 k) so that (28) 0 ≤n0 h,k(a)≤Nmax(k)∀a∈ Nh,k∇n0 h,kk ≤ Cstabk∇n0 kk, (29) n0 h,k →n0 kin H1(Ω)-strongly as h→0. There is an additional technicality regarding the sequence of initial data that we must consider: (H4) Assume {nh,k}h,k>0to be such that (30) −(∇Ih(n0 h,k)k,∇nh)−ν(∇n0 h,k,∇nh)+(G(p(n0 h,k))n0 h,k, nh)h≥0∀nh∈Nhwith nh≥0. 8 F. GUILL´ EN-GONZ ´ ALEZ AND J. V. GUTI´ ERREZ-SANTACREU Remark 2.2. This last condition is related to imposing ∂tnh,k(0) ≥0which is crucial to prove the k→+∞limit. The existence and uniqueness of a solution to scheme (24) may be readily justified by Picard’s theorem. To be more precise, one may prove that there exists a time interval [0, Th) for which problem (24) is uniquely solvable. As a consequence of a priori energy estimates, which we shall prove in the next section, one deduces that Th=Tfor all h > 0. 2.3. Main result. We now are ready to state our main result of this paper. We shall prove that scheme (24) produces a sequence of discrete solutions which satifies a priori energy bounds uniform with respect to (h, k) allowing us to pass to the limit as (h, k)→(0,+∞) towards weak solutions of the Hele–Shaw-like system (10)-(13). Theorem 2.2. Assume that (H1)-(H3) hold. Then the discrete solution {(nh,k, ph,k)}h,k of (24) satisfies the following estimates, for all a∈ Nhand t∈[0, T]: 0≤nh,k(a, t)≤Nmax(k), 0≤p(nh,k(a, t)) ≤Pmax, ∂tnh,k(a, t)≥0, ∂tp(nh,k(a, t)) ≥0. Furthermore, {nh,k,Ih((nh,k)k)}h,k converges towards weak solutions (n∞, p∞)of problem (10)-(13) in the sense that nh,k →n∞in L∞(0, T;H1(Ω))-weakly-?and in Lp((0, T)×Ω)-strongly, and Ih((nh,k)k)→p∞in L∞(0, T;H1(Ω))-weakly-?and in Lp((0, T)×Ω)-strongly, for any 1<p<∞provided that (H5) k h →0as (h, k)→(0,+∞). 3. Proof of Theorem 2.2 3.1. A priori energy estimates. Our goal is to prove a priori energy estimates for the discrete solution nh,k of (24) independent of (h, k). This first lemma will be focused on proving a discrete maximum principle for nh,k based on the hypothesis of right-angled triangulations. Moreover, some a priori energy estimates are obtained. Lemma 3.1. Assume that (H1)-(H3) hold. Then the solution nh,k of scheme (24)satisfies (31) 0 ≤nh,k(a, t)≤Nmax(k)∀a∈ Nhand ∀t≥0, and (32) knh,kkL∞(0,T ;L2(Ω)) +knh,kkL2(0,T ;H1(Ω)) ≤C, where C > 0is independent of (h, k). Proof. We first proceed to verify (31). In doing so, we introduce a modification to scheme (25) which truncates the nonlinear diffusion term as follows: (33) (∂tnh,k, nh)h+ (D([nh,k]T)∇nh,k,∇nh) + ν(∇nh,k,∇nh) = (G(p([nh,k]T))nh,k, nh)h, where [nh,k]Tis the usual truncation of nh,k from below by 0 and from above by Nmax(k). Again, by means of Picard’s theorem, one has the existence and uniqueness of a solution nh,k to (33). Let nmin h,k =Ih(n− h,k)∈Nhbe defined as nmin h,k =X a∈Nh n− h,k(a)ϕa, APPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 9 where n− h,k(a) = min{0, nh,k(a)}. Analogously, one defines nmax h,k =Ih(n+ h,k)∈Nhas nmax h,k =X a∈Nh n+ h,k(a)ϕa, where n+ h,k(a) = max{0, nh,k(a)}. Notice that nh,k =nmin h,k +nmax h,k . On choosing nh=nmin h,k in (33), it follows that (34) 1 2 d dtknmin h,k k2 h+ (D([nh,k]T)∇nh,k,∇nmin h,k ) + ν(∇nh,k,∇nmin h,k ) = kG(p([nh,k]T))1/2nmin h,k k2 h ≤G(0)knmin h,k k2 h. Next observe that (D([nh,k]T)∇nh,k,∇nmin h,k ) = (D([nh,k]T)∇nmin h,k ,∇nmin h,k )+(D([nh,k]T)∇nmax h,k ,∇nmin h,k ) =kD([nh,k]T)1/2∇nmin h,k k2+X a6=˜ a∈Nh n− h,k(a)n+ h,k(e a)(D([nh,k]T)∇ϕa,∇ϕe a). Then, using the fact that n− h,k(a)n+ h,k(e a)≤0 if a6=˜ aand that D([nh,k]T) is a nonnegative diagonal matrix function, one deduces, from (22), that D([nh,k]T)∇ϕa·∇ϕe a≤0∀a6=e a∈ Nh and thereby (35) (D([nh,k]T)∇nh,k,∇nmin h,k )≥ kD([nh,k]T)1/2∇nmin h,k k2. Analogously, one obtains (36) ν(∇nh,k,∇nmin h,k )≥νk∇nmin h,k k2, where we have used again (22) but now for D=Id, with Idbeing the d×dunit matrix. Inserting (35) and (36) into (34) yields 1 2 d dtknmin h,k k2 h+kD([nh,k]T)1/2∇nmin h,k k2+νk∇nmin h,k k2≤G(0)knmin h,k k2 h. By Gr¨onwall’s lemma, we have nmin h,k (t)≡0 in Ω, for any t≥0, since nmin h,k (0) ≡0 in Ω; thereby this implies 0 ≤nh,k in (31). For the other inequality nh,k ≤Nmax(k) in (31), we proceed in a similar fashion. In this case, one chooses nh= (nh,k −Nmax(k))max in (33) and takes into account that G(p([nh,k]T))nh,k(nh,k −Nmax(k))max ≡0 due to p([nh,k]T) = Pmax if nh,k ≥Nmax(k). It should be noted that any solution nh,k of the modified scheme (33) satisfies the discrete maximum principle (31), and consequently [nh,k]T≡nh,k; hence nh,k satisfies the non-truncated scheme (24) as well. Finally, by uniqueness of solutions for scheme (24), the solution of (24) takes values between 0 and Nmax(k); that is (31). Now selecting nh=nh,k in (25) and invoking Gr¨onwall’s lemma, the following energy estimate holds, for all t∈[0, T]: (37) 1 2knh,k(t)k2 h+ZT 0kD(nh,k)1/2∇nh,kk2+νZT 0k∇nh,kk2≤exp(2G(0)T)1 2kn0 h,kk2 h. Then the weak estimates (32) are deduced from (37) and (17).  A discrete maximum principle for (nh,k)k−1and (nh,k)kfollows as a direct consequence of (31). 16 F. GUILL´ EN-GONZ ´ ALEZ AND J. V. GUTI´ ERREZ-SANTACREU On testing (57) against Qh(Ih(nk h,k)ψ) with ψ∈C∞ c(Ω ×[0, T]) such that ψ≥0, it follows that (59) ZT 0 (∂tnh,k ∗ρε,Qh(Ih(nk h,k)ψ))h=ZT 0 ((G(p(nh,k))nh,k)∗ρε,Qh(Ih(nk h,k)ψ))h −ZT 0 (∇(Σ(nh,k)∗ρε),∇Qh(Ih(nk h,k)ψ)). Since (∂tnh,k ∗ρε,Qh(Ih(nk h,k)ψ))h≥0, we obtain (60) 0 ≤ZT 0 ((G(p(nh,k))nh,k)∗ρε,Qh(Ih(nk h,k)ψ))h−ZT 0 (∇(Σ(nh,k)∗ρε),∇Qh(Ih(nk h,k)ψ)). Taking the limit as (h, k)→(0,∞) yields (61) ZT 0 ((G(p(nh,k))nh,k)∗ρε,Qh(Ih(nk h,k)ψ))hdt→ZT 0 ((G(p∞)n∞)∗ρε, p∞ψ)dt and (62) ZT 0 (∇(Σ(nh,k)∗ρε),∇Qh(Ih(nk h,k)ψ))dt→ZT 0 (∇((p∞+νn∞)∗ρε),∇(p∞ψ))dt. In order to prove (61), we use the decomposition (uh, vh)h= (uh, vh) + (Ih(uhvh)−uhvh,1) for uh= (G(p(nh,k))nh,k)∗ρεand vh=Qh(Ih(nk h,k)ψ) to write ZT 0 ((G(p(nh,k))nh,k)∗ρε,Qh(Ih(nk h,k)ψ))h=ZT 0 ((G(p(nh,k))nh,k)∗ρε,Qh(Ih(nk h,k)ψ)) +ZT 0 (Ih((G(p(nh,k)))nh,k)∗ρεQh(Ih(nk h,k)ψ)−(G(p(nh,k))nh,k)∗ρεQh(Ih(nk h,k)ψ),1). Then, it follows from (21), (49) and (52) that the first term converges to RT 0((G(p∞)n∞)∗ ρε, p∞ψ)dt, and, on noting that k∇Qh(Ih(nk h,k)ψ)k ≤ Ck∇Ih(nk h,k)kkψkL∞+CkIh(nk h,k)kL∞k∇ψk from (19), and on recalling (18) and (46), the second term converges to zero; thereby (61) holds. In order to prove (3.2.4), we write ZT 0 (∇(Σ(nh,k)∗ρε),∇Qh(Ih(nk h,k)ψ)) = ZT 0 (∇(Σ(nh,k)∗ρε),∇(Ih(nk h,k)ψ)) −ZT 0 (∇(Σ(nh,k)∗ρε),∇(Ih(nk h,k)ψ−Qh(Ih(nk h,k)ψ))). Then, it follows from (58), (48) and (51) that the first term converges to RT 0(∇((p∞+νn∞)∗ ρε),∇(p∞ψ))dt, and on noting that k∇(Ih(nk h,k)ψ−Qh(Ih(nk h,k)ψ))k ≤ hk∇Ih(nk h,k)kkψkW2,∞(Ω) from (21) and on recalling (46), the second term converges to zero; thereby (3.2.4) holds. Thus, by applying the previous convergences (61) and to (60), we arrive at 0≤ZT 0 (G(p∞)n∞∗ρε, p∞ψ)−(∇((p∞+νn∞)∗ρε),∇(p∞ψ))dt, and finally (55) holds by taking the limit as ε→0. APPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 17 •We proceed to prove (56). Write the first term on the right-hand side of (59) as (63) ZT 0 (∂tnh,k ∗ρε,Qh(Ih(nk h,k)ψ))h=ZT 0 (∂tnh,k ∗ρε,Ih(nk h,k)ψ)h +ZT 0 (∂tnh,k ∗ρε,Qh(Ih(nk h,k)ψ)−Ih(nk h,k)ψ)h. These two terms are handled as follows. For the second term of (63), we have, by (17), (21) and (41), that ZT 0 (∂tnh,k ∗ρε,Qh(Ih(nk h,k)ψ)−Ih(nk h,k)ψ)hds→0 as (h, k)→(0,∞). For the first term of (63), we have that, for each a∈ Nh, (∂tnh,k(a, t)∗ρε)nk h,k(a, t) = nk h,k(a, t)ZR ∂tnh,k(a, s)ρε(t−s)ds =ZR nk h,k(a, s)∂tnh,k(a, s)ρε(t−s)ds +ZR (nk h,k(a, t)−nk h,k(a, s))∂tnh,k(a, s)ρε(t−s)ds. On integrating by parts in time and using (31) and (39), we obtain ZR nk h,k(a, s)∂tnh,k(a, s)ρε(t−s)ds=1 k+ 1 ZR nk+1 h,k (a, s)∂tρε(t−s)ds→0 as (h, k)→(0,∞). Furthermore, for s > t, we have that nk h,k(a, t)−nk h,k(a, s)≤0 owing to (40). Then, if we choose supp(ρε)⊂(−ε, 0), then ZR (nk h,k(a, t)−nk h,k(a, s))∂tnh,k(a, s)ρε(t−s)ds≤0. Letting first (h, k)→(0,∞) in (63) and then ε→0, we obtain (56) by repeating the arguments that led to (55). As a result of (55) and (56), we note that (64) ZT 0 (G(p∞)n∞, p∞ψ)−(∇(p∞+νn∞),∇(p∞ψ))ds= 0 is satisfied for all ψ∈C∞ c(Ω ×[0, T]) with ψ≥0, and therefore it also holds for all ψ∈C∞ c(Ω × [0, T]). From the fact that p∞∇n∞= 0 and p∞≥0 a.e. in Ω×(0, T), we also deduce that ∇p∞·∇n∞= 0 a.e. in Ω ×(0, T). As a consequence, the above variational equation (64) is equivalent to ZT 0 (G(p∞)n∞, p∞ψ)−(∇p∞,∇(p∞ψ))ds= 0 which, taking into account (53), implies (13) in the distributional sense. 4. An algorithm on unstructured meshes In order to avoid using structured meshes, we propose the following scheme. Find nh,k ∈ C1([0, T]; Nh) such that (65) (∂tnh,k, nh)h+k((nh,k)k−1∇nh,k,∇nh) + ν(∇nh,k,∇nh) = (G(p(nh,k))nh,k, nh)h∀nh∈Nh, nh,k(0) = n0 h,k. 18 F. GUILL´ EN-GONZ ´ ALEZ AND J. V. GUTI´ ERREZ-SANTACREU Equivalently, we may write (65)1as (∂tnh,k, nh)h+ (nh,k∇p(nh,k),∇nh) + ν(∇nh,k,∇nh) = (G(p(nh,k))nh,k, nh)h. Here the finite-element space Nhis constructed over a family of triangulations {Th}h>0of Ω being shape-regular, quasi-uniform and with acute angles. This acuteness property implies (22) for the particular case where Dis the d×didentity matrix [7]. We summarize the properties of scheme (65) in the following theorem. Theorem 4.1. Suppose that (H1)-(H4) are satisfied. Then scheme (65)satisfies the following properties. For all a∈ Nhand t≥0, we have: 0≤nh,k(a, t)≤Nmax(k) 0≤nk h,k(a, t)≤PmaxNmax(k), ∂tnh,k(a, t)≥0, ∂tnk h,k(a, t)≥0, and the a priori estimates: knh,kkL∞(0,T ;L2(Ω))∩L2(0,T;H1(Ω)) ≤C, k∂tnh,kkL∞(0,T ;L1(Ω)) +k∂tnk h,kkL1(0,T ;L1(Ω)) ≤C, with C > 0being a constant independent of (h, k). Proof. Full details of the proof are left to the interested reader since it follows mutatis mutandis the same arguments as for scheme (24).  Corollary 4.1. Under hypotheses (H1)-(H4), it follows that (66) X K∈ThZK>|∂xink h,k(x)|2+ZK<|∂xiIhnk h,k(x)|2dx≤C, where K>=(x∈K:nk−1 k,h (ξi) nk−1 h,k (x)>1) and K<=(x∈K:nk−1 k,h (ξi) nk−1 h,k ,(x)<1), with C > 0being a constant independent of (h, k). Proof. Choose ¯nh=Ih(nk h,k) to get (67) (∂tnh,k,Ih(nk h,k))h+((nh,k)k−1∇nh,k,∇Ih(nk h,k))+ν(∇nh,k,∇Ih(nk h,k)) = (G(p(nh,k))nh,k,Ih(nk h,k))h. It follows immediately from (39) and (41) that (68) (∂tnh,k,Ih(nk h,k))h≥0, and from (26) that (69) ν(∇nh,k,∇Ih(nk h,k)) = ν(D(nh,k)∇nh,k,∇nh,k)≥0. Combining (67)-(69) yields on noting (31) and (39) that ((nh,k)k−1∇nh,k,∇Ih(nk h,k)) ≤G(0)|Ω|Nmax(k)2Pmax. APPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 19 Finally, we invoke again (26) and recall (27) to set k((nh,k)k−1∇nh,k,∇Ih(nk h,k)) = (∇nk h,k,∇Ih(nk h,k)) =X K∈ThZK nk−1 k,h (ξi) nk−1 h,k (x)|∂xink h,k(x)|2+nk−1 h,k (x) nk−1 k,h (ξi)|∂xiIhnk h,k(x)|2!dx. This completes the proof via the definitions of K<and K>. Remark 4.1. Unfortunately, convergence for scheme (65)is not clear because estimate (66)does not provide enough control over the gradient of {nk h,k}h,k or {Ihnk h,k}h,k in order to obtain compactness and therefore to pass to the limit as (k, h)→(0,+∞). 5. Numerical simulation 5.1. Temporal integration. It is assumed here for simplicity that we have a uniform partition of [0, T] into Mpieces, with time step size τ=T/M and the time values (tm=mτ)M m=0. To simplify the notation let us denote δtnm+1 =nm+1 −nm τ. First we present a first-order time integration for scheme (65). Algorithm 1: Linear semi-implicit time-stepping scheme Step (m+ 1): Given nm h,k ∈Nh, find nm+1 h,k ∈Nhsolving the algebraic linear system (70) (δtnm+1 h,k , nh)h+k((nm h,k)k−1∇nm+1 h,k ,∇nh) +ν(∇nm+1 h,k ,∇nh) = (G(p(nm h,k))nm h,k, n)h, for all nh∈Nh. 5.2. Computational experiments. In this section, we present several numerical experiments to test the algorithm presented herein. To do this, we consider the evolution of problem (1)-(3) with n0(x, y) = α e−(x2+y2) on the computational domain Ω = (−10,10) ×(−10,10) with α > 0. In the numerical setting, we construct a structured triangulation partitioning the edges of square into 100 subintervals, corresponding with the mesh size h= 0.12582843 and the time step size is τ= 10−5. The choice of the time step τis such that it helps to mitigate the possibly numerical deviation of the d-simplexes K∈ Thfrom the right-angled structure. The resulting matrix is strictly diagonally dominant. Our intention is to illustrate the behavior exhibited by the solution to problem (1)-(3) when the diffusion coefficient ν, the parameter kand the homeostatic pressure Pmax vary. We will set αand Pmax to be 1 and kto be 100 if not stated otherwise. Moreover, we consider G(p) = 200 πarctan(4(Pmax −p)+), 5.2.1. Analysis of the effect of α(contraction/dilation coefficient of the initial datum). In this test we choose α= 0.5 and 1. We are interested in comparing the evolution of the density nand the pressure p(n) when the maximum of the initial density takes different values. In particular, we have for α= 0.5 that the maximum value of n0occurs only at the point (0,0) and is 0.5, hence Nmax(k) remains below of 1. We thus observe that the maximum increases without modifying essentially the exponential shape of the initial datum n0until reaches Nmax(k) = 1. Once the density takes the value 1 at t= 0.01583 the measure of points at which the density reaches the maximum grows radially around (0,0) due to the fact that the pressure starts increasing and pushes forward the 20 F. GUILL´ EN-GONZ ´ ALEZ AND J. V. GUTI´ ERREZ-SANTACREU tumor cells. Then the exponential structure of the initial datum n0becomes a traveling wave shape which moves outwards as tincreases. This behavior causes that the evolution of the interface is delayed concerning the case α= 1 as shown in Figures 1and 2since the maximum value 1 is reached from the beginning. Figure 3represents the difference between the density and the pressure at times t= 0.1, 0.2, 0.3 and 0.4, and indicates that the pressure is responsible for the advance of the tumor cells which is deduced from the annulus shape of the difference. Figure 1. Evolution of the density at times t= 0.1,0.2,0.3,0.4 for α= 0.5 (top) and 1 (bottom). Figure 2. Evolution of the pressure at times t= 0.1,0.2,0.3,0.4 for α= 0.5 (top) and α= 1 (bottom). 5.2.2. Analysis of the effect of ν(active motion coefficient). Now we set Pmax = 1 and take different values of ν= 0, 0.5 and 1. The evolution of the density nh,k is shown in Figure 4where we see that the velocity of propagation of the tumor cells increases with respect to νas noted for times t= 0.1, 0.2, 0.3 and 0.4. Moreover, no particular differences have been observed in the width of the interface between the tumor and pre-tumor cells for the different values of ν. 5.2.3. Analysis of the effect of k.In this simulation we select k= 10 and 1000. The first thing we have noted is that there is a dependence between kand τwhich has been taken 0.5·10−5. As can be seen in Figure 5, there are no particular differences for k= 10 and 1000 at times t= 0.1, 0.2, 0.3 and 0.4. APPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 21 Figure 3. Evolution of the difference between the density and pressure at times t= 0.1,0.2,0.3,0.4 for α= 0.5 (top) and α= 1 (bottom). Figure 4. Comparison of the density at times t= 0.1,0.2,0.3,0.4 for different ν= 0 (top), 0.5 (middle), 1 (bottom). Figure 5. Comparison of the density at times t= 0.1,0.2,0.3,0.4 for different k= 10 (top) and 1000 (bottom). 22 F. GUILL´ EN-GONZ ´ ALEZ AND J. V. GUTI´ ERREZ-SANTACREU 5.2.4. Analysis of the effect of Pmax.Let us take Pmax = 10 and 30. Figure 6shows that the dynamics is sensitive to the different values for the homeostatic pressure. We highlight that, for Pmax = 30, the evolution of the interphase is faster than the one for Pmax = 10. Moreover, the shape of the interphase seems different as depicted in Figure 6for times t= 0.1, 0.2, 0.3 and 0.4. Figure 6. Comparison of the density at times t= 0.1,0.2,0.3,0.4 for different Pmax = 10 (top) and 30 (bottom). References [1] Becker, R.; Feng, X.; Prohl, A. Finite element approximations of the Ericksen-Leslie model for nematic liquid crystal flow. SIAM J. Numer. Anal. 46 (2008), no. 4, 1704–1731. [2] Bertoluzza, S The discrete commutator property of approximation spaces. C. R. Acad. Sci. Paris S´er. I Math. 329 (1999), no. 12, 1097–1102. [3] Betteridge, R.; Owen, M. R.; Byrne, H. M.; Alarc´ on, T.; Maini, P. K.The impact of cell crowding and active cell movement on vascular tumour growth. Netw. Heterog. Media 1 (2006), no. 4, 515?535. [4] Brenner, S. C.; Scott, L. R.,The mathematical theory of finite element methods, Third edition. Texts in Applied Mathematics, 15. Springer, New York, 2008. [5] Br´ u, A.; Albertos, S.; Subiza, J. L.; Asenjo, J. A.; Broe, I. The universal dynamics of tumor growth. Biophys. J. 85 (2003), no. 5, 2948–2961. [6] Byrne, H. M.; Drasdo, D.,Individual-based and continuum models of growing cell populations: a comparison., J. Math. Biol. 58 (2009), no. 4-5, 657–687. [7] Ciarlet, P.G.; Raviart, P.-A. Maximum principle and uniform convergence for the finite element method, Comput. Methods Appl. Mech. Engrg. 2 (1973) 17–31. [8] Drasdo, D.; Hoehme, S. Modeling the impact of granular embedding media, and pulling versus pushing cells on growing cell clones. New J. Phys. 14 (2012) 055025. [9] Ern, A; Guermond, J.-L.,Theory and practice of finite elements, Applied Mathematical Sciences, 159. Springer-Verlag, New York, 2004. [10] Girault, V.; Lions, J.-L.,Two-grid finite-element schemes for the transient Navier-Stokes problem. M2AN Math. Model. Numer. Anal. 35 (2001), no. 5, 945–980. [11] Perthame, B.; Quir´ os, F.; Tang, M.; Vauchelet, N.,Derivation of a Hele-Shaw type system from a cell model with active motion, Interfaces Free Bound. 16 (2014), no. 4, 489–508. [12] Perthame, B.; Quir´ os, F.; V´ azquez, J. L.,The Hele-Shaw asymptotics for mechanical models of tumor growth, Arch. Ration. Mech. Anal. 212 (2014), no. 1, 93–127. [13] J. Ranft, M. Basan, J. Elgeti, J.-F. Joanny, J. Prost and F. J¨ ulicher,Fluidization of tissues by cell division and apoptosis, Proceedings of the National Academy of Sciences, 107 (2010), no. 49, 20863–20868. [14] Saut, O.; Lagaert, J.-B.; Colin, T.; Fathallah-Shaykh, H. M. A multilayer grow-or-go model for GBM: effects of invasive cells and anti-angiogenesis on growth. Bull. Math. Biol. 76 (2014), no. 9, 2306–2333. [15] Scott, L.R.; Zhang, S. Finite element interpolation of non-smooth functions satisfying boundary conditions. Math. Comp. 54 (1990) 483–493. [16] Simon, J.Compact sets in the space Lp(0, T ;B). Ann. Mat. Pura Appl. (4) 146 (1987), 65–96. APPROXIMATION OF A HELE-SHAW-LIKE SYSTEM 23 (†)Dpto. E.D.A.N. and IMUS, Universidad de Sevilla, Aptdo. 1160, 41080 Sevilla, Spain. E-mail: [email protected] (‡)Dpto. de Matem´ atica Aplicada I, E. T. S. I. Inform´ atica, Universidad de Sevilla. Avda. Reina Mercedes, s/n. E-41012 Sevilla, Spain. E-mail: [email protected]