Impact Factor 1.895 | CiteScore 2.24
More on impact ›

Original Research ARTICLE

Front. Phys., 09 July 2019 |

Numerical Solutions for Time-Fractional Cancer Invasion System With Nonlocal Diffusion

J. Manimaran1, L. Shangerganesh1, Amar Debbouche2* and Valery Antonov3
  • 1Department of Applied Sciences, National Institute of Technology Goa, Goa, India
  • 2Department of Mathematics, Guelma University, Guelma, Algeria
  • 3Department of Mathematics, Peter the Great Saint Petersburg Polytechnic University, Saint Petersburg, Russia

This article studies the existence and uniqueness of a weak solution of the time-fractional cancer invasion system with nonlocal diffusion operator. Existence and uniqueness results are ensured by adapting the Faedo-Galerkin method and some a priori estimates. Further, finite element numerical scheme is implemented for the considered system. Finally, various numerical computations are performed along with the convergence analysis of the scheme.

1. Introduction

In the past few decades, a large number of mathematical models have been applied for biological studies. In addition, mathematical models give a deeper conceptual understanding of behavioral dynamics of complex systems. Some of the advantages of mathematical models include cost efficient experiments, which can be performed speedily without disturbing biological variants. Cancer is a disease defined by a normal cell which starts replicating out-of-control. Over the years, cancer modeling has gained popularity with applied mathematicians because of its challenges, resulting in numerous research findings on the dynamics of tumor invasion. Some of these propositions are available in the literature to acquaint oneself with the developments in cancer modeling (see for instance [14]).

On the other hand, fractional differential equations (FDEs) have been extensively used for constructing biological models and other areas of science and engineering. We refer the following monographs [57] and research articles [812] that have explored recent developments using FDEs. Biological phenomena have an anomalous diffusion property which includes heterogeneous systems that are witnessed in porous materials (see for example [13, 14]). Linear and nonlinear models of anomalous diffusion, which have been experimented by researchers could not do justice to the biological phenomena. But the fractional models have contributed to replicate the biological phenomena with a greater accuracy. Diffusions in biological tissue has characterized as anomalous. Therefore, it has been shown to be best described using fractional calculus tools. It means equations involving non-integer derivatives and integrals. Cancer models adapting fractional differential equations are studied by Ahmed et al. [15] and Iyiola and Zaman [16] and also see the references there in.

Theoretical and numerical analysis of fractional partial differential equations (FPDEs), concerned with only very few articles, are available in the literature. Alikhanov [17] applied the method of energy inequalities to obtain the existence of solutions for a time-fractional boundary value problem of the diffusion-wave equation. Jiao and Zhou [18] established an existence result for a fractional boundary value problem with the application of the critical point theory. Further, Zhou et al. [19] analyzed the time-fractional reaction-diffusion equation under nonlocal boundary condition. Zhou and Peng [11] established the existence, uniqueness and regularity of time-fractional Navier-Stokes equations. Further, they ensured the existence of weak solutions and also provided the sufficient conditions for optimal controls in Zhou and Peng [12]. Finite volume method [2022], meshless method [23, 24], finite difference method [2528], finite element method [2931] and the spectral method [3235] are widely preferred numerical methods in the literature to solve fractional partial differential equations. Finite element methods have become popular for numerical simulations of time-fractional diffusion equations due to their good approximation and feasibility to work with any domains. Recently, Esen et al. [36] studied the numerical solutions of time-fractional diffusion equations and diffusion-wave equations using Galerkin finite element method. Jin et al. [29] analyzed the numerical solutions of multiple time-fractional derivative using the Galerkin finite element method. Wang et al. [37] combined second-order time approximation with the finite element method to solve nonlinear fractional Cable equation. Liu et al. [38] used fully discrete mixed finite element scheme to study the second order convergence for nonlinear time-fractional diffusion problem with fourth-order derivative term. Jin et al. [39] solved proposed Crank-Nicolson-Galerkin finite element scheme to solve the linear time FPDEs. Kumar et al. [40] proposed Crank-Nicolson-Galerkin finite element scheme to solve the time-fractional nonlinear diffusion equation using Newton's algorithm. However, according to author's knowledge there is no paper available in the literature to study the fractional order cancer invasion system using finite element method.

Recently fractional reaction-diffusion systems are applied for many applications in science and engineering. Fractional models are proposed and used in chemical reactions, propagation phenomena, transport systems, pattern formation processes and spatiotemporal distribution of species [4145] and references therein. In this connection, we are interested to study and analyze the time-fractional cancer invasion model with nonlocal diffusion operator. Existence and uniqueness of a weak solution and various numerical simulations are presented for the below considered model. We considered a mathematical model proposed in Solis and Delgadillo [46] with four unknown variables namely two cancer cells density, normal cells density and acidification medium concentration. Further, we extend the same model for fractional differential equations and we show the importance of fractional derivatives using numerical simulations. The dynamics of cancer invasion system with time-fractional is governed by the following nonlocal diffusion system:

tαu1d1(l(u1))Δu1=            u1(1u1)β1u1u2ρu1γ1u1u3                 in QT,tαu2d2(l(u2))Δu2=            r2u2(1u2)β2u1u2+ρu1δ1u2u3           in  QT,tαu3=r3u3(1u3)            γ2u1u3δ2u2u3σu3u4                                   in QT, tαu4d4(l(u4))Δu4=             ξ(u1+u2u4)                                                                 in QT,}    (1.1)

with initial and boundary conditions

uj(x,0)=uj,0(x), j=1,2,3,4 in Ω, ui(x,t)=0, i=1,2,4 in ΣT,

where QT = Ω × (0, T), ΣT = ∂Ω × (0, T), T > 0 is final time and α ∈ (0, 1]. Here, Ω is a bounded domain in ℝN with smooth boundary ∂Ω. The unknown functions u1(x, t) and u2(x, t), respectively, describe the density of two types of cancer cells. Further, u3(x, t) and u4(x, t), respectively, represent the density of normal cells and medium acidification concentration due to excess H+ ions. The constants β1 and β2, respectively, denote the rates of interaction and the positive constant ρ delineates the intrinsic mutation rate of cancer cells. Furthermore, γ1 and δ1 represent the rate of consumption of cancer cell populations. The proliferation rate of cancer cells is given by r2 ≥ 0 and r3 ≥ 0. Here, ξ represents the production rate of the H+ions. Moreover, γ2 and δ2, respectively, denote the interaction rate of two types of cancer cells with normal cells and σ denotes the degradation rate of normal cells due to acidification. In (1.1), the diffusion rates di: ℝ → ℝ are the Lipschitz continuous functions with di(ξ) ≥ mi > 0 where i = 1, 2, 4. Further, di, i = 1, 2, 4 are taken to be depend on the whole of each population in the domain rather than on the local density. From a physical point of view of biological models, especially migration of cancer cells through normal cell is more like movement in a porous medium. Therefore, we consider the cell random motility to be a function of unknowns, see for example Szymańska et al. [47]. Therefore, it is more realistic to work with density dependent diffusion like nonlocal diffusion instead of linear diffusion function. Further, we assume that linear continuous nonlocal operator l(s) ∈ xs(L2(Ω))′ where s ∈ ℝ. This work investigates existence and uniqueness of a weak solution and numerical solutions for the time-fractional cancer invasion system (1.1).

