Numerical solution of a malignant invasion model using some finite difference methods

Appanah Rao Appadu; Gysbert Nicolaas de Waal

doi:10.1515/dema-2022-0244

Article Open Access

Numerical solution of a malignant invasion model using some finite difference methods

Appanah Rao Appadu and Gysbert Nicolaas de Waal

Published/Copyright: July 28, 2023

Published by

Become an author with De Gruyter Brill

Author Information Explore this Subject

From the journal Demonstratio Mathematica Volume 56 Issue 1

Abstract

In this article, one standard and four nonstandard finite difference methods are used to solve a cross-diffusion malignant invasion model. The model consists of a system of nonlinear coupled partial differential equations (PDEs) subject to specified initial and boundary conditions, and no exact solution is known for this problem. It is difficult to obtain theoretically the stability region of the classical finite difference scheme to solve the set of nonlinear coupled PDEs, this is one of the challenges of this class of method in this work. Three nonstandard methods abbreviated as NSFD1, NSFD2, and NSFD3 are considered from the study of Chapwanya et al., and these methods have been constructed by the use of a more general function replacing the denominator of the discrete derivative and nonlocal approximations of nonlocal terms. It is shown that NSFD1, which preserves positivity when used to solve classical reaction-diffusion equations, does not inherit this property when used for the cross-diffusion system of PDEs. NSFD2 and NSFD3 are obtained by appropriate modifications of NSFD1. NSFD2 is positivity-preserving when the functional relationship [ ψ ( h ) ] 2 = 2 ϕ ( k ) holds, while NSFD3 is unconditionally dynamically consistent with respect to positivity. First, we show that NSFD2 and NSFD3 are not consistent methods. Second, we tried to modify NSFD2 in order to make it consistent but we were not successful. Third, we extend NSFD3 so that it becomes consistent and still preserves positivity. We denote the extended version of NSFD3 as NSFD5. Finally, we compute the numerical rate of convergence in time for NSFD5 and show that it is close to the theoretical value. NSFD5 is consistent under certain conditions on the step sizes and is unconditionally positivity-preserving.

Keywords: standard finite difference method; nonstandard finite difference method; consistency; positivity preserving; cross-diffusion

MSC 2010: 35K55; 65M06; 92-10

1 Introduction

Fundamental understanding of many natural phenomena such as the propagation of heat and sound requires working with partial differential equations (PDEs); some literature can be found in [1,2]. In real world, models consisting of PDEs often possess no exact or analytical solutions; this requires us to obtain numerical approximations to solve and understand such models [3]. A useful and well-known method is the finite difference method (FDM); the method discretises the spatial and time domains, and a solution can be approximated at the discrete points using finite difference approximations [4]. A problem arises when using standard finite difference (SFD) schemes; they do not inherently convey the properties of the exact solution to the numerical solution such as positivity-preserving solutions for certain PDEs [5].

Ronald E. Mickens is the pioneer of nonstandard finite difference (NSFD) methods. He started work on this class of methods around 1990 [6]. NSFD schemes are constructed by using discrete models. The numerical solutions can preserve properties from the exact solution when solving a differential equation; this makes NSFD methods appealing [6]. One advantage of NSFD over classical methods for reaction-diffusion equations is that classical methods can show blow up at long propagation times [7,8]. Studying the stability is not easy for methods discretising nonlinear PDEs. The freezing coefficient technique and von Neumann stability analysis can be used, but freezing coefficient is not an accurate technique. For NSFD methods discretising reaction-diffusion equations, we can obtain conditions for which the methods are positivity-preserving and we can prove boundedness [9]. Positivity and boundedness will ensure stability [9]. The construction of NSFD schemes is based on some fundamental principles [10]:

The denominator of the discrete derivative must be replaced by a more general function, for example, ∂ u ∂ t ≈ u m n + 1 − u m n ϕ ( k ) , where ϕ ( k ) = exp ( k ) − 1 .
In general, we use nonlocal representation of nonlinear terms [11], for example, ( u m ) 2 ≈ u m u m + 1 and ( u m ) 3 ≈ 2 ( u m ) 3 − ( u m ) 2 u m + 1 .
The difference equation should have the same order as the original equation. In general, when the order of the difference equation is larger than the order of the differential equation, spurious solutions might appear as discussed in [12].
The discrete approximation should preserve some important properties of the corresponding differential equation. Properties such as boundedness and positivity should be preserved [13].

A scheme is called a nonstandard method if at least one of the first two principles mentioned earlier is satisfied [14].

Some previous work on advection-diffusion and reaction-diffusion equations can be found in [15–17]. Diffusion equations have been widely applied and studied in the modelling of biological processes such as the spread of diseases and biofilm growth [18,19]. Reaction-diffusion terms are often present in these models to study biological system, these systems involve, most commonly, chemical substances that experience local reaction, and the substances are transformed into each other and spread out over a surface via diffusion [20]. An interesting natural phenomenon due to reaction-diffusion is the formation and spread of patterns such as spots and stripes over the surface of animals through the chemical interaction between cells [20].

In population dynamics, cross-diffusion equations are often present and require positivity-preserving solutions; these equations are strongly coupled nonlinear parabolic systems making them difficult to work with. Cross-diffusion occurs when the concentration gradient of one species induces the flux of another species, briefly discussed in [5,18]. These equations are at the core of modelling several natural processes such as cancer growth [21] and population dynamics [22] via Volterra-Lotka cross-diffusion systems and chemotaxis [23].

Mathematical analysis for cross-diffusion equations is a challenge, which is largely underdeveloped [18]. From a theoretical point of view, cross-diffusion equations are challenging mainly because they are strongly coupled nonlinear parabolic systems, which do not enjoy the max principle, and thus, deriving appropriate estimates and proving the existence of positive solutions is not easy [18]. However, some results on global and local existence of solutions as well as on their long-time behaviour have been established in [24,25].

We consider a diffusion matrix of a system consisting of two species where both species concentration gradients induce a flux on each other and themselves:

D = D 11 D 12 D 21 D 22 ,

where D m n are the cross-diffusion coefficients [26], and D m n is the diffusion of species m related to the concentration gradient of species n . Unlike reaction-diffusion coefficients, the cross-diffusion coefficients can be negative when species m increases the gradient concentration of species n [26]. The diffusion matrix is not strictly diagonal and not symmetric positive, which distinguishes it from reaction-diffusion systems [18].

The design of reliable numerical methods that produce non-negative solutions for cross-diffusion equations has been an open problem until some interesting results from Chapwanya et al. (see [18]).

The article is organised as follows. In Section 2, we describe the malignant invasion model to be solved and give some discussion on the numerical rate of convergence. Section 3 is dedicated to the classical scheme SFD1 to solve the malignant invasion model. In Section 3.2, we present the numerical results using SFD1 for the malignant invasion model. In Section 4, we describe NSFD1 scheme for the malignant invasion model, and we also check the consistency. We present some numerical results in Section 4.1. We also consider another case of initial conditions in Section 4.2 and present some results to show that NSFD1 does not always preserve positivity. In Sections 5 and 6, we describe NSFD2 and NSFD3 schemes for the malignant invasion model and check the consistency for both schemes. We also present some numerical results for both schemes in Sections 5.2 and 6.2. In Section 7, we construct a new method termed as NSFD4 by modifying NSFD2 and check the consistency. In Section 8, we modify NSFD3 in order to obtain NSFD5 and present some numerical results in Sections 8.2 and 8.3. Finally, in Section 9, we give some concluding remarks.

2 Malignant invasion model

The mathematical model for cancer growth proposed in [21] is described as follows:

(1) ∂ u ∂ t = u ( 1 − u ) − ∂ ∂ x u ∂ c ∂ x ,

(2) ∂ c ∂ t = − p c ,

(3) ∂ p ∂ t = ε − 1 ( u c − p ) ,

where u ( x , t ) , c ( x , t ) , and p ( x , t ) describe the concentration level of invasive cells, connective tissue, and protease, respectively.

The system considered is a cross-diffusion system, since in equation (1), the concentration gradient of the invasive cells induces a flux on the concentration gradient of connective tissue. The parameter ε > 0 chosen to be small is a parameter that signifies the fact that protease enzymes (proteins) are much smaller than those of connective tissue and invasive cells.

We consider two cases with initial conditions described by [18]:

(4) ( i ) u ( x , 0 ) = exp ( − x 2 ) ; c ( x , 0 ) = 1 − 0.5 exp ( − x 2 ) ; p ( x , 0 ) = 0.5 exp ( − x 2 )

(5) ( ii ) u ( x , 0 ) = ( x − 2 ) 2 1 + ( x − 2 ) 2 ; c ( x , 0 ) = 1 1 + x 2 ; p ( x , 0 ) = 0.5 u ( x , 0 ) .

For both cases, the space domain is x ∈ [ 0 , 20 ] and time domain is t ∈ [ 0 , 50 ] .

The second case of initial conditions is used to computationally show that NSFD1 does not always replicate the positivity property of solutions of the continuous model.

For appropriate boundary conditions, the zero-flux boundary condition can be considered at the left boundary and Dirichlet’s or zero-flux boundary conditions can be considered at the right boundary [27]. We will consider zero-flux conditions at both the right boundary and the left boundary.

As discussed in [18], the solutions have the property:

u ( x , t ) ≥ 0 , c ( x , t ) ≥ 0 , p ( x , t ) ≥ 0 ,

where c decreases in time whenever the initial conditions are non-negative. It follows that numerical solution should also be non-negative.