It should be remarked that throughout the paper, we use the Caputo sense fractional derivatives for time. The main advantage the Caputo derivative, we can use initial conditions as in integer order derivatives. However, for more details we refer the interested readers to the book [48].

The rest of the manuscript is arranged as follows. In section 2, we present some preliminaries of fractional calculus and existence and uniqueness of weak solution of (1.1) using the Faedo-Galerkin approximation method and priori estimates. In section 3, we give the variational formulation of (1.1), finite element discretization and temporal discretization. Finally, in section 4, we present the convergence study of the numerical scheme and some computations with various numerical experiments.

2. Existence and Uniqueness

The goal of this section is to prove existence and uniqueness of a weak solution of nonlocal density dependent diffusion cancer invasion parabolic system with time-fractional derivative (1.1). By adopting the Faedo-Galerkin approximation method and deriving uniform a priori estimates for approximation solution, we show the existence of a weak solution in appropriate solution space. Here, we use the same notations and definitions as in Zhou et al. [19, 49]. Further, in order to avoid too many notations, we use a generic constant C instead of different constants.

Theorem 2.1 ([6]). Consider the fractional ordinary differential equation (FODE)

D0Cxαu(x)=f(x,u(x)), 0<x<T,u(0)          =b0,}    (2.1)

where 0 < α < 1 and b0 ∈ ℝ be a given constant. Suppose U be an open and connected set inand Ω = [0, T] × U. Assume that f(x, y):(0, T) × U → ℝ be a continuous function satisfying the Lipschitz condition. Then for any (x0, y0) ∈ Ω, there exists h > 0 such that the real interval [x0h, x0 + h] ⊂ (0, T) and there exists a unique solution u(x): [x0h, x0 + h] → U for (2.1) such that u(x)C[x0h, x0 + h], and for any x ∈ [x0h, x0 + h].

Lemma 2.1. Suppose u:[0, T] → X where X: = L2(Ω) is a real Hilbert space. Assume that there exists fractional derivative of u in the Caputo sense, then the following inequality holds true:


Lemma 2.2. Let α ∈ (0, 1) and a non-negative integrable function c1(t) for t ∈ [0, T] satisfies the inequality

D0Ctαu(t)c1(t),    (2.2)

for almost all t ∈ [0, T]. Then

u(t)u(0)+1Γ(α)0t(ts)α1c1(s)ds.    (2.3)

Lemma 2.3 ([50], p.9). Let α ∈ (0, 1). Suppose u, v are two integrable functions, v is nondecreasing and g is a continuous function in [a, b]. If




where Eα(·) is one parameter Mittag-Leffler function.

Lemma 2.4. Let α ∈ (0, 1). Suppose u(·) is a non-negative, absolute continuous function on [0, T], which satisfies for a.e. t the following differential inequality

D0Ctαu(t)Cu(t),    (2.4)

for constant C ≥ 0. Then


Proof. From (2.4),


where we use the Lemma 2.2 Using the Lemma 2.3, we get


Assume that X0, X1, X are Hilbert spaces. The Fourier transform of u : ℝ → X1 is defined by û(τ)=-e-2iθtτu(t)dt (See [51]). Then, we have


For 0 < α ≤ 1, define a Hilbert space


endowed with the norm


For any set J ⊂ ℝ, define a subspace WJα of Wα (see p. 274, [49]) as with support contained in J:


Further, we use the space


throughout the article.

Theorem 2.2 ([51]). Assume that X0XX1 is continuous and X0X is compact. Then for any bounded set J and α > 0, WJα(,X0,X1)L2(,X) is compact.

Without loss of generality, we rewrite the nonlocal density dependent diffusion cancer invasion parabolic system (1.1) with time-fractional derivative in the following form:

tαu1d1(l(u1))Δu1+G1(x,t,u1,u2,u3)=(1ρ)u1in QT,tαu2d2(l(u2))Δu2+G2(x,t,u1,u2,u3)=r2u2+ρu1in QT,tαu3+G3(x,t,u1,u2,u3,u4)=r3u3in QT,tαu4d4(l(u4))Δu4=ξ(u1+u2u4)in QT,}    (2.5)


G1(x,t,u1,u2,u3)        =u1(u1+β1u2+γ1u3),G2(x,t,u1,u2,u3)        =u2(β2u1+r2u2+δ1u3),G3(x,t,u1,u2,u3,u4)=u3(γ2u1+δ2u2+r3u3+σu4).

Theorem 2.3. Suppose the initial conditions uj,0, j = 1, 2, 3, 4 are in L2(Ω). Then the system (2.5) admits a weak solution u1, u2, u3, u4, which satisfies the following conditions:


such that for every ϕjL2(0,T;H01(Ω)),j=1,2,3,4,

0Ttα(u1,ϕ1)dt+d1(l(u1))QTu1ϕ1dxdt           +QTG1(x,t,u1,u2,u3)ϕ1dxdt           =(1ρ)QTu1ϕ1dxdt,0Ttα(u2,ϕ2)dt+d2(l(u2))QTu2ϕ2dxdt           +QTG2(x,t,u1,u2,u3)ϕ2dxdt           =QT(r2u2+ρu1)ϕ2dxdt,0Ttα(u3,ϕ3)dt+QTG3(x,tu1,u2,u3,u4)ϕ3dxdt           =QTr3u3ϕ3dxdt,0Ttα(u4,ϕ4)dt+d4(l(u4))QTu4ϕ4dxdt           =QTξ(u1+u2u4)ϕ4dxdt.    (2.6)

Now, we use the following regularized system in order to find weak solutions of the system (2.5). For ϵ > 0,

tαu1ϵd1(l(u1ϵ))Δu1ϵ+G1,ϵ(x,t,u1ϵ,u2ϵ,u3ϵ)=u1ϵρu1ϵin QT,tαu2ϵd2(l(u2ϵ))Δu2ϵ+G2,ϵ(x,t,u1ϵ,u2ϵ,u3ϵ)=r2u2ϵ+ρu1ϵin QT,tαu3ϵ+G3,ϵ(x,t,u1ϵ,u2ϵ,u3ϵ,u4ϵ)=r3u3ϵin QT,tαu4ϵd4(l(u4ϵ))Δu4ϵ=ξ(u1ϵ+u2ϵu4ϵ)in QT,}    (2.7)

with initial and boundary conditions

ujϵ(x,0)=uj,0ϵ(x),j=1,2,3,4 inΩ,uiϵ(x,t)=0,i=1,2,4in ΣT,

where Gj,ϵ=Gj1+ϵ|Gj|,j=1,2,3.

We apply the Faedo-Galerkin method to solve (2.7). Let {el} be a denumerable orthogonal base of H01(Ω) orthonormal with respect to L2(Ω). We consider the sequence of finite dimensional spaces Sn = span {el, ln}. Let uj,nϵ(x,t)=l=1ncj,n,l(t)el(x),j=1,2,3,4, be the weak solution of system (2.7), where cj,n,l(t), j = 1, 2, 3, 4 are unknowns of the nonlinear FODE system,

Ωtαu1,nϵemdx+d1(l(u1,nϵ))Ωu1,nϵemdx                   +ΩG1,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ)emdx=(1ρ)Ωu1,nϵemdx,Ωtαu2,nϵemdx+d2(l(u2,nϵ))Ωu2,nϵemdx                   +ΩG2,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ)emdx                   =Ω(r2u2,nϵ+ρu1,nϵ)emdx,Ωtαu3,nϵemdx+ΩG3,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ,u4,nϵ)emdx                   =Ωr3u3,nϵemdx,Ωtαu4,nϵemdx+d4(l(u4,nϵ))QTu4,nϵemdx                   =Ωξ(u1,nϵ+u2,nϵu4,nϵ)emdx,    (2.8)

for all emSn. Thus, we get

D0Ctαc1,n,m(t)=d1(l(u1,nϵ))Ωu1,nϵemdx            ΩG1,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ)emdx+(1ρ)Ωu1,nϵemdx,                          =:F1m(t,{c1,n,l}l=1n,{c2,n,l}l=1n,{c3,n,l}l=1n,{c4,n,l}l=1n), D0Ctαc2,n,m(t)=d2(l(u2,nϵ))Ωu2,nϵemdx            ΩG2,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ)emdx+Ω(r2u2,nϵ+ρu1,nϵ)ϕ2,ndx,                          =:F2m(t,{c1,n,l}l=1n,{c2,n,l}l=1n,{c3,n,l}l=1n,{c4,n,l}l=1n),D0Ctαc3,n,m(t)=ΩG3,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ,u4,nϵ)emdx     +Ωr3u3,nϵemdx,                            =:F3m(t,{c1,n,l}l=1n,{c2,n,l}l=1n,{c3,n,l}l=1n,{c4,n,l}l=1n),D0Ctαc4,n,m(t)=d4(l(u4,nϵ))Ωu4,nϵemdx+      Ωξ(u1,nϵ+u2,nϵu4,nϵ)emdx,                         =:F4m(t,{c1,n,l}l=1n,{c2,n,l}l=1n,{c3,n,l}l=1n,{c4,n,l}l=1n).    (2.9)

Further, it is east to see that all Fjm,j=1,2,3,4 are continuous functions of cj,n,l(t). From the Theorem 2.1, the system (2.9) has a local solution cj,n,l(t), j = 1, 2, 3, 4 on some interval [0, tn), 0 < tn < T. The extension of these solutions to the whole interval [0, T] is a consequence of the following apriori estimates.

Lemma 2.5. Assume the hypothesis of Theorem 2.3. Then there exists a constant C > 0 is independent on n such that

(u1,nϵ,u2,nϵ,u3,nϵ,u4,nϵ)L(0,T;L2(Ω))C, (u1,nϵ,u2,nϵ,u4,nϵ)L2(0,T;H01(Ω))C,G1,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ)u1,nϵL1(QT)+G2,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ)u2,nϵL1(QT)+G3,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ,u4,nϵ)u3,nϵL1(QT)C.}    (2.10)

Proof. Now, we set ϕj,n(x,t)=l=1nbj,n,l(t)el(x), j=1,2,3,4. the coefficients {b j, n, l}, j = 1, 2, 3, 4. are absolutely continuous functions. Then, from (2.8), the Faedo-Galerkin approximation solution satisfy the following weak formulation

Ωtαu1,nϵϕ1,ndx+d1(l(u1,nϵ))Ωu1,nϵϕ1,ndx         +ΩG1,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ)ϕ1,ndx         =(1-ρ)Ωu1,nϵϕ1,ndx,Ωtαu2,nϵϕ2,ndx+d2(l(u2,nϵ))Ωu2,nϵϕ2,ndx         +ΩG2,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ)ϕ2,ndx            =Ω(r2u2,nϵ+ρu1,nϵ)ϕ2,ndx,Ωtαu3,nϵϕ3,ndx+ΩG3,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ,u4,nϵ)ϕ3,ndx         =Ωr3u3,nϵϕ3,ndx,Ωtαu4,nϵϕ4,ndx+d4(l(u4,nϵ))QTu4,nϵϕ4,ndx         =Ωξ(u1,nϵ+u2,nϵ-u4,nϵ)ϕ4,ndx.    (2.11)

Choosing ϕj,n=uj,nϵ,j=1,2,3,4,, respectively, in the above system (2.11) and summing the resulting terms, we get

12D0CtαΩj=14|uj,nϵ|2dx+Ω(m1|u1,nϵ|2+m2|u2,nϵ|2+m4|u4,nϵ|2)dx     +ΩG1,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ)u1,nϵdx+ΩG2,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ)u2,nϵdx                +ΩG3,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ,u4,nϵ)u3,nϵdxCΩj=14|uj,nϵ|2dx,

where C is the positive constant which does not depend on n and t. Further, using the Lemma 2.3, we have,

Ωj=14|uj,nϵ|2dxC,    (2.12)

for some positive constant C depending only on the given data and independent of n. From (2.12), we get

(u1,nϵ,u2,nϵ,u3,nϵ,u4,nϵ)L(0,T;L2(Ω))C,(u1,nϵ,u2,nϵ,u4,nϵ)L2(0,T;H01(Ω))C.    (2.13)

Using the results (2.12) and (2.13), we get

G1,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ)u1,nϵL1(QT)+G2,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ)u2,nϵL1(QT)+G3,ϵ(x,t,u1,nϵ,u2,nϵ,u3,nϵ,u4,nϵ)u3,nϵL1(QT)C.    (2.14)

Lemma 2.6. Assume the hypothesis of Theorem 2.3. Then {ũ1,nϵ,ũ2,nϵ,ũ4,nϵ} is bounded set of Wα(,H01(Ω),L2(Ω)), where

z˜ϵ,n(x,t)={zϵ,n(x,t)t[0,T],0\[0,T].    (2.15)

Proof. In order to prove the result, Theorem 2.2 with Lemma 2.5, we have to show that

-|τ|2γ|u^j,nϵ(τ)|dτC,j=1,2,4,    (2.16)

for some γ > 0, where ûj,nϵ, j=1,2,4 denote the Fourier transform of ũj,nϵ,j=1,2,4.

Now, rewritten is defined as in (2.7),

(D0Ctαu˜1,nϵ,ϕj)=F˜1n,ϕj+(u1,n(0),ϕj)It1αδ0                                  (u1,n(T),ϕj)It1αδT,(D0Ctαu˜2,nϵ,ϕj)=F˜2n,ϕj+(u2,n(0),ϕj)It1αδ0                                  (u2,n(T),ϕj)It1αδT,(D0Ctαu˜4,nϵ,ϕj)=F˜4n,ϕj+(u4,n(0),ϕj)It1αδ0                                  (u4,n(T),ϕj)It1αδT,}    (2.17)

where F~in is defined as in (2.15),


Here δ0 and δT denote the Dirac distribution at 0 and T. Indeed, it is classical that since ũm has discontinuities at 0 and T, the Caputo derivatives of ũn is given by Zhou and Peng [12] and Zhou et al. [19]. Substitute ϕj,n=ûj,nϵ j=1,2,4,, respectively, in (2.11), and using the Fourier transform, we get

(2iπ τ)α|u^1,nϵ(τ)|2=F^1n,u^1,nϵ+(u1,n(0),u^1,nϵ(τ))(2iπ τ)α1                                                   (u1,n(T),u^1,nϵ(τ))e2iπ Tτ,(2iπ τ)α|u^2,nϵ(τ)|2=F^2n,u^2,nϵ+(u2,n(0),u^2,nϵ(τ))(2iπ τ)α1                                                   (u2,n(T),u^2,nϵ(τ))e2iπ Tτ,(2iπ τ)α|u^4,nϵ(τ)|2=F^4n,u^4,nϵ+(u4,n(0),u^4,nϵ(τ))(2iπ τ)α1                                                    (u4,n(T),u^4,nϵ(τ))e2iπ Tτ.}    (2.18)