Taking the right-hand side to be zero, the system described by equations (1), (2), and (3) has three types of constant steady-state solutions E = ( u , c , p ) :

E t = ( 0 , 0 , 0 ) trivial equilibrium,

E m = ( 1 , 0 , 0 ) fully malignant equilibrium,

E n = ( 0 , c , 0 ) normal healthy equilibrium, where c > 0 is any constant.

The solution domain is ( x , t ) ∈ [ 0 , x max ] × [ 0 , T ] , partitioning the interval x 0 < x 1 < ⋯ < x N with x m = m h for m = 0 , 1 , 2 , … , N gives the space step size h = x max ∕ N , equally dividing the spatial domain [ 0 , x max ] . Similarly, the time step size is given by k = T ∕ M , equally dividing the temporal domain [ 0 , T ] , with t n = n k for n = 0 , 1 , 2 , … , M [5].

As we do not have an exact solution, we obtain the numerical rate of convergence in time by using

(6) R T = ln E k E k 2 ln ( 2 ) ,

where E k = ∣ ∣ U k − U 2 k ∣ ∣ and E k 2 = ∣ ∣ U k 2 − U k ∣ ∣ are discrete maximum norm errors [16,17].

All numerical simulations are done in MATLAB using an Intel Core i5-10600k with 16GB RAM.

3 Standard method to solve the malignant invasion model

3.1 Derivation of SFD1

We first consider the standard finite difference scheme constructed in [5], which we denote as SFD1. Khalsaraei et al. [5] discretised equations (1), (2), and (3) as follows:

(7) U m n + 1 − U m n k = U m n ( 1 − U m n ) − U m n − U m − 1 n h C m n − C m − 1 n h − U m n C m + 1 n − 2 C m n + C m − 1 n h 2 ,

(8) C m n + 1 − C m n k = − P m n C m n + 1 ,

(9) P m n + 1 − P m n k = ε − 1 ( U m n C m n + 1 − P m n ) .

Rewriting equations (7), (8), and (9), they obtained the explicit form as [5]:

(10) U m n + 1 = ( 1 + k ( 1 − U m n ) ) U m n − k h 2 ( U m n ( C m + 1 n − C m n ) + U m − 1 n ( C m − 1 n − C m n ) ) ,

(11) C m n + 1 = C m n ( 1 + k P m n ) ,

and

(12) P m n + 1 = P m n + ε − 1 k ( U m n C m n + 1 − P m n ) .

We would like to mention that it appears that equations (7) and (9) from [5] have typing errors which we corrected in order to obtain equations (10) and (12) as given in [5].

To check for consistency and to obtain the order of accuracy of SFD1, we obtain the Taylor series expansion of equations (10)–(12) about ( t n , x m ) . Dividing throughout by k and after some re-arrangement, we obtain

(13) ∂ U ∂ t − U ( 1 − U ) + ∂ U ∂ x ∂ C ∂ x + U ∂ 2 C ∂ x 2 = − k 2 ∂ 2 U ∂ t 2 − k 2 6 ∂ 3 U ∂ t 3 − k 3 24 ∂ 4 U ∂ t 4 + h 2 ∂ U ∂ x ∂ 2 C ∂ x 2 − h 2 6 ∂ U ∂ x ∂ 3 C ∂ x 3 + h 3 12 ∂ 2 U ∂ x 2 ∂ 3 C ∂ x 3 + O ( k 4 ) + O ( h 4 ) ,

(14) ∂ C ∂ t + P C = − k 2 ∂ 2 C ∂ t 2 − k 2 6 ∂ 3 C ∂ t 3 − k 3 24 ∂ 4 C ∂ t 4 − k ∂ C ∂ t P − k 2 2 ∂ 2 C ∂ t 2 P − k 3 6 ∂ 3 C ∂ t 3 P + O ( k 4 ) ,

and

(15) ∂ P ∂ t − ε − 1 U C + k U ∂ C ∂ t + k 2 2 U ∂ 2 C ∂ t 2 + k 3 6 U ∂ 3 C ∂ t 3 − P = − k 2 ∂ 2 P ∂ t 2 − k 2 6 ∂ 3 P ∂ t 3 − k 3 24 ∂ 4 P ∂ t 4 + O ( k 4 ) .

As k , h → 0 , we have

∂ U ∂ t − U ( 1 − U ) + ∂ U ∂ x ∂ C ∂ x + U ∂ 2 C ∂ x 2 = 0 ,

∂ C ∂ t + P C = 0 ,

and

∂ P ∂ t − ε − 1 ( U C − P ) = 0 .

Hence, SFD1 is consistent with the PDEs. The scheme is first-order accurate both in time and in space.

3.2 Numerical results using SFD1

Stability analysis using Von Neumann conditions cannot be used to obtain the stability region as the equations are coupled. So we fix h = 1.0 , and by running some experiments with different values of k , we try to obtain a range of values of k such that reasonable numerical solutions are obtained. This is not an ideal method but we do not think there is an alternative method. We therefore fix h = 1.0 and decrease k 50 100 , 50 101 , 50 102 , … , until we obtain reasonable numerical solutions. Reasonable solutions here refer to solutions that are bounded with possibly some dispersive oscillations at some values of x and t . Standard FDMs generally generate unphysical solutions for some PDEs, especially with advection and reactive terms present.

Figure 1 shows results when h = 1.0 and k = 0.5 , and we observe unbounded solutions. Figure 2 shows results when h = 1.0 and k = 50 143 , and we observe that the profiles are not very smooth, as clearly seen in Figure 2(e). Reasonable solutions are obtained when h = 1.0 with k = 50 166 or k = 50 200 , as displayed in Figures 3 and 4.

Figure 1

Results for cancer growth model using SFD1 with k = 0.5 and h = 1.0 . (a) Plot of numerical solutions vs t at x = 5 , (b) plot of numerical solution for conc. of IC vs x vs t , (c) plot of numerical solution for conc. of CT vs x vs t , and (d) plot of numerical solution for conc. of protease vs x vs t .

Figure 2

Results for cancer growth model using SFD1 with k = 0.3497 and h = 1.0 . (a) Plot of numerical solutions vs t at x = 5 , (b) plot of numerical solutions for conc. of IC vs x vs t , (c) plot of numerical solutions for conc. of CT vs x vs t , (d) plot of numerical solutions for conc. of protease vs x vs t and (e) plot of numerical solutions for conc. of protease vs t at x = 1 .

Figure 3

Results for cancer growth model using SFD1 with k = 0.3012 and h = 1.0 . (a) Plot of numerical solutions vs t at x = 5 , (b) plot of numerical solution for conc. of IC vs x vs t , (c) plot of numerical solution for conc. of CT vs x vs t , and (d) plot of numerical solution for conc. of protease vs x vs t .

Figure 4

Results for cancer growth model using SFD1 with k = 0.25 and h = 1.0 . (a) Plot of numerical solutions vs t at x = 5 , (b) plot of numerical solution for conc. of IC vs x vs t , (c) plot of numerical solution for conc. of CT vs x vs t , and (d) plot of numerical solution for conc. of protease vs x vs t .

4 NSFD1 scheme to solve the malignant invasion model

The logistic equation is the space-independent case of equation (1), described by:

(16) d u d t = u ( 1 − u ) .

An exact scheme constructed in [13] to solve equation (16) is

(17) U n + 1 − U n ϕ ( k ) = U n ( 1 − U n + 1 ) ,

where ϕ ( k ) = exp ( k ) − 1 , with k → 0 .

The initial value problem of equation (16) is a solvable separable differential equation, and it can be seen from the solution that in using the exponential denominator function ϕ ( k ) , there is no error between the continuous and discrete solutions of the initial value problem of equation (16), unlike the case when using the denominator k [11,18].

Using the nonstandard approach, Anguelov et al. [15] discretised the Fisher-Kolmogorov-Petrovsky-Piskunov (Fisher-KPP) equation, and using these ideas, equation (1) was discretised by Chapwanya et al. [18] as follows:

(18) U m n + 1 − U m n ϕ ( k ) = U m n ( 1 − U m n + 1 ) − U m n − U m − 1 n ψ ( h ) C m n − C m − 1 n ψ ( h ) − U m n C m + 1 n − 2 C m n + C m − 1 n [ ψ ( h ) ] 2 ,

where ϕ ( k ) = exp ( k ) − 1 and ψ ( h ) = h . Equation (18) can be rewritten as:

(19) U m n + 1 − U m n ϕ ( k ) = U m n ( 1 − U m n + 1 ) − U m n C m + 1 n − ( U m n + U m − 1 n ) C m n + U m − 1 n C m − 1 n [ ψ ( h ) ] 2 .

Chapwanya et al. [18] chose the functional relation

(20) ϕ ( k ) [ ψ ( h ) ] 2 = 1 2 .

This gives

U m n + 1 = U m n + ϕ ( k ) U m n ( 1 − U m n + 1 ) − 1 2 [ U m n C m + 1 n − ( U m n + U m − 1 n ) C m n + U m − 1 n C m − 1 n ] .

On re-arrangement, the following explicit NSFD method was obtained in [18], namely

(21) U m n + 1 = 2 U m n ( 1 + ϕ ( k ) ) − ( U m n C m + 1 n − ( U m n + U m − 1 n ) C m n + U m − 1 n C m − 1 n ) 2 ( 1 + ϕ ( k ) U m n ) .

Equation (21) can be rewritten as:

(22) U m n + 1 = U m n ( 2 + 2 ϕ ( k ) + C m n − C m + 1 n ) + U m − 1 n ( C m n − C m − 1 n ) 2 ( 1 + ϕ ( k ) U m n ) .

From equation (22), we see that preservation of the positivity of solutions for invasive cells is guaranteed if C m n ≥ C m − 1 n , which is a very strict restriction. This means that there is a restriction on the type of profile at a given time and on the initial profile. The complication arises from equation (1) with the addition of the cross-diffusion term and will require appropriate manipulation to obtain a scheme that enjoys positivity-preserving solutions [18].

Using the usual rules of the nonstandard approach, namely the non-local approximation of non-linear terms and the use of complex functions as denominators, the following discretisations are used to solve the PDEs described by equations (2) and (3) [18]. We have

(23) C m n + 1 − C m n ϕ ( k ) = − P m n C m n + 1 ,

where ϕ ( k ) = exp ( k ) − 1 , and making C m n + 1 the subject of the equation, we obtain

(24) C m n + 1 = C m n 1 + ϕ ( k ) P m n .

Also,

(25) P m n + 1 − P m n ε ϕ ( ε − 1 k ) = ε − 1 ( U m n C m n + 1 − P m n + 1 ) ,

which gives

(26) P m n + 1 = P m n + ϕ ( ε − 1 k ) U m n C m n + 1 1 + ϕ ( ε − 1 k ) .

In an extreme situation with initial values between 0 and 1, the maximum value of U m n + 1 is

1 ( 2 + 2 ϕ ( k ) + 1 − 0 ) + 1 ( 1 ) 2 ( 1 + ϕ ( k ) ( 0 ) ) = ( 2 + ϕ ( k ) ) ,

while the minimum value of U m n + 1 is

− 1 2 .

Hence, the numerical solutions from NSFD1 might not always be bounded between 0 and 1.

We next check if NSFD1 is consistent. Taking the Taylor series expansion of equation (22) and after some re-arrangement using k ϕ ( k ) ≈ 1 and 1 2 ϕ ( k ) = 1 h 2 , we obtain

(27) ∂ U ∂ t − U ( 1 − U ) + ∂ U ∂ x ∂ C ∂ x + U ∂ 2 C ∂ x 2 = − k 2 ∂ 2 U ∂ t 2 + 2 U ∂ U ∂ t − k 2 6 ∂ 3 U ∂ t 3 + 3 U ∂ 2 U ∂ t 2 − k 3 24 ∂ 4 U ∂ t 4 + 4 U ∂ 3 U ∂ t 3 + h 2 ∂ U ∂ x ∂ 2 C ∂ x 2 + ∂ 2 U ∂ x 2 ∂ C ∂ x − h 2 2 1 3 ∂ U ∂ x ∂ 3 C ∂ x 3 + 1 2 ∂ 2 U ∂ x 2 ∂ 2 C ∂ x 2 + 1 3 ∂ 3 U ∂ x 3 ∂ C ∂ x + 1 6 U ∂ 4 C ∂ x 4 + h 3 12 1 2 ∂ U ∂ x ∂ 4 C ∂ x 4 + ∂ 2 U ∂ x 2 ∂ 3 C ∂ x 3 + ∂ 3 U ∂ x 3 ∂ 2 C ∂ x 2 + 1 2 ∂ 4 U ∂ x 4 ∂ C ∂ x + O ( k 4 ) + O ( h 4 ) .

Similarly, the Taylor series expansions of equations (24) and (26) about ( t n , x m ) give

(28) ∂ C ∂ t + P C = − k 2 ∂ 2 C ∂ t 2 + 2 ∂ C ∂ t P − k 2 6 ∂ 3 C ∂ t 3 + 3 ∂ 2 C ∂ t 2 P − k 3 24 ∂ 4 C ∂ t 4 + 4 ∂ 3 C ∂ t 3 P − k 4 24 ∂ 4 C ∂ t 4 + O ( k 5 )

and

(29) ∂ P ∂ t − ε − 1 ( U C − P ) = k 2 2 ε − 1 U ∂ C ∂ t − ∂ P ∂ t − ∂ 2 P ∂ t 2 + k 2 6 3 ε − 1 U ∂ 2 C ∂ t 2 − ∂ 2 P ∂ t 2 − ∂ 3 P ∂ t 3 + k 3 24 4 ε − 1 U ∂ 3 C ∂ t 3 − ∂ 3 P ∂ t 3 − ∂ 4 P ∂ t 4 + k 4 24 ε − 1 U ∂ 4 C ∂ t 4 − ∂ 4 P ∂ t 4 + O ( k 5 ) .

As k , h → 0 , we have

∂ U ∂ t − U ( 1 − U ) + ∂ U ∂ x ∂ C ∂ x + U ∂ 2 C ∂ x 2 = 0 ,

∂ C ∂ t + P C = 0 ,

and

∂ P ∂ t − ε − 1 ( U C − P ) = 0 .

Hence, NSFD1 scheme is consistent with the PDEs. The scheme is first-order accurate both in time and in space.

4.1 Numerical results for NSFD1 (case 1)

In Figures 5 and 6, we display the plot of the numerical solutions for concentration of IC, CT, and protease using NSFD1 with k = 1.0 , h = 1.8538 , and k = 0.4055 , h = 1.0 , respectively.

$Figure 5 Results for cancer growth model using NSFD1 with k = 1.0 k=1.0 and h = 1.8538 h=1.8538 (case 1). (a) Plot of numerical solutions vs t t at x ( 4 ) = 5.5614 x\left(4)=5.5614 , (b) plot of numerical solution for conc. of IC vs x x vs t t , (c) plot of numerical solution for conc. of CT vs x x vs t t , and (d) plot of numerical solution for conc. of protease vs x x vs t t .$

Figure 5

Results for cancer growth model using NSFD1 with k = 1.0 and h = 1.8538 (case 1). (a) Plot of numerical solutions vs t at x ( 4 ) = 5.5614 , (b) plot of numerical solution for conc. of IC vs x vs t , (c) plot of numerical solution for conc. of CT vs x vs t , and (d) plot of numerical solution for conc. of protease vs x vs t .

$Figure 6 Results for cancer growth model using NSFD1 with k = 0.4055 k=0.4055 and h = 1.0 h=1.0 (case 1). (a) Plot of numerical solutions vs t t at x ( 6 ) = 5 x\left(6)=5 , (b) plot of numerical solution for conc. of IC vs x x vs t t , (c) plot of numerical solution for conc. of CT vs x x vs t t , and (d) plot of numerical solution for conc. of protease vs x x vs t t .$

Figure 6

Results for cancer growth model using NSFD1 with k = 0.4055 and h = 1.0 (case 1). (a) Plot of numerical solutions vs t at x ( 6 ) = 5 , (b) plot of numerical solution for conc. of IC vs x vs t , (c) plot of numerical solution for conc. of CT vs x vs t , and (d) plot of numerical solution for conc. of protease vs x vs t .

Figure 7 gives the plot of the numerical solutions for concentration of IC, CT, and protease vs x at some values of t using NSFD1 with k = 1.0 and h = 1.8538 .

Figure 7

Plot of numerical solutions vs x at some values of t using NSFD1 with k = 1.0 , h = 1.8538 (case 1). (a) Plot of numerical solutions vs t at x = 5 , (b) plot of numerical solutions vs x at t = 15 , (c) plot of numerical solutions vs x at t = 25 , and (d) plot of numerical solutions vs x at t = 35 .

From Figures 5 to 7, we observe that the numerical solutions using NSFD1 are non-negative and bounded between 0 and 1 for case 1.

4.2 Numerical results using NSFD1 (case 2)

To computationally show that the scheme in equation (22) is not always positivity-preserving, we consider the case of initial conditions given by equation (5). Clearly, C m 0 ≥ C m − 1 0 is not satisfied, as shown in Figures 8(c), 9(c), and 10(a).

$Figure 8 Results for cancer growth model using NSFD1 with k = 0.4055 k=0.4055 and h = 1.0 h=1.0 (case 2). (a) Plot of numerical solutions vs t t at x ( 6 ) = 5 x\left(6)=5 , (b) plot of numerical solution for conc. of IC vs x x vs t t , (c) plot of numerical solution for conc. of CT vs x x vs t t , and (d) plot of numerical solution for conc. of protease vs x x vs t t .$

Figure 8

Results for cancer growth model using NSFD1 with k = 0.4055 and h = 1.0 (case 2). (a) Plot of numerical solutions vs t at x ( 6 ) = 5 , (b) plot of numerical solution for conc. of IC vs x vs t , (c) plot of numerical solution for conc. of CT vs x vs t , and (d) plot of numerical solution for conc. of protease vs x vs t .

$Figure 9 Results for cancer growth model using NSFD1 with k = 1.0 k=1.0 and h = 1.8538 h=1.8538 (case 2). (a) Plot of numerical solutions vs t t at x ( 4 ) = 5.5614 x\left(4)=5.5614 , (b) plot of numerical solution for conc. of IC vs x x vs t t , (c) plot of numerical solution for conc. of CT vs x x vs t t , and (d) plot of numerical solution for conc. of protease vs x x vs t t .$

Figure 9

Results for cancer growth model using NSFD1 with k = 1.0 and h = 1.8538 (case 2). (a) Plot of numerical solutions vs t at x ( 4 ) = 5.5614 , (b) plot of numerical solution for conc. of IC vs x vs t , (c) plot of numerical solution for conc. of CT vs x vs t , and (d) plot of numerical solution for conc. of protease vs x vs t .