It is obvious that the boundedness of solution shows that,

supτF^in(τ)H-1(Ω)C, n, i=1,2,4,

where C > 0 is a positive constant. On account of {(u1,nϵ,u2,nϵ,u4,nϵ)} is bounded set in L(0, T; L2(Ω)), we get

(|u1,n(0)|, |u2,n(0)|, |u4,n(0)|)C and (|u1,n(T)|, |u2,n(T)|, |u4,n(T)|)C,

and from (2.18), we get

|τ|α|u^i,nϵ(τ)|2C max{1,|τ|α-1} u^i,nϵ(τ)|H01(Ω),i=1,2,4.

For γ fixed, γ<α4, we see that



-|τ|2γ|u^i,nϵ(τ)|2dτC(γ)-1+|τ|α1+|τ|α-2γ|u^i,nϵ|2dτ                                                               C1(γ)-u^i,nϵH01(Ω)2dτ                                                               +C2(γ)-|τ|α-1u^i,nϵ(τ)H01(Ω)1+|τ|α-2γdτ.

Applying the Parseval identity, the first integral is bounded as n → ∞, and we have to prove that

-|τ|α-1u^i,nϵ(τ)H01(Ω)1+|τ|α-2γdτC.    (2.19)

From the Schwarz inequality, we prove (2.19) as follows.

-|τ|α-1u^i,nϵ(τ)H01(Ω)1+|τ|α-2γdτ(-dτ(1+|τ|α-2γ)2)12                                                                        (-|τ|2α-2u^i,nϵ(τ)H01(Ω)2dτ),

the first integral is finite due to γ<14. On the other hand, it follows from the Parseval identity that

|τ|2α2u^i,nϵ(τ)H01(Ω)2dτ=(It1αu˜i,nϵ(t)H01(Ω))2dt                                                                                 =0T0It1αui,nϵ(t)H01(Ω)2dt                                                                                 (T1αΓ(2α))20Tui,nϵ(t)H01(Ω)2dt.

From the above integral, we understand that (2.19) is true. Thus, {ũ1,nϵ,ũ2,nϵ,ũ4,nϵ} is bounded set of Wα(,H01(Ω),L2(Ω)).

Theorem 2.4. Suppose the hypotheses of Theorem 2.3 hold true, then the regularized system (2.7) possesses a weak solution (u1ϵ,u2ϵ,u3ϵ,u4ϵ), which satisfies the following conditions:


such that for every ϕjL2(0,T;H01(Ω)),j=1,2,3,4,

0Ttα(u1ϵ,ϕ1)dt+d1(Ωu1ϵdx)QTu1ϵϕ1dxdt      +QTG1,ϵ(x,t,u1ϵ,u2ϵ,u3ϵ)ϕ1dxdt=(1ρ)QTu1ϵϕ1dxdt,0Ttα(u2ϵ,ϕ2)dt+d2(Ωu2ϵdx)QTu2ϵϕ2dxdt      +QTG2,ϵ(x,t,u1ϵ,u2ϵ,u3ϵ)ϕ2dxdt=QT(r2u2ϵ+ρu1ϵ)ϕ2dxdt,0Ttα(u3ϵ,ϕ3)dt+QTG3(x,t,u1ϵ,u2ϵ,u3ϵ,u4ϵ)ϕ3dxdt      =QTr3u3ϵϕ3dxdt,0Ttα(u4ϵ,ϕ4)dt+d4(Ωu4ϵdx)QTu4ϵϕ4dxdt      =QTξ(u1ϵ+u2ϵu4ϵ)ϕ4dxdt.

Proof. By Lemma 2.5, Lemma 2.6 and Theorem 2.2, we can extract the subsequences of (u1,nϵ,u2,nϵ,u3,nϵ,u4,nϵ) such that as n → ∞, we get

(u1,nϵ,u2,nϵ,u3,nϵ,u4,nϵ)(u1ϵ,u2ϵ,u3ϵ,u4ϵ) weak  in L(0,T;L2(Ω)),(u1,nϵ,u2,nϵ,u4,nϵ)(u1ϵ,u2ϵ,u4ϵ) weakly in L2(0,T;H01(Ω)).

It is enough, we show that

dj(l(uj,nϵ))dj(l(ujϵ)) in L2(0,T),T>0,j=1,2,4.

Since dj is continuous functions, it is enough to show that

l(uj,nϵ)l(ujϵ) strongly in L2(0,T).



This result concludes that

dj(l(uj,nϵ))dj(l(ujϵ)) in L2(0,T),T>0,j=1,2,4.

Substitute ϕ1,n = v and integrate the first equation of (2.11) from 0 to t and 0 to t0, respectively. Then, subtract the resulting equation, we get


Using the Lemma 2.6 and the Lebesgue dominated convergence theorem, we can prove that (u1ϵ(t),v)-(u1ϵ(t0),v)0 as tt0 as in Zhou and Peng [12] and Zhou et al. [19]. Similarly, it can be show that (u2ϵ(t),v)-(u2ϵ(t0),v)0 as tt0, (u3ϵ(t),v)-(u3ϵ(t0),v)0 as tt0 and (u4ϵ(t),v)-(u4ϵ(t0),v)0 as tt0.

Now, we prove the Theorem 2.3 and some priori estimates and compactness results.

Lemma 2.7. Assume that the assumptions of Theorem 2.3 are satisfied, then

u1ϵ,u2ϵ,u3ϵ,u4ϵL(0,T;L2(Ω))C,u1ϵ,u2ϵ,u4ϵL2(0,T;H01(Ω))C,G1,ϵ(x,t,u1ϵ,u2ϵ,u3ϵ)u1ϵL1(QT)+G2,ϵ(x,t,u1ϵ,u2ϵ,u3ϵ)u2ϵL1(QT)+G3,ϵ(x,t,u1ϵ,u2ϵ,u3ϵ,u4ϵ)u3ϵL1(QT)C.}    (2.20)

Proof. Let us presume uj-ϵ=sup(0,-ujϵ),j=1,2,3,4. Multiplying the first equation (2.7) by u1-ϵand integrating over Ω, we get


Using the assumptions of the given data, we get


This concludes that the solution u1ϵ is a nonnegative solution. Similarly, we can claim that solution u2ϵ,u3ϵ,u4ϵ are nonnegative. We prove (2.20) as in Lemma 2.5, by replacing ui,nϵ by ujϵ,j=1,2,3,4.

Lemma 2.8. Assume the hypothesis of Theorem 2.3 are hold true. Then {ũ1ϵ,ũ2ϵ,ũ4ϵ} is bounded set of Wα(,H01(Ω),L2(Ω)), where


Proof. As similar as proof of the Lemma 2.6.

Proof of Theorem 2.3: From Lemma 2.7, Lemma 2.8 and Theorem 2.2, we understand that sequences have convergent subsequences which are still denoted by (u1ϵ,u2ϵ,u3ϵ,u4ϵ). Then there exists (u1, u2, u3, u4) as n → ∞, we get