Figure 10

Plot of numerical solutions vs x at some values of t using NSFD1 with k = 1.0 and h = 1.8538 (case 2). (a) Plot of numerical solutions vs x at t = 0 , (b) plot of numerical solutions vs x at t = 4 , (c) plot of numerical solutions vs x at t = 8 , and (d) plot of numerical solutions vs x at t = 12 .

Figures 8, 9(b) and (d) illustrate that NSFD1 does not always preserve the positivity of solutions.

Figure 10 shows that the numerical solution for the concentration of invasive cells converges to 1 as we advance in time and space. In Figure 10(b), we observe that the numerical solutions for the concentration of invasive cells and protease are larger than 1.

5 Positivity-preserving NSFD methods

5.1 NSFD2 to solve the malignant invasion model

Chapwanya et al. [18] proposed a modification by replacing U m n on the left-hand side of equation (19) by U m n + U m − 1 n 2 [18]. This gives

(30) U m n + 1 − U m n + U m − 1 n 2 ϕ ( k ) = U m n ( 1 − U m n + 1 ) − [ U m n C m + 1 n − ( U m n + U m − 1 n ) C m n + U m − 1 n C m − 1 n ] [ ψ ( h ) ] 2 ,

where ϕ ( k ) = exp ( k ) − 1 and ψ ( h ) = h .

Choosing ϕ ( k ) [ ψ ( h ) ] 2 = 1 2 and upon some re-arrangement give the first equation for NSFD2 as:

(31) U m n + 1 = U m n ( 2 ϕ ( k ) + 1 + C m n − C m + 1 n ) + U m − 1 n ( 1 − C m − 1 n + C m n ) 2 ( 1 + ϕ ( k ) U m n ) .

For NSFD1 and NSFD2, equations (2) and (3) are discretised in a similar manner, i.e.

C m n + 1 = C m n 1 + ϕ ( k ) P m n

and

P m n + 1 = P m n + ϕ ( ε − 1 k ) U m n C m n + 1 1 + ϕ ( ε − 1 k ) ,

where ϕ ( k ) = exp ( k ) − 1 .

Theorem

[18] For U m 0 ≥ 0 , C m 0 ≥ 0 , and P m 0 ≥ 0 , we have C m n ≥ 0 and P m n ≥ 0 , with the sequence { C m n } being decreasing in n. Moreover if C m 0 ≤ 1 , then U m n ≥ 0 . NSFD2 scheme when used to discretise equation (1) is positivity-preserving if C m 0 ≤ 1 .

To check for consistency of NSFD2 scheme, we consider equation (31). From the Taylor series expansion and after some re-arrangement using k ϕ ( k ) ≈ 1 and 1 2 ϕ ( k ) = 1 h 2 , we obtain

(32) ∂ U ∂ t − U ( 1 − U ) + ∂ U ∂ x ∂ C ∂ x + U ∂ 2 C ∂ x 2 = − 1 h ∂ U ∂ x + 1 2 ∂ 2 U ∂ x 2 − h 6 ∂ 3 U ∂ x 3 + h 2 24 ∂ 4 U ∂ x 4 − k 2 ∂ 2 U ∂ t 2 + 2 U ∂ U ∂ t − k 2 6 ∂ 3 U ∂ t 3 + 3 U ∂ 2 U ∂ t 2 − k 3 24 ∂ 4 U ∂ t 4 + 4 U ∂ 3 U ∂ t 3 + h 2 ∂ U ∂ x ∂ 2 C ∂ x 2 + ∂ 2 U ∂ x 2 ∂ C ∂ x − h 2 2 1 3 ∂ U ∂ x ∂ 3 C ∂ x 3 + 1 2 ∂ 2 U ∂ x 2 ∂ 2 C ∂ x 2 + 1 3 ∂ 3 U ∂ x 3 ∂ C ∂ x + 1 6 U ∂ 4 C ∂ x 4 + h 3 12 1 2 ∂ U ∂ x ∂ 4 C ∂ x 4 + ∂ 2 U ∂ x 2 ∂ 3 C ∂ x 3 + ∂ 3 U ∂ x 3 ∂ 2 C ∂ x 2 + 1 2 ∂ 4 U ∂ x 4 ∂ C ∂ x + O ( k 4 ) + O ( h 4 ) .

Hence, as k , h → 0 , we have

(33) ∂ U ∂ t − U ( 1 − U ) + ∂ U ∂ x ∂ C ∂ x + U ∂ 2 C ∂ x 2 = − 1 h ∂ U ∂ x + 1 2 ∂ 2 U ∂ x 2 .

From equation (33), we conclude that NSFD2 scheme is not consistent with the PDE for invasive cell concentration. Hence, NSFD2 is not a reliable method to solve the malignant invasion model. We note that this is a novel result obtained; the consistency of NSFD2 was not checked in [18].

5.2 Numerical results using NSFD2

Figures 11 and 12 display the plot of the numerical solutions for concentration of IC, CT, and protease using NSFD2 with k = 1.0 , h = 1.8538 , and k = 0.4055 , h = 1.0 , respectively.

$Figure 11 Results for cancer growth model using NSFD2 with k = 1.0 k=1.0 and h = 1.8538 h=1.8538 (case 1). (a) Plot of numerical solutions vs t t at x ( 4 ) = 5.5614 x\left(4)=5.5614 , (b) plot of numerical solution for conc. of IC vs x x vs t t , (c) plot of numerical solution for conc. of CT vs x x vs t t , and (d) plot of numerical solution for conc. of protease vs x x vs t t .$

Figure 11

Results for cancer growth model using NSFD2 with k = 1.0 and h = 1.8538 (case 1). (a) Plot of numerical solutions vs t at x ( 4 ) = 5.5614 , (b) plot of numerical solution for conc. of IC vs x vs t , (c) plot of numerical solution for conc. of CT vs x vs t , and (d) plot of numerical solution for conc. of protease vs x vs t .

Figure 12

Results for cancer growth model using NSFD2 with k = 0.4055 and h = 1.0 (case 1). (a) Plot of numerical solutions vs t at x = 5 , (b) plot of numerical solution for conc. of IC vs x vs t , (c) plot of numerical solution for conc. of CT vs x vs t , and (d) plot of numerical solution for conc. of protease vs x vs t .

Figures 11 and 12 show that numerical solutions obtained using NSFD2 are non-negative.

The results using NSFD2 are very different from those when using NSFD1. This is expected as NSFD2 is not a consistent method.

We observe that the initial change in concentration gradients occurs rather quickly, as seen in Figures 11 and 12. This initial change should occur at t = ± 5 , as in Figures 3–6.

6 NSFD3 to solve the malignant invasion model

6.1 Construction and analysis of NSFD3

By decomposing the cross-diffusion term in equation (19) into its positive and negative parts and multiplying the negative part by U m n + 1 U m n , another scheme is constructed in [18] as follows:

(34) U m n + 1 − U m n ϕ ( k ) = U m n ( 1 − U m n + 1 ) + ( U m n + U m − 1 n ) C m n [ ψ ( h ) ] 2 − ( U m n C m + 1 n + U m − 1 n C m − 1 n ) [ ψ ( h ) ] 2 × U m n + 1 U m n ,

where ϕ ( k ) = exp ( k ) − 1 and ψ ( h ) = h .

By taking R = ϕ ( k ) [ ψ ( h ) ] 2 > 0 and making U m n + 1 the subject of the equation, we obtain the explicit NSFD scheme, which is denoted by NSFD3:

(35) U m n + 1 = ( 1 + ϕ ( k ) ) U m n + R ( U m n + U m − 1 n ) C m n 1 + ϕ ( k ) U m n + R ( U m n C m + 1 n + U m − 1 n C m − 1 n ) U m n .

Theorem

The scheme given by equation (35) is unconditionally dynamically consistent with respect to positivity [18].

Next, we check for consistency of NSFD3, and therefore, consider equation (35). From the Taylor series expansion and after some re-arrangement using k ϕ ( k ) ≈ 1 , we obtain

(36) ∂ U ∂ t − U ( 1 − U ) + ∂ U ∂ x ∂ C ∂ x + U ∂ 2 C ∂ x 2 = − k 2 ∂ 2 U ∂ t 2 + 2 U ∂ U ∂ t + 2 y h 2 U ∂ U ∂ t − k 2 6 ∂ 3 U ∂ t 3 + 3 U ∂ 2 U ∂ t 2 + 3 y h 2 U ∂ 2 U ∂ t 2 − k 3 24 ∂ 4 U ∂ t 4 + 4 U ∂ 3 U ∂ t 3 + 4 y h 2 U ∂ 3 U ∂ t 3 + h 2 ∂ U ∂ x ∂ 2 C ∂ x 2 + ∂ 2 U ∂ x 2 ∂ C ∂ x − h 2 2 1 3 ∂ U ∂ x ∂ 3 C ∂ x 3 + 1 2 ∂ 2 U ∂ x 2 ∂ 2 C ∂ x 2 + 1 3 ∂ 3 U ∂ x 3 ∂ C ∂ x + 1 6 U ∂ 4 C ∂ x 4 + h 3 12 1 2 ∂ U ∂ x ∂ 4 C ∂ x 4 + ∂ 2 U ∂ x 2 ∂ 3 C ∂ x 3 + ∂ 3 U ∂ x 3 ∂ 2 C ∂ x 2 + 1 2 ∂ 4 U ∂ x 4 ∂ C ∂ x + O ( k 4 ) + O ( h 4 ) .

Here, we denote y by:

y = 2 U C − h ∂ U ∂ x C + h 2 2 2 U ∂ 2 C ∂ x 2 + 2 ∂ U ∂ x ∂ C ∂ x + ∂ 2 U ∂ x 2 C − h 3 2 ∂ U ∂ x ∂ 2 C ∂ x 2 + ∂ 2 U ∂ x 2 ∂ C ∂ x + 1 3 ∂ 3 U ∂ x 3 C .

As k , h → 0 , we have

(37) ∂ U ∂ t − U ( 1 − U ) + ∂ U ∂ x ∂ C ∂ x + U ∂ 2 C ∂ x 2 = − C 2 k h 2 ∂ U ∂ t + k 2 h 2 ∂ 2 U ∂ t 2 + k 3 3 h 2 ∂ 3 U ∂ t 3 .

From equation (37), we conclude that NSFD3 scheme is not consistent with the PDE for invasive cell concentration. We note that the consistency of NSFD3 was not checked in [18].

6.2 Numerical results using NSFD3

In Figures 13 and 14, we present plots of the numerical solutions for concentration of IC, CT, and protease using NSFD3 with k = 1.0 , h = 1.8538 , and k = 0.4055 , h = 1.0 , respectively.

$Figure 13 Results for cancer growth model using NSFD3 with k = 1.0 k=1.0 and h = 1.8538 h=1.8538 , i.e. R = 0.5 R=0.5 (case 1). (a) Plot of numerical solutions vs t t at x ( 4 ) = 5.5614 x\left(4)=5.5614 , (b) plot of numerical solution for conc. of IC vs x x vs t t , (c) plot of numerical solution for conc. of CT vs x x vs t t , and (d) plot of numerical solution for conc. of protease vs x x vs t t .$

Figure 13

Results for cancer growth model using NSFD3 with k = 1.0 and h = 1.8538 , i.e. R = 0.5 (case 1). (a) Plot of numerical solutions vs t at x ( 4 ) = 5.5614 , (b) plot of numerical solution for conc. of IC vs x vs t , (c) plot of numerical solution for conc. of CT vs x vs t , and (d) plot of numerical solution for conc. of protease vs x vs t .

Figure 14

Results for cancer growth model using NSFD3 with k = 0.4055 and h = 1.0 , i.e. R = 0.5 (case 1). (a) Plot of numerical solutions vs t at x = 5 , (b) plot of numerical solution for conc. of IC vs x vs t , (c) plot of numerical solution for conc. of CT vs x vs t , and (d) plot of numerical solution for conc. of protease vs x vs t .

Figures 13 and 14 show that numerical solutions obtained using NSFD3 are non-negative.

The profiles from NSFD2 and NSFD3 are very different from each other though they are obtained on a slight modification of NSFD1.

From Figures 13 and 14, we see that the initial change in concentration gradient is quite delayed, as previously mentioned, this initial change should occur at t approximately equal to 5.

7 NSFD4 to solve the malignant invasion model

We consider NSFD2 given by equation (30) in a general setting without fixing the value of R . We rewrite equation (30) as:

(38) U m n + 1 = 0.5 ( U m n + U m − 1 n ) + U m n ( ϕ ( k ) − R C m + 1 n + R C m n ) + R U m − 1 n ( C m n − C m − 1 n ) ( 1 + ϕ ( k ) U m n ) ,

where ϕ ( k ) = exp ( k ) − 1 , ψ ( h ) = h , and R = ϕ ( k ) [ ψ ( h ) ] 2 .

We can rewrite equation (38) as:

(39) U m n + 1 = U m n ( 0.5 + ϕ ( k ) − R C m + 1 n + R C m n ) + U m − 1 n ( 0.5 + R C m n − R C m − 1 n ) ( 1 + ϕ ( k ) U m n ) .

From equation (39), we see that in order for NSFD4 to preserve positivity, we require that R ≤ 1 2 + ϕ ( k ) and R ≥ 1 2 .

To check for consistency of NSFD4 scheme, we consider equation (39). We obtain the Taylor series expansion and after some re-arrangement using k ϕ ( k ) ≈ 1 , we obtain the modified equation as:

(40) ∂ U ∂ t − U ( 1 − U ) + ∂ U ∂ x ∂ C ∂ x + U ∂ 2 C ∂ x 2 = − h 2 ϕ ( k ) ∂ U ∂ x + h 2 4 ϕ ( k ) ∂ 2 U ∂ x 2 − h 3 12 ϕ ( k ) ∂ 3 U ∂ x 3 − k 2 ∂ 2 U ∂ t 2 + 2 U ∂ U ∂ t − k 2 6 ∂ 3 U ∂ t 3 + 3 U ∂ 2 U ∂ t 2 − k 3 24 ∂ 4 U ∂ t 4 + 4 U ∂ 3 U ∂ t 3 + O ( k 4 ) + O ( h 4 ) .

After taking k , h → 0 , we have

(41) ∂ U ∂ t − U ( 1 − U ) + ∂ U ∂ x ∂ C ∂ x + U ∂ 2 C ∂ x 2 = − h 2 ϕ ( k ) ∂ U ∂ x − h 2 ∂ 2 U ∂ x 2 + h 2 6 ∂ 3 U ∂ x 3 .

From equation (41), we conclude that NSFD4 scheme is not consistent with the PDE for invasive cell concentration. Hence, NSFD4 is not a reliable method to solve the malignant invasion model. We were not successful to construct a scheme that is consistent and that preserves the positivity of solutions by working with the general form of NSFD2 scheme.

8 NSFD5 to solve the malignant invasion model

8.1 Construction and analysis of NSFD5

In this section, we try to make NSFD3 consistent. We therefore consider

U m n + 1 = ( 1 + ϕ ( k ) ) U m n + R ( U m n + U m − 1 n ) C m n 1 + ϕ ( k ) U m n + R ( U m n C m + 1 n + U m − 1 n C m − 1 n ) U m n .

We rewrite the scheme as:

(42) U m n + 1 = ( 1 + ϕ ( k ) ) ( U m n ) 2 + R ( U m n + U m − 1 n ) C m n U m n U m n + ϕ ( k ) ( U m n ) 2 + R ( U m n C m + 1 n + U m − 1 n C m − 1 n ) .

We consider

U m n + 1 [ U m n + ϕ ( k ) ( U m n ) 2 + R ( U m n C m + 1 n + U m − 1 n C m − 1 n ) ] = ( 1 + ϕ ( k ) ) ( U m n ) 2 + R ( U m n + U m − 1 n ) C m n U m n .

Taylor series expansion about ( t n , x m ) gives

(43) ∂ U ∂ t − U ( 1 − U ) + ∂ U ∂ x ∂ C ∂ x + U ∂ 2 C ∂ x 2 = − k 2 ∂ 2 U ∂ t 2 + 2 U ∂ U ∂ t − k 2 6 ∂ 3 U ∂ t 3 + 3 U ∂ 2 U ∂ t 2 − k 3 24 ∂ 4 U ∂ t 4 + 4 U ∂ 3 U ∂ t 3 − k 2 4 h 2 ∂ U ∂ t C − 2 h U ∂ U ∂ t ∂ U ∂ x C + 1 U ∂ U ∂ t 2 U ∂ 2 C ∂ x 2 + 2 ∂ U ∂ x ∂ C ∂ x + ∂ 2 U ∂ x 2 C + k 2 h U ∂ U ∂ t ∂ U ∂ x ∂ 2 C ∂ x 2 + ∂ 2 U ∂ x 2 ∂ C ∂ x + 1 3 ∂ 3 U ∂ x 3 − k 2 6 6 h 2 ∂ 2 U ∂ t 2 C − 3 h U ∂ 2 U ∂ t 2 ∂ U ∂ x C + 3 2 U ∂ 2 U ∂ t 2 2 U ∂ 2 C ∂ x 2 + 2 ∂ U ∂ x ∂ C ∂ x + ∂ 2 U ∂ x 2 C + k 2 6 3 h 2 U ∂ 2 U ∂ t 2 ∂ U ∂ x ∂ 2 C ∂ x 2 + ∂ 2 U ∂ x 2 ∂ C ∂ x + 1 3 ∂ 3 U ∂ x 3 − k 3 24 8 h 2 ∂ 3 U ∂ t 3 C − 4 h U ∂ 3 U ∂ t 3 ∂ U ∂ x C + 2 U ∂ 3 U ∂ t 3 2 U ∂ 2 C ∂ x 2 + 2 ∂ U ∂ x ∂ C ∂ x + ∂ 2 U ∂ x 2 C + k 3 24 2 h U ∂ 3 U ∂ t 3 ∂ U ∂ x ∂ 2 C ∂ x 2 + ∂ 2 U ∂ x 2 ∂ C ∂ x + 1 3 ∂ 3 U ∂ x 3 + h 2 ∂ U ∂ x ∂ 2 C ∂ x 2 + ∂ 2 U ∂ x 2 ∂ C ∂ x − h 2 2 1 3 ∂ U ∂ x ∂ 3 C ∂ x 3 + 1 2 ∂ 2 U ∂ x 2 ∂ 2 C ∂ x 2 + 1 3 ∂ 3 U ∂ x 3 ∂ C ∂ x + 1 6 U ∂ 4 C ∂ x 4 + h 3 12 1 2 ∂ U ∂ x ∂ 4 C ∂ x 4 + ∂ 2 U ∂ x 2 ∂ 3 C ∂ x 3 + ∂ 3 U ∂ x 3 ∂ 2 C ∂ x 2 + 1 2 ∂ 4 U ∂ x 4 ∂ C ∂ x + O ( k 4 ) + O ( h 4 ) .