(u1ϵ,u2ϵ,u3ϵ,u4ϵ)(u1,u2,u3,u4) weak * in L(0,T;L2(Ω)),(u1ϵ,u2ϵ,u4ϵ)(u1,u2,u4) weakly in L2(0,T;H01(Ω)).

This concludes the proof of the Theorem 2.3.

Theorem 2.5. The solution (u1, u2, u3, u4) of system (1.1) is unique.

Proof. Let (v1, v2, v3, v4) and (ṽ1, ṽ2, ṽ3, ṽ4) be any two solutions of (1.1). Now, we consider ui = vi − ṽ i, i = 1, 2, 3, 4 and choose ϕi = ui, i = 1, 2, 3, 4 in (2.6), we get


Using the Lipschitz assumptions of di, i = 1, 2, 4, the non-negativity and boundedness of solutions of (2.7) with the Young inequality, we obtain


Summing up the above inequalities and using the Lemma 2.3, we get


This concludes the proof of the theorem.

3. Finite Element Scheme

In this section, we present a finite element scheme for the considered model (1.1). Here, first we present a weak formulation for time-fractional cancer invasion system (1.1), where the time-fractional derivative is the Caputo derivative. Further, the spatial and temporal discretization are presented. Furthermore, an iteration fixed point type is proposed to handle the nonlinear terms of the system as in Ganesan and Shangerganesh [52, 53] and Ganesan and Tobiska [54].

3.1. Finite Element Semi-discretization

Let Th be a partition of Ω into non overlapping triangles TkTh cells, and use the piecewise linear (P1) finite elements on each cell. Now, let VhV be a conforming finite element subspace of V with basis function {ϕk}k=1N such that Vh = span{ϕk}, where N is the number of degrees of freedom. Then the semidiscrete problem, based on (2.6) is to find a solution whVh such that for a.e t ∈ (0, T), for all ψ ∈ Vh

(tαwh,ψ)+B(wh,ψ)=F(wh,ψ),    (3.1)

where B(wh,ψ)=(bu1(u1,h;u1,h;u2,h;u3,h,ψ1)bu2(u2,h;u1,h;u2,h;u3,h,ψ2)bu3(u3,h;u1,h;u2,h;u3,h;u4,h,ψ3)bu4(u4,h;u4,h,ψ4)) & F(wh,ψ)=(0f1(u1,h,ψ2)0f2(u1,h,u2,h,ψ4)).


           bu1(u1,h;u1,h;u2,h;u3,h,ψ1)=Ωd1(l(u1,h))u1,hψ1dx                                                                                             +Ωu1,h(u1,h+β1u2,h                                                                                              +ρ+γ1u3,h1)ψ1dx,           bu2(u2,h;u1,h;u2,h;u3,h,ψ2)=Ωd2(l(u2,h))u2,hψ2dx                                                                                            +Ωu2,h(r2u2,h+β2u1,h                                                                                            +δ1u3,hr2)ψ2dx,  bu3(u3,h;u1,h;u2,h;u3,h;u4,h,ψ3)=Ωu3,h(r3u3,h+γ2u1,h                                                                                            +δ2u2,h+σu4,hr3)ψ3dx,           bu4(u4,h;u1,h;u2,h;u4,h,ψ4)=Ωd4(l(u4,h))u4,hψ4dx                                                                                            +Ωξu4,hψ4dx,                                                          f1(u1,h,ψ2)=Ωρu1,hψ2dx,f2(u1,h,u2,h,ψ4)​​​                                                                                        =Ωξ(u1,h+u2,h)ψ4dx.

Furthermore, using the basis {ϕk}k=1N, we describe the discrete solution whVh, in terms of the basis of Vh as


where uj,k(t), j = 1, 2, 3, 4, t ∈ [0, T], are the unknown coefficients to be determined. Thus, we have

M0CDtαω(t)+Aω=F,    (3.2)

where ω=(u1,1,u2,2,,u1,N,u2,1,u2,2,,u2,N,u3,1,u3,2,,u3,N,u4,1,u4,2,,u4,N) is unknown vector. Further, mass matrix M, stiffness matrices A and known vectors F1, F2 are given by

M=[M0000M0000M0000M],(M)pq=Thϕp(x)ϕq(x),A=[A(U1)A(U2)A(U3)A(U4)],(A(U1))pq=Thd1(l(k=1Nu1,k(t)ϕk(x)))ϕp(x)ϕq(x)dx                             +Th(k=1Nu1,k(t)ϕk(x)+β1k=1Nu2,k(t)ϕk(x)                              +ρ+γ1k=1Nu3,k(t)ϕk(x)1)ϕp(x)ϕq(x)dx,(A(U2))pq=Thd2(l(k=1Nu2,k(t)ϕk(x)))                              ϕp(x)ϕq(x)dx                             +Th(r2k=1Nu2,k(t)ϕk(x)+β2k=1Nu1,k(t)ϕk(x)                              +δ1k=1Nu3,k(t)ϕk(x)r2)ϕp(x)ϕq(x)dx,(A(U3))pq=Th(r3k=1Nu3,k(t)ϕk(x)+                              γ2k=1Nu1,k(t)ϕk(x)+δ2k=1Nu2,k(t)ϕk(x)                             +σk=1Nu4,k(t)ϕk(x)r3)ϕp(x)ϕq(x)dx,(A(U4))pq=Thd1(l(k=1Nu4,k(t)ϕk(x)))                             ϕp(x)ϕq(x)dx,                              F=[0F10F2], (F1)=Thρk=1Nu1,k(t)ϕk(x)ϕq(x)dx,                              (F2)=Thξ(k=1Nu1,k(t)ϕk(x)ϕq(x)                            +k=1Nu2,k(t)ϕk(x)ϕq(x))dx.

The system (3.2) is a FODE system. Solvability of (3.2) can be achieved as in Theorem 2.3 and Theorem 2.5.

3.1.1. Temporal Discretization and Linearization

We present now the time discretization of (3.2). Let 0 = t0 < t1 < t2 < ··· < t𝔑 = T be a decomposition of the considered time interval [0, T] and δt = t𝔯t𝔯−1, 𝔯 = 1, 2, 3, …, 𝔑 represents the uniform time step. We discretize the Caputo fractional time derivative be using finite difference scheme as in Lin and Xu [55] and Sun and Wu [56]. The approximation is given by

tαu(x,t𝔯)=δt-αΓ(2-α)𝔩=0𝔯a𝔩(u(x,t𝔯+1-𝔩)-u(x,t𝔯-𝔩)),    (3.3)

a𝔩=(𝔩+1)α-1-(𝔩)α-1, 𝔩=1,2,3,,𝔯.

Substitute (3.3) in (3.1), we have

(δt-αΓ(2-α)𝔩=0𝔯a𝔩(w𝔯+1-𝔩-w𝔯-𝔩),ϕ)+B(w𝔯+1,ϕ)=F(w𝔯+1,ϕ),    (3.4)

where a𝔩=(𝔩+1)α-1-(𝔩)α-1, 𝔯 = 1, 2, 3, …, 𝔑 and ψ ∈ V. Moreover, we use fixed point iteration technique to control the nonlinear terms of the given system. Initiate with uj,n+10=uj,n,j=1,2,3,4, then the nonlinear terms in (3.4) can be written as