If we choose k → 0 , h → 0 , and k ≪ h , it is possible to make NSFD3 consistent.

Choosing different combinations of k and h results in the following modified equations:

For k = 0.001 and h = 1.0 , we have

(44) ∂ U ∂ t − U ( 1 − U ) + ∂ U ∂ x ∂ C ∂ x + U ∂ 2 C ∂ x 2 = − 1 0 − 3 2 ∂ 2 U ∂ t 2 + 2 U ∂ U ∂ t − 1 0 − 6 6 ∂ 3 U ∂ t 3 + 3 U ∂ 2 U ∂ t 2 − 1 0 − 9 24 ∂ 4 U ∂ t 4 + 4 U ∂ 3 U ∂ t 3 − 1 0 − 3 2 4 ∂ U ∂ t C − 2 U ∂ U ∂ t ∂ U ∂ x C + 1 U ∂ U ∂ t 2 U ∂ 2 C ∂ x 2 + 2 ∂ U ∂ x ∂ C ∂ x + ∂ 2 U ∂ x 2 C + 1 0 − 3 2 1 U ∂ U ∂ t ∂ U ∂ x ∂ 2 C ∂ x 2 + ∂ 2 U ∂ x 2 ∂ C ∂ x + 1 3 ∂ 3 U ∂ x 3 − 1 0 − 6 6 6 ∂ 2 U ∂ t 2 C − 3 U ∂ 2 U ∂ t 2 ∂ U ∂ x C + 3 2 U ∂ 2 U ∂ t 2 2 U ∂ 2 C ∂ x 2 + 2 ∂ U ∂ x ∂ C ∂ x + ∂ 2 U ∂ x 2 C + 1 0 − 6 6 3 2 U ∂ 2 U ∂ t 2 ∂ U ∂ x ∂ 2 C ∂ x 2 + ∂ 2 U ∂ x 2 ∂ C ∂ x + 1 3 ∂ 3 U ∂ x 3 − 1 0 − 9 24 8 ∂ 3 U ∂ t 3 C − 4 U ∂ 3 U ∂ t 3 ∂ U ∂ x C + 2 U ∂ 3 U ∂ t 3 2 U ∂ 2 C ∂ x 2 + 2 ∂ U ∂ x ∂ C ∂ x + ∂ 2 U ∂ x 2 C + 1 0 − 9 24 2 U ∂ 3 U ∂ t 3 ∂ U ∂ x ∂ 2 C ∂ x 2 + ∂ 2 U ∂ x 2 ∂ C ∂ x + 1 3 ∂ 3 U ∂ x 3 + 1 2 ∂ U ∂ x ∂ 2 C ∂ x 2 + ∂ 2 U ∂ x 2 ∂ C ∂ x − 1 2 1 3 ∂ U ∂ x ∂ 3 C ∂ x 3 + 1 2 ∂ 2 U ∂ x 2 ∂ 2 C ∂ x 2 + 1 3 ∂ 3 U ∂ x 3 ∂ C ∂ x + 1 6 U ∂ 4 C ∂ x 4 + 1 12 1 2 ∂ U ∂ x ∂ 4 C ∂ x 4 + ∂ 2 U ∂ x 2 ∂ 3 C ∂ x 3 + ∂ 3 U ∂ x 3 ∂ 2 C ∂ x 2 + 1 2 ∂ 4 U ∂ x 4 ∂ C ∂ x + …

For k = 0.00025 , h = 0.5 , we have

(45) ∂ U ∂ t − U ( 1 − U ) + ∂ U ∂ x ∂ C ∂ x + U ∂ 2 C ∂ x 2 = − 2.5 × 1 0 − 4 2 ∂ 2 U ∂ t 2 + 2 U ∂ U ∂ t − 6.25 × 1 0 − 8 6 ∂ 3 U ∂ t 3 + 3 U ∂ 2 U ∂ t 2 − 1.5625 × 1 0 − 11 24 ∂ 4 U ∂ t 4 + 4 U ∂ 3 U ∂ t 3 − 2.5 × 1 0 − 4 2 16 ∂ U ∂ t C − 4 U ∂ U ∂ t ∂ U ∂ x C + 1 U ∂ U ∂ t 2 U ∂ 2 C ∂ x 2 + 2 ∂ U ∂ x ∂ C ∂ x + ∂ 2 U ∂ x 2 C + 2.5 × 1 0 − 4 2 1 2 U ∂ U ∂ t ∂ U ∂ x ∂ 2 C ∂ x 2 + ∂ 2 U ∂ x 2 ∂ C ∂ x + 1 3 ∂ 3 U ∂ x 3 − 6.25 × 1 0 − 8 6 24 ∂ 2 U ∂ t 2 C − 6 U ∂ 2 U ∂ t 2 ∂ U ∂ x C + 3 2 U ∂ 2 U ∂ t 2 2 U ∂ 2 C ∂ x 2 + 2 ∂ U ∂ x ∂ C ∂ x + ∂ 2 U ∂ x 2 C + 6.25 × 1 0 − 8 6 3 4 U ∂ 2 U ∂ t 2 ∂ U ∂ x ∂ 2 C ∂ x 2 + ∂ 2 U ∂ x 2 ∂ C ∂ x + 1 3 ∂ 3 U ∂ x 3 − 1.5625 × 1 0 − 11 24 32 ∂ 3 U ∂ t 3 C − 8 U ∂ 3 U ∂ t 3 ∂ U ∂ x C + 2 U ∂ 3 U ∂ t 3 2 U ∂ 2 C ∂ x 2 + 2 ∂ U ∂ x ∂ C ∂ x + ∂ 2 U ∂ x 2 C + 1.5625 × 1 0 − 11 24 4 U ∂ 3 U ∂ t 3 ∂ U ∂ x ∂ 2 C ∂ x 2 + ∂ 2 U ∂ x 2 ∂ C ∂ x + 1 3 ∂ 3 U ∂ x 3 + 1 4 ∂ U ∂ x ∂ 2 C ∂ x 2 + ∂ 2 U ∂ x 2 ∂ C ∂ x − 1 8 1 3 ∂ U ∂ x ∂ 3 C ∂ x 3 + 1 2 ∂ 2 U ∂ x 2 ∂ 2 C ∂ x 2 + 1 3 ∂ 3 U ∂ x 3 ∂ C ∂ x + 1 6 U ∂ 4 C ∂ x 4 + 1 96 1 2 ∂ U ∂ x ∂ 4 C ∂ x 4 + ∂ 2 U ∂ x 2 ∂ 3 C ∂ x 3 + ∂ 3 U ∂ x 3 ∂ 2 C ∂ x 2 + 1 2 ∂ 4 U ∂ x 4 ∂ C ∂ x + …

From equations (44) and (45), we can conclude that choices such as k = 0.001 , h = 1.0 and k = 0.00025 , h = 0.5 appear to be suitable in order to make NSFD3 consistent.

In an extreme situation, the maximum value of U m n + 1 is

( 1 + ϕ ( k ) ( 1 ) ) + ϕ ( k ) [ ψ ( h ) ] 2 ( 1 + 1 ) ( 1 ) ( 1 + ϕ ( k ) ( 0 ) ) + ϕ ( k ) [ ψ ( h ) ] 2 ( 0 ) ∕ 1 = 1 + ϕ ( k ) + 2 ϕ ( k ) [ ψ ( h ) ] 2 ,

as k → 0 , h → 0 , and k h 2 → 0 , we have that U m n + 1 → 1 .

Error estimates: Let ( u , c , p ) ∈ C 2 , 2 ( θ ) , where θ = { ( x , t ) ∣ 0 ≤ x ≤ 1 , 0 < t ≤ T } for some T > 0 . If h and k are such that

ϕ ( k ) [ ψ ( h ) ] 2 = exp ( k ) − 1 h 2 → 0 ,

then there exists constants C u , C > 0 , independent of h and k , such that

max ( x m , t n ) ∈ θ ∣ u m n − U m n ∣ ≤ max 0 ≤ m ≤ N + 1 ∣ u m 0 − U m 0 ∣ ( ϕ ( k ) + 1 ) ( ∣ u m 0 − U m 0 ∣ ϕ ( k ) + 1 ) T k + C u M , max ( x m , t n ) ∈ θ ∣ c m n − C m n ∣ ≤ max 0 ≤ m ≤ N + 1 ∣ c m 0 − C m 0 ∣ ( ∣ c m 0 − C m 0 ∣ ϕ ( k ) + 1 ) T k + C M c ,

and

max ( x m , t n ) ∈ θ ∣ p m n − P m n ∣ ≤ max 0 ≤ m ≤ N + 1 ∣ p m 0 − P m 0 ∣ + C M p ,

where

M = max ( x m , t n ) ∈ θ ∂ 2 U ∂ t 2 ( x , t ) + 2 u ( x , t ) ∂ U ∂ t ( x , t ) , ∂ U ∂ x ( x , t ) ∂ 2 C ∂ x 2 ( x , t ) + ∂ 2 U ∂ x 2 ( x , t ) ∂ C ∂ x ( x , t ) ,

M c = max ( x m , t n ) ∈ θ ∂ 2 C ∂ t 2 ( x , t ) + 2 p ( x , t ) ∂ C ∂ t ( x , t ) ,

and