bu1(u1,𝔯+1𝔪;u1,n+1𝔪;u2,𝔯+1𝔪;u3,𝔯+1𝔪,ϕ)bu1                    (u1,𝔯+1𝔪;u1,𝔯+1𝔪-1;u2,𝔯+1𝔪-1;u3,𝔯+1𝔪-1,ϕ),bu2(u2,𝔯+1𝔪;u1,𝔯+1𝔪;u2,𝔯+1𝔪;u3,𝔯+1𝔪,ϕ)bu2                    (u1,𝔯+1𝔪;u1,𝔯+1𝔪-1;u2,𝔯+1𝔪-1;u3,𝔯+1𝔪-1,ϕ),bu3(u3,𝔯+1𝔪;u1,𝔯+1𝔪;u2,𝔯+1𝔪;u3,𝔯+1𝔪;u3,𝔯+1𝔪,ϕ)bu3                    (u3,𝔯+1𝔪;u1,𝔯+1𝔪-1;u2,𝔯+1𝔪-1;u3,𝔯+1𝔪-1;u4,𝔯+1𝔪-1,ϕ),

for 𝔪 = 1, 2, 3, …. In addition, we iterate until the residual is less than the prescribed threshold value (10−10) or the given maximal number of iterations is reached. In computations, the fixed point iteration converges within six or seven iterations for the prescribed residual value. Furthermore, the number of fixed point iteration increases when δt is increased.

4. Numerical Experiment

In this section, we perform series of numerical computations to understand the impact of α in cancer invasion system. Here, all numerical computations are performed in the unit square domain Ω = [0, 1] × [0, 1]. We used Freefem++ [57] for finite element scheme and UMFPACK [58, 59] is used to solve the resulting algebraic system. All computations are carried out by using Intel (R) Core (TM) i7-7700 CPU with 3.60GHz and 8 GB RAM.

4.1. Convergence Study

We consider the cancer invasion model (2.5)

tαu1d1(l(u1))Δu1u1(1u1)+β1u1u2+ρu1+γ1u1u3=fu1,tαu2d2(l(u2))Δu2r2u2(1u2)+β2u1u2ρu1+δ1u2u3=fu2,tαu3r3u3(1u3)+γ2u1u3+δ2u2u3+σu3u4=fu3,tαu4d4(l(u4))Δu4ξ(u1+u2u4)=fu4,}    (4.1)

where fu1, fu2, fu3, and fu4 are forcing terms. They are chosen such that following smooth solutions are satisfies (4.1).


Moreover, we fixed di(l(ui)) = Di sin (l(ui)) where l(u)=Ωudx. Further all other parameters of the model are chosen as

D1=0.0035,D2=0.035,D4=0.0002,β1=0.0015,β2=0.0015, γ1=0.0015,γ2=0.003,r2=0.0012,r3=0.001, ξ=0.1,δ1=0.25,δ2=0.35,ρ=0.0015,σ=0.001.

A set of finite element computations on uniformly refined meshes with δt=h22-α are performed. In order to compare the discretization errors at different mesh levels and verify the order of convergence of numerical scheme, the following errors are computed for each unknowns uj, j = 1, ··· , 4 of the system.


First, for α = 0.4 then the obtained numerical errors and corresponding convergence rates are depicted in Table 1. Then, for α = 0.7 obtained results are shown in Table 2. Tables 1, 2 clearly shows that existence of second order convergence for the errors E1 and E2, respectively, for all unknowns of the system.


Table 1. Errors and order of convergence when α = 0.4.


Table 2. Errors and order of convergence when α = 0.7.

4.2. Numerical Results and Discussion

We understand the influence of fractional derivative on cancer invasion system (1.1) by performing numerical simulations with different values of α. All computations are performed until at end time T = 20. Further, uniform time step size δt is taken as 0.1. We discretize the unit square domain using triangular elements with characteristic element length 140 × 140 and a uniform mesh size h = 0.0101015. We used 19881 degrees of freedom for each unknown in all computations with total 79,524 degrees of freedom. We assumed the homogeneous Dirichlet boundary conditions for all unknowns with the following initial conditions.

u1(x,0)=1.01exp(-(x-0.5)2-(y-0.5)2ϵ1), u2(x,0)=0,u3(x,0)=1-0.99exp(-(x-0.5)2-(y-0.5)2ϵ1),u4(x,0)=1.01exp(-(x-0.5)2-(y-0.5)2ϵ2),

where ϵ1 = 0.005, and ϵ2 = 0.075. We performed simulations for α = 0.1, 0.4, 0.7, 1 and all other parameters are taken of the model of the previous section 4.1.

Now numerical simulations are carried out to analyse the influence of fractional parameter α on the cancer invasion system (2.5). The first two rows of Figure 1 show the comparison between the fractional derivatives when α = 0.1&0.7 and the integer order derivative for cancer density u1 at time T = 10&20. Similarly rows (iii) & (iv) of Figure 1 show the effects on cancer density u2 at time T = 10&20. Differences on the evolution of u1 and u2 can be observed depending on the value of α. We observed huge morphological changes in the invasion of cancer cells with fractional derivatives than the integer order derivative. We note that investigations with different fractional order derivatives suggest that they have relatively little impact on general properties of the cancer invasion system. Fractional derivative equips the distinct sequence of cancer cells (type I and II) migration toward the normal cell domain, see Figure 1 when α = 0.1&0.7. Similar pattern changes also observed in the evolution of normal cells density (u3) and acidification concentration (u4) due to the influence of fractional derivatives, see Figure 2 (i) − (iv).


Figure 1. Evaluation of u1 and u2 at T = 10, 20 with the different values of α = 0.1, 0.7, and 1.


Figure 2. Evaluation of u3 and u4 at T = 10, 20 with the different values of α = 0.1, 0.7, and 1.

Further, the influence of fractional derivatives α = 0.1, 0.4&0.7 compared with α = 1.0 on the evolution of cancer density u1&u2, normal density u3 and acidification concentration u4 are discussed along the y = x. Numerical results are depicted in Figures 36 at time T = 5, 10, 15, and 20. From the these figures (Figures 36), it is clear that density of cancer cells (both type I and II) increasing when α decreases. At the same instance, acidification concentration (due to H+ions) u4 increases, when α decreases. By comparing all these numerical results, we understand that fractional derivatives increase the population of cancer cells at some position of the domain. Therefore, by comparing all the simulation results in Figures 16, it is not difficult to find the that fractional derivatives change the invasion of cancer cells toward normal cells by comparing with integer order derivatives. Therefore, we conclude that proposed computational model can be employed to foresee the location and the shape of the tumor at a particular instance during cancer growth and invasion.


Figure 3. Evaluation of u1, u2, u3, and u4 at t = 5 along the line y = x.


Figure 4. Evaluation of u1, u2, u3, and u4 at t = 10 along the line y = x.


Figure 5. Evaluation of u1, u2, u3, and u4 at t = 15 along the line y = x.


Figure 6. Evaluation of u1, u2, u3, and u4 at t = 20 along the line y = x.

Data Availability

All datasets generated for this study are included in the manuscript and the supplementary files.

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


1. Araujo RP, MceLwain DLS. A history of the study of solid tumour growth: the contribution of mathematical modelling. Bull Math Biol. (2004) 66:1039–91. doi: 10.1016/j.bulm.2003.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Clatz O, Sermesant M, Bondiau PY, Delingette H, Warfield SK, Malandain G, et al. Realistic simulation of the 3D growth of brain tumors in MR images coupling diffusion with biomechanical deformation. IEEE Trans Med Imaging. (2005) 24:1334–46. doi: 10.1109/TMI.2005.857217

CrossRef Full Text | Google Scholar

3. Enderling H, Chaplain MAJ. Mathematical modeling of tumor growth and treatment. Curr Pharm Des. (2014) 20:4934–40. doi: 10.2174/1381612819666131125150434

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Tracqui P. Biophysical models of tumor growth. Rep Prog Phys. (2009) 72:056701. doi: 10.1088/0034-4885/72/5/056701

CrossRef Full Text | Google Scholar

5. Hilfer R. Applications of Fractional Calculus in Physics. Singapore: World Scientific (2000).

6. Kilbas AA, Srivastava HM, Trujillo JJ. Theory and Applications of the Fractional Differential Equations. Amsterdam: Elsevier (2006).

Google Scholar

7. Zhou Y. Basic Theory of Fractional Differential Equations. Singapore: World Scientific (2014).

8. Henry BI, Wearne SL. Fractional reaction-diffusion. Phys A. (2000) 276:448–55. doi: 10.1016/S0378-4371(99)00469-0

CrossRef Full Text | Google Scholar

9. Henry BI, Langlands T, Wearne S. Turing pattern formation in fractional activator-inhibitor systems. Phys Rev E. (2005) 72:026101. doi: 10.1103/PhysRevE.72.026101

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Magin RL. Fractional calculus models of complex dynamics in biological tissues. Comput Math Appl. (2010) 59:1586–93. doi: 10.1016/j.camwa.2009.08.039

CrossRef Full Text | Google Scholar

11. Zhou Y, Peng L. On the time-fractional Navier-Stokes equations. Comput Math Appl. (2017) 73:874–91. doi: 10.1016/j.camwa.2016.03.026

CrossRef Full Text | Google Scholar

12. Zhou Y, Peng L. Weak solutions of the time-fractional Navier-Stokes equations and optimal control. Comput Math Appl. (2017) 73:1016–27. doi: 10.1016/j.camwa.2016.07.007

CrossRef Full Text | Google Scholar

13. Berkowitz B, Scher H. Anomalous transport in random fracture networks. Phys Rev Lett. (1997) 79:4038–41. doi: 10.1103/PhysRevLett.79.4038

CrossRef Full Text | Google Scholar

14. Caputo M. Linear models of dissipation whose q is almost frequency independent. J Geophys. (1967) 13:529–39. doi: 10.1111/j.1365-246X.1967.tb02303.x

CrossRef Full Text | Google Scholar

15. Ahmed E, Hashis AH, Rihan FA. On fractional order cancer model. J Fract Calc Appl. (2012) 3:1–6.

Google Scholar

16. Iyiola OS, Zaman FD. A fractional diffusion equation model for cancer tumor. Am Inst Phys Adv. (2014) 4:107121. doi: 10.1063/1.4898331

CrossRef Full Text | Google Scholar

17. Alikhanov AA. A priori estimates for solutions of boundary value problems for fractional-order equations. Differ Equ. (2010) 46:660–6. doi: 10.1134/S0012266110050058

CrossRef Full Text | Google Scholar

18. Jiao F, Zhou Y. Existence results for fractional boundary value problem via critical point theory. Int J Bifurcation Chaos. (2012) 22:1–17. doi: 10.1142/S0218127412500861

CrossRef Full Text | Google Scholar

19. Zhou Y, Shangerganesh L, Manimaran J, Debbouche A. A class of time fractional reaction-diffusion equation with nonlocal boundary condition. Math Meth Appl Sci. (2018) 41:2987–99. doi: 10.1002/mma.4796

CrossRef Full Text | Google Scholar

20. Hejazi H, Moroney T, Liu F. Stability and convergence of a finite volume method for the space fractional advection-dispersion equation. J Comput Appl Math. (2014) 255:684–97. doi: 10.1016/

CrossRef Full Text | Google Scholar

21. Jia JH, Wang H. A preconditioned fast finite volume scheme for a fractional differential equation discretized on a locally refined composite mesh. J Comput Phys. (2015) 299:842–62. doi: 10.1016/

CrossRef Full Text | Google Scholar

22. Liu F, Zhuang P, Turner I, Burrage K, Anh V. A new fractional finite volume method for solving the fractional diffusion equation. Appl Math Model. (2014) 38:3871–8. doi: 10.1016/j.apm.2013.10.007

CrossRef Full Text | Google Scholar

23. Dehghan M, Abbaszadeh M, Mohebbi A. Analysis of a meshless method for the time fractional diffusion-wave equation. Numer Algor. (2016) 73:445–76. doi: 10.1007/s11075-016-0103-1

CrossRef Full Text | Google Scholar

24. Hosseini VR, Shivanian E, Chen W. Local radial point interpolation (MLRPI) method for solving time fractional diffusion-wave equation with damping. J Comput Phys. (2016) 312:307–32. doi: 10.1016/

CrossRef Full Text | Google Scholar

25. Alikhanov AA. A new difference scheme for the time fractional diffusion equation. J Comput Phys. (2015) 280:424–38. doi: 10.1016/

CrossRef Full Text | Google Scholar

26. Li CP, Zeng FH. Finite difference methods for fractional differential equations. Int J Bifur Chaos. (2012) 22:427–32. doi: 10.1142/S0218127412300145

CrossRef Full Text | Google Scholar

27. Sun H, Sun ZZ, Gao GH. Some temporal second order difference schemes for fractional wave equations, Numer. Methods Partial Differ. Equ. (2016) 32:970–1001. doi: 10.1002/num.22038

CrossRef Full Text | Google Scholar

28. Zhang YN, Sun ZZ, Liao HL. Finite difference methods for the time fractional diffusion equation on non-uniform meshes. J Comput Phys. (2014) 265:195–210. doi: 10.1016/

CrossRef Full Text | Google Scholar

29. Jin B, Lazarov R, Liu Y, Zhou Z. The Galerkin finite element method for a multi-term time-fractional diffusion equation. J Comput Phys. (2015) 281:825–43. doi: 10.1016/

CrossRef Full Text | Google Scholar

30. Jin B, Lazarov R, Pasciak J, Zhou Z. Error analysis of semidiscrete finite element methods for inhomogeneous time-fractional diffusion. IMA J Numer Anal. (2015) 35:561–82. doi: 10.1093/imanum/dru018

CrossRef Full Text | Google Scholar

31. Jin B, Lazarov R, Zhou Z. Error estimates for a semidiscrete finite element method for fractional order parabolic equations. SIAM J Numer Anal. (2013) 51:445–66. doi: 10.1137/120873984

CrossRef Full Text | Google Scholar

32. Le KN, McLean W, Mustapha K. Numerical solution of the time-fractional Fokker-Planck equation with general forcing. SIAM J Numer Anal. (2016) 54:1763–84. doi: 10.1137/15M1031734

CrossRef Full Text | Google Scholar

33. Li XJ, Xu CJ. A space-time spectral method for the time fractional diffusion equation. SIAM J Numer Anal. (2009) 47:2108–31. doi: 10.1137/080718942

CrossRef Full Text | Google Scholar

34. Zheng M, Liu F, Turner I, Anh V. A novel high order space-time spectral method for the time fractional Fokker-Planck equation. SIAM J Scient Comput. (2015) 37:A701–24. doi: 10.1137/140980545