M p = max ( x m , t n ) ∈ θ ∂ 2 P ∂ t 2 ( x , t ) + 2 ε − 1 ∂ P ∂ t ( x , t ) − u ( x , t ) ∂ C ∂ t ( x , t ) for m = 1 , 2 , 3 , … , N .

8.2 Numerical results using NSFD5 (case 1)

Figures 15 and 16 give the plot of the numerical solutions for concentration of IC, CT, and protease using NSFD5 with h = 0.5 , k = 0.00025 ( R = 0.001 ) and h = 1 , 0 , k = 0.001 ( R = 0.001 ), respectively.

$Figure 15 Plot of numerical solutions using NSFD5 with k = 0.00025 k=0.00025 and h = 0.5 ( R = 0.001 ) h=0.5\hspace{0.33em}\left(R=0.001) (case 1). (a) Plot of numerical solutions vs t t at x ( 11 ) = 5 x\left(11)=5 , (b) plot of numerical solution for conc. of IC vs x x vs t t , (c) plot of numerical solution for conc. of CT vs x x vs t t , and (d) plot of numerical solution for conc. of protease vs x x vs t t .$

Figure 15

Plot of numerical solutions using NSFD5 with k = 0.00025 and h = 0.5 ( R = 0.001 ) (case 1). (a) Plot of numerical solutions vs t at x ( 11 ) = 5 , (b) plot of numerical solution for conc. of IC vs x vs t , (c) plot of numerical solution for conc. of CT vs x vs t , and (d) plot of numerical solution for conc. of protease vs x vs t .

Figure 16

Plot of numerical solutions vs x at some values of t using NSFD5 with k = 0.00025 and h = 0.5 (case 1). (a) Plot of numerical solutions vs x at t = 5, (b) plot of numerical solution vs x at t = 15, (c) plot of numerical solutions vs x at t = 25, and (d) plot of numerical solutions vs x at t = 35.

Figures 17 and 18 depict the plot of the numerical solutions of concentration of IC, CT, and protease vs x at some values of t using NSFD5 with h = 0.5 , k = 0.00025 and h = 1 , 0 , k = 0.001 , respectively.

$Figure 17 Plot of numerical solutions using NSFD5 with k = 0.001 k=0.001 and h = 1.0 ( R = 0.001 ) h=1.0\hspace{0.33em}\left(R=0.001) (case 1). (a) Plot of numerical solutions vs t at x(6) = 5, (b) plot of numerical solution for conc. of IC vs x vs t, (c) plot of numerical solution for conc. of CT vs x vs t, and (d) plot of numerical solution for conc. of protease vs x vs t.$

Figure 17

Plot of numerical solutions using NSFD5 with k = 0.001 and h = 1.0 ( R = 0.001 ) (case 1). (a) Plot of numerical solutions vs t at x(6) = 5, (b) plot of numerical solution for conc. of IC vs x vs t, (c) plot of numerical solution for conc. of CT vs x vs t, and (d) plot of numerical solution for conc. of protease vs x vs t.

Figure 18

Plot of numerical solutions vs x at some values of t using NSFD5 with k = 0.001 and h = 1.0 (case 1). (a) Plot of numerical solutions vs x at t = 5 , (b) plot of numerical solution vs x at t = 15 , (c) plot of numerical solutions vs x at t = 25 , and (d) plot of numerical solutions vs x at t = 35 .

Tables 1 and 2 give the maximum norm errors E k and the numerical rate of convergence in time R T for NSFD5.

Table 1

E k errors and the numerical rate of convergence in time using NSFD5 with h = 1.0 at time t = 1.0 and x ∈ [ 0 , 20 ] (case 1)

k	E k (IC)	E k (CT)	E k (protease)	R T (IC)	R T (CT)	R T (protease)
1.0 × 1 0 ‒ 3	1.8803 × 1 0 ‒ 3	6.8828 × 1 0 ‒ 4	1.5879 × 1 0 ‒ 3
5.0 × 1 0 ‒ 4	9.4271 × 1 0 ‒ 4	3.4723 × 1 0 ‒ 4	7.9974 × 1 0 ‒ 4	0.9961	0.9871	0.9895
2.5 × 1 0 ‒ 4	4.7198 × 1 0 ‒ 4	1.7439 × 1 0 ‒ 4	4.0123 × 1 0 ‒ 4	0.9981	0.9935	0.9951
1.25 × 1 0 ‒ 4	2.3615 × 1 0 ‒ 4	8.7392 × 1 0 ‒ 5	2.0093 × 1 0 ‒ 4	0.9990	0.9968	0.9977

Table 2

E k errors and the numerical rate of convergence in time using NSFD5 with h = 0.5 at time t = 1.0 and x ∈ [ 0 , 20 ] (case 1)

k	E k (IC)	E k (CT)	E k (protease)	R T (IC)	R T (CT)	R T (protease)
2.5 × 1 0 ‒ 4	1.6899 × 1 0 ‒ 3	4.8855 × 1 0 ‒ 4	1.3425 × 1 0 ‒ 3
1.25 × 1 0 ‒ 4	8.4874 × 1 0 ‒ 4	2.4553 × 1 0 ‒ 4	6.7460 × 1 0 ‒ 4	0.9936	0.9926	0.9929
6.25 × 1 0 ‒ 5	4.2531 × 1 0 ‒ 4	1.2308 × 1 0 ‒ 4	3.3813 × 1 0 ‒ 4	0.9968	0.9963	0.9964
3.125 × 1 0 ‒ 5	2.1289 × 1 0 ‒ 4	6.1618 × 1 0 ‒ 5	1.6928 × 1 0 ‒ 4	0.9984	0.9981	0.9982

Figures 15 to 18 demonstrate that numerical solutions are non-negative using NSFD5.

The profiles obtained in Figures 15 and 16 for NSFD5 share some similarities to those of NSFD1 in Figures 5 and 6 and are quite similar to those of SFD1 in Figures 3 and 4.

From Figures 17 and 18, we observe that the numerical solutions of concentration of IC, CT, and protease vs x at some values of t using NSFD5 share some resemblance to those of NSFD1 in Figure 7. The numerical waves appear to travel at the same speeds.

We also note that the profiles in Figure 18 resemble those in Figure 5 from [28], where a finite elemental method is used to solve a cancer growth model.

8.3 Numerical results using NSFD5 (case 2)

Figures 19 and 20 display the plot of the numerical solutions for the concentration of IC, CT, and protease using NSFD5 for case 2 with h = 0.5 , k = 0.00025 and h = 1 , 0 , k = 0.001 , respectively.

$Figure 19 Plot of numerical solutions using NSFD5 with k = 0.00025 k=0.00025 and h = 0.5 ( R = 0.001 ) h=0.5\hspace{0.33em}\left(R=0.001) (case 2). (a) Plot of numerical solutions vs t t at x ( 11 ) = 5.5614 x\left(11)=5.5614 , (b) plot of numerical solution for conc. of IC vs x x vs t t , (c) plot of numerical solution for conc. of CT vs x x vs t t , and (d) plot of numerical solution for conc. of protease vs x x vs t t .$

Figure 19

Plot of numerical solutions using NSFD5 with k = 0.00025 and h = 0.5 ( R = 0.001 ) (case 2). (a) Plot of numerical solutions vs t at x ( 11 ) = 5.5614 , (b) plot of numerical solution for conc. of IC vs x vs t , (c) plot of numerical solution for conc. of CT vs x vs t , and (d) plot of numerical solution for conc. of protease vs x vs t .

Figure 20

Plot of numerical solutions vs x at some values of t using NSFD5 with k = 0.00025 and h = 0.5 (case 2). (a) Plot of numerical solutions vs x at t = 0, (b) plot of numerical solution vs x at t = 4, (c) plot of numerical solutions vs x at t = 8, and (d) plot of numerical solutions vs x at t = 24.

Figures 21 and 22 give the plot of the numerical solutions of concentration of IC, CT, and protease vs x at some values of t using NSFD5 and case 2 with h = 0.5 , k = 0.00025 and h = 1 , 0 , k = 0.001 , respectively.

Figure 21

Plot of numerical solutions using NSFD5 and case 2 with k = 0.001 and h = 1.0 (case 2). (a) Plot of numerical solutions vs t at x = 5, (b) plot of numerical solution for conc. of IC vs x vs t, (c) plot of numerical solution for conc. of CT vs x vs t, and (d) plot of numerical solution for conc. of protease vs x vs t.

Figure 22

Plot of numerical solutions vs x at some values of t using NSFD5 with k = 0.001 and h = 1.0 (case 2). (a) Plot of numerical solutions vs x at t = 0 , (b) plot of numerical solution vs x at t = 4 , (c) plot of numerical solutions vs x at t = 8 , and (d) plot of numerical solutions vs x at t = 24 .

Figures 19–22 demonstrate that numerical solutions are non-negative using NSFD5 and case 2.

Figures 21 and 22 illustrate that the numerical solution for concentration of invasive cells converges to 1 as we progress in time and space.

9 Conclusion