CrossRef Full Text | Google Scholar

35. Zheng M, Liu F, Anh V, Turner I. A high-order spectral method for the multi-term time-fractional diffusion equations. Appl Math Model. (2016) 40:4970–85. doi: 10.1016/j.apm.2015.12.011

CrossRef Full Text | Google Scholar

36. Esen A, Ucar Y, Yagmurlu N, Tasbozan O. A Galerkin finite element method to solve fractional diffusion and fractional diffusion-wave equations. Math Model Anal. (2013) 18:260–73. doi: 10.3846/13926292.2013.783884

CrossRef Full Text | Google Scholar

37. Wang YJ, Liu Y, Li H, Wang JF. Finite element method combined with second-order time discrete scheme for nonlinear fractional Cable equation. Eur Phys J Plus. (2016) 131:61. doi: 10.1140/epjp/i2016-16061-3

CrossRef Full Text | Google Scholar

38. Liu N, Liu Y, Li H, Wang J. Time second-order finite difference/finite element algorithm for nonlinear time-fractional diffusion problem with fourth-order derivative term. Comput Math Appl. (2018) 75:3521–36. doi: 10.1016/j.camwa.2018.02.014

CrossRef Full Text | Google Scholar

39. Jin B, Li B, Zhou Z. An analysis of the Crank Nicolson method for subdiffusion. IMA J Numer Anal. (2017) 38:518–41. doi: 10.1093/imanum/drx019

CrossRef Full Text | Google Scholar

40. Kumar D, Chaudhary S, Kumar VVKS. Galerkin finite element schemes with fractional Crank-Nicolson method for the coupled time-fractional nonlinear diffusion system. Comput Appl Math. (2019) 38:123. doi: 10.1007/s40314-019-0889-2

CrossRef Full Text | Google Scholar

41. Owolabi KM. Mathematical analysis and numerical simulation of patterns in fractional and classical reaction-diffusion systems. Chaos Solit Fract. (2016) 93:89–98. doi: 10.1016/j.chaos.2016.10.005

CrossRef Full Text | Google Scholar

42. Owolabi KM, Atangana A. Numerical approximation of nonlinear fractional parabolic differential equations with Caputo-Fabrizio derivative in Riemann-Liouville sense. Chaos Solit Fract. (2017) 99:171–9. doi: 10.1016/j.chaos.2017.04.008

CrossRef Full Text | Google Scholar

43. Owolabi KM. Robust and adaptive techniques for numerical simulation of non-linear partial differential equations of fractional order. Commun Nonlin Sci Numer Simul. (2017) 44:304–17. doi: 10.1016/j.cnsns.2016.08.021

CrossRef Full Text | Google Scholar

44. Owolabi KM. Mathematical modelling and analysis of two-component system with Caputo fractional derivative order. Chaos Solit Fract. (2017) 103:544–54. doi: 10.1016/j.chaos.2017.07.013

CrossRef Full Text | Google Scholar

45. Owolabi KM, Atangana A. Analysis and application of new fractional Adams Bashforth scheme with Caputo-Fabrizio derivative. Chaos Solit Fract. (2017) 105:111–9. doi: 10.1016/j.csfx.2019.100007

CrossRef Full Text | Google Scholar

46. Solis FJ, Delgadillo SE. Evolution of a mathematical model of an aggressive-invasive cancer under chemotherapy. Comput Math Appl. (2015) 69:545–58. doi: 10.1016/j.camwa.2015.01.013

CrossRef Full Text | Google Scholar

47. Szymańska Z, Morales C, Lachowicz M, Chaplain MAJ. Mathematical modelling of cancer invasion of tissue: the role and effect of nonlocal interactions. Math Models Methods Appl Sci. (2009) 19:257–81. doi: 10.1142/S0218202509003425

CrossRef Full Text | Google Scholar

48. Podlubny I. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications, Vol.198. Elsevier; Academic Press (1999).

Google Scholar

49. Zhou Y, Manimaran J, Shangerganesh L, Debbouche A. Weakness and Mittag-Leffler stability of solutions for time-fractional Keller-Segel models. Int J Nonlin Sci Num. (2018) 19:753–61. doi: 10.1515/ijnsns-2018-0035

CrossRef Full Text | Google Scholar

50. Almeida R. A Gronwall inequality for a general Caputo fractional operator. arXiv: 1705.10079 (2017). doi: 10.7153/mia-2017-20-70

CrossRef Full Text | Google Scholar

51. Temam R. Navier-Stokes Equations: Theory and Numerical Analysis. Amsterdam: North-Holland (1979). doi: 10.1115/1.3424338

CrossRef Full Text | Google Scholar

52. Ganesan S, Shangerganesh L. Galerkin finite element method for cancer invasion mathematical model. Comput Math Appl. (2017) 73:2603–17. doi: 10.1016/j.camwa.2017.04.006

CrossRef Full Text | Google Scholar

53. Ganesan S, Shangerganesh L. A biophysical model for tumor invasion. Commun Nonlin Sci Numer Simul. (2017) 46:135–52. doi: 10.1016/j.cnsns.2016.10.013

CrossRef Full Text | Google Scholar

54. Ganesan S, Tobiska L. An accurate finite element scheme with moving meshes for computing 3d-axisymmetric interface flows. Int J Numer Meth Fluids. (2008) 57:119–138. doi: 10.1002/fld.1624

CrossRef Full Text | Google Scholar

55. Lin Y, Xu C. Finite difference/spectral approximations for the time-fractional diffusion equation. J Comput Phys. (2007) 225:1533–52. doi: 10.1016/

CrossRef Full Text | Google Scholar

56. Sun ZZ, Wu X. A fully discrete scheme for a diffusion wave system. Appl Numer Math. (2006) 56:193–209. doi: 10.1016/j.apnum.2005.03.003

CrossRef Full Text | Google Scholar

57. Hecht F. New development in freefem++. J Numer Math. (2012) 20:251–65. doi: 10.1515/jnum-2012-0013

CrossRef Full Text | Google Scholar

58. Davis TA. Algorithm 832: UMFPACK V4.3-an unsymmetric-pattern multifrontal method. ACM Trans Math Softw. (2004) 30:196–9. doi: 10.1145/992200.992206

CrossRef Full Text | Google Scholar

59. Davis TA. A column pre-ordering strategy for the unsymmetric-pattern multifrontal method. ACM Trans Math Softw. (2004) 30:167–95. doi: 10.1145/992200.992205

CrossRef Full Text

Keywords: cancer invasion dynamic system, fractional differential equations, reaction-diffusion system, weak solution, numerical solution

Citation: Manimaran J, Shangerganesh L, Debbouche A and Antonov V (2019) Numerical Solutions for Time-Fractional Cancer Invasion System With Nonlocal Diffusion. Front. Phys. 7:93. doi: 10.3389/fphy.2019.00093

Received: 25 March 2019; Accepted: 19 June 2019;
Published: 09 July 2019.

Edited by:

Mustafa Inc, Firat University, Turkey

Reviewed by:

Guo-Cheng Wu, Neijiang Normal University, China
Praveen Agarwal, Anand International College of Engineering, India
Carla M. A. Pinto, Instituto Superior de Engenharia do Porto (ISEP), Portugal

Copyright © 2019 Manimaran, Shangerganesh, Debbouche and Antonov. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Amar Debbouche,