In this work, one SFD and a few NSFD methods were used to solve a continuous malignant invasion model. Stability of SFD methods is not easy to obtain theoretically, and the solutions do not always preserve positivity; this is one of the challenges of classical schemes in this work. Three NSFD methods from [18] were considered to solve the malignant invasion model, namely, NSFD1, NSFD2, and NSFD3. We show that NSFD1 is consistent with the continuous model, and it is positivity-preserving if C m n ≥ C m − 1 n ; therefore, the choice of profile determines whether the solutions remain non-negative as we progress in time (Figures 6 and 8). We found that NSFD2 is not consistent with the PDE for invasive cell concentration, though it is positivity-preserving when ϕ ( k ) [ ψ ( h ) ] 2 = 1 2 ; NSFD2 is not useful in this work. Similarly, we found that NSFD3 is not consistent with the PDE for invasive cell concentration, though it is unconditionally positivity-preserving. The profiles obtained from NSFD2 and NSFD3 were very different from each other (Figures 11–14), likely due to them being inconsistent methods. NSFD3 being unconditionally dynamically consistent with respect to positivity proved to be useful, and using this property, we extend NSFD3 by appropriately choosing the step sizes ( k , h ) to obtain NSFD5. We show that NSFD5 preserves positivity and is consistent if, for instance, k , h is chosen such that k = 0.001 , h = 1.0 and k = 0.00025 , h = 0.5 . We also show that using these k , h values, the numerical solutions more accurately represent the expected solutions of the continuous malignant invasion model.

Acknowledgements

The authors would like to thank the anonymous reviewers for the constructive comments and suggestions that have enabled the authors to significantly improve this manuscript. G.N. de Waal is grateful to the Department of Mathematics and Applied Mathematics for some funding in 2021 and 2022. A.R. Appadu is grateful to Nelson Mandela University and the National Research Foundation of South Africa (Grant Number: 136161), which allowed this work to be carried out.

Funding information: Appadu is grateful to Nelson Mandela University to allow him to access research funds for the open access publication fees.
Author contributions: The topic of the research was provided by Appadu. Derivation of methods and writing up of this manuscript were done by both authors. Coding was mainly done by de Waal under the supervision of Appadu.
Conflict of interest: The authors state no conflict of interest.

References

[1] A. M. Wazwaz, Partial Differential Equations and Solitary Wave Theory, Springer, Berlin, 200910.1007/978-3-642-00251-9Search in Google Scholar

[2] A. Jefferey, Applied Partial Differential Equations: An Introduction, Academic Press, USA, 2002, 978-0123822529, 389. Search in Google Scholar

[3] A. D. Polyanin, Partial Differential Equations, 2008, [online]. Available from: http://www.scholarpedia.org/article/Partial_differential_equation [ accessed 13 September 2021]. 10.4249/scholarpedia.4605Search in Google Scholar

[4] J. C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, 2nd edition, SIAM, Wisconsin, 2004, 978-0898715675, 431. 10.1137/1.9780898717938Search in Google Scholar

[5] M. M. Khalsaraei, Sh. Heydari, and L. D. Algoo, Positivity preserving nonstandard finite difference schemes applied to cancer growth model, J. Cancer Treat. Res 4 (2017), no. 4, 27–33, DOI: https://10.11648/j.jctr.20160404.11. Search in Google Scholar

[6] R. E. Mickens, Advances in the Application of Nonstandard Finite Difference Schemes, World Scientific, Singapore, 2005, 978-9812564047, 664. 10.1142/5884Search in Google Scholar

[7] G. F. Sun, G. R. Liu, and M. Li, An efficient explicit finite-difference scheme for simulating coupled biomass growth on nutritive substrates, Math. Probl. Eng. 2015 (2015), 1–17, DOI: https://doi.org/10.1155/2015/708497. 10.1155/2015/708497Search in Google Scholar

[8] G. F. Sun, G. R. Liu, and M. Li, A novel explicit positivity-preserving finite-difference scheme for simulating bounded growth of biological films, Int. J. Comput. Methods 13 (2016), no. 2, 1640014, DOI: https://doi.org/10.1142/S0219876216400132. 10.1142/S0219876216400132Search in Google Scholar

[9] Y. Yu, W. Deng, and Y. Wu, Positivity and boundedness preserving schemes from space-time fractional predator-prey reaction-diffusion model, Comp. Math. Appl. 69 (2015), no. 8, 743–759, DOI: https://doi.org/10.1016/j.camwa.2015.02.024. 10.1016/j.camwa.2015.02.024Search in Google Scholar

[10] R. E. Mickens, Exact solutions to a finite-difference model of a nonlinear reaction-advection equation: Implications for numerical analysis, Numer. Methods Partial Differ. Equ. 5 (1989), no. 4, 313–325, DOI: https://doi.org/10.1002/num.1690050404. 10.1002/num.1690050404Search in Google Scholar

[11] R. Anguelov and J. M-S. Lubuma, Contributions to the mathematics of the nonstandard finite difference method and applications, Numer. Methods Partial Differ. Equ. 17 (2001), no. 5, 518–543, DOI: https://doi.org/10.1002/num.1025. 10.1002/num.1025Search in Google Scholar

[12] F. B. Hildebrand, Finite-Difference Equations and Simulations, Prentice-Hall, New Jersey, 1968, 978-0133172300, 338. Search in Google Scholar

[13] R. E. Mickens, Nonstandard Finite Difference Models of Differential Equations, World Scientific, Singapore, 1994, 978-9810214586, 264. 10.1142/2081Search in Google Scholar

[14] R. Anguelov and J. M-S. Lubuma, Nonstandard finite difference method by nonlocal approximation, Math. Comput. Simul. 61 (2003), no. 2003, 465–475, DOI: https://doi.org/10.1016/S0378-4754(02)00106-4. 10.1016/S0378-4754(02)00106-4Search in Google Scholar

[15] R. Anguelov, P. Kama, and J. M-S. Lubuma, On non-standard finite difference models of reaction-diffusion equations, J. Comput. Appl. Math. 175 (2005), no. 1, 11–29, DOI: https://doi.org/10.1016/j.cam.2004.06.002. 10.1016/j.cam.2004.06.002Search in Google Scholar

[16] H. P. Bhatt and A. Q. M. Khaliq, Fourth-order compact schemes for the numerical simulation of coupled Burgers’ equation. Comput. Phys. Commun. 200 (2015), 117–138. 10.1016/j.cpc.2015.11.007Search in Google Scholar

[17] Y. O. Tijani, A. R. Appadu, and A. A. Aderogba, Some finite difference methods to model biofilm growth and decay: classical and non-standard, Comput. 9 (2021), no. 11, 123, DOI: https://doi.org/10.3390/computation9110123. 10.3390/computation9110123Search in Google Scholar

[18] M. Chapwanya, J. M-S. Lubuma, and R. E. Mickens, Positivity-preserving nonstandard finite difference schemes for cross-diffusion equations in biosciences, Comput. Math. Appl. 68 (2014), no. 9, 1071–1082, DOI: https://doi.org/10.1016/j.camwa.2014.04.021. 10.1016/j.camwa.2014.04.021Search in Google Scholar

[19] N. Agoulmine, Autonomic Network Management and Principles from Concepts to Applications, Elsevier, California, 2009, 978-0123821904, 281. Search in Google Scholar

[20] V. Volpert and S. Petrovskii, Reaction-diffusion waves in biology, Phys. Life Rev. 6 (2009), no. 4, 267–310, DOI: https://doi.org/10.1016/j.plrev.2009.10.002. 10.1016/j.plrev.2009.10.002Search in Google Scholar PubMed

[21] B. P. Marchant, J. Norbury, and A. J. Perumpanani, Traveling shock waves arising in a model of malignant invasion, SIAM J. Appl. Math. 60 (2000), no. 2, 463–476, DOI: https://doi.org/10.1137/S0036139998328034. 10.1137/S0036139998328034Search in Google Scholar

[22] J. D. Murray, Mathematical Biology II: Spatial Models and Biomedical Applications, 3rd ed., Springer, New York, 2003, 978-038795228, 814. Search in Google Scholar

[23] J. D. Murray, Mathematical Biology I: An Introduction, 3rd ed., Springer, New York, 2002, 978-0387952239, 551. Search in Google Scholar

[24] L. Chen and A. Jü ngel, Analysis of a parabolic cross-diffusion population model without self-diffusion, J. Differential Equations 60 (2006), no. 1, 39–59, DOI: https://doi.org/10.1016/j.jde.2005.08.002. 10.1016/j.jde.2005.08.002Search in Google Scholar

[25] D. Le, Cross-diffusion equations on n spatial dimensional domains, Electronic Journal of Differential Equations, Fifth Mississippi State Conference on Differential Equations and Computational Simulations, vol. 10, Mississippi, 2003, pp. 193–210. Search in Google Scholar

[26] J. Shi, Z. Xie, and K. Little, Cross-diffusion induced instability and stability in reaction-diffusion systems, J. Appl. Anal. Comput. 24 (2010), no. 3, 95–119. https://jxshix.people.wm.edu/shi/2011-Shi-Xie-Little.pdf. 10.11948/2011007Search in Google Scholar

[27] M. Kovel and B. Zubik-Kowal, Numerical solutions for a model of tissue invasion and migration of tumour cells, Comput. Math. Methods Med. 2011 (2010), 1–16, DOI: https://doi.org/10.1155/2011/452320. 10.1155/2011/452320Search in Google Scholar PubMed PubMed Central

[28] M. Fuest, S. Heydari, P. Knobloch, J. Lankeit, and T. Wick, Global Existence of Classical Solutions and Numerical Simulations of a Cancer Invasion Model, 2022, https://arxiv.org/pdf/2205.08168v1.pdf. Search in Google Scholar

Received: 2022-03-30

Revised: 2023-02-26

Accepted: 2023-05-16

Published Online: 2023-07-28

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/dema-2022-0244

Keywords for this article

standard finite difference method; nonstandard finite difference method; consistency; positivity preserving; cross-diffusion

Creative Commons

BY 4.0