Convergence and Numerical Solution of a Model for Tumor Growth

: In this paper, we show the application of the meshless numerical method called “General-ized Finite Diference Method” (GFDM) for solving a model for tumor growth with nutrient density, extracellular matrix and matrix degrading enzymes, [recently proposed by Li and Hu]. We derive the discretization of the parabolic–hyperbolic–parabolic–elliptic system by means of the explicit formulae of the GFDM. We provide a theoretical proof of the convergence of the spatial–temporal scheme to the continuous solution and we show several examples over regular and irregular distribution of points. This shows the feasibility of the method for solving this nonlinear model appearing in Biology and Medicine in complicated and realistic domains.


Introduction
Cancer includes a large group of diseases that can start in almost any tissue or organ of the body when a cluster of cells starts to grow out of control, going beyond their usual boundaries to invade adjoining areas of the body and spread to other parts of the body. The latter process is known as metastasis and is the main cause of death from cancer. In this process, the ECM (extracellular matrix) composition plays an essential role [1,2]. The ECM is mostly made up of water, fibrous proteins such as collagen, elastin, laminin and proteoglycans. Collagen structure is very important to maintain the ECM function and its degradation by proteases can facilitate the invasion of cancer cells through the basal membrane [3]. Matrix-degrading enzymes (MDE) that break down the ECM can lead to invasion, where either individual or collections of cancerous cells can escape or separate from the tumor and migrate through the surrounding tissue. Thus these cells can also enter through the blood or lymphatic vessels and travel to distant locations which can develop into metastasis. Tumor cells also modify the behavior of the neighboring cells, so that they can serve as nutrient donors for cancer cell proliferation and spreading. Hence, an alteration in the nutrient density within the tumor microenvironment has been associated with malignancy progression and death [4].
Understanding how the mechanical microenvironment regulates cancer cell development that ends in metastasis represents a newly developing study model. The research related to this new model could lead to discovering better anti-metastasis therapeutics.
Over the last few years, mathematical modeling has played an essential role in Cancer research. Particularly, an increasing number of mathematical models describing tumour growth have been developed.
To describe the process and the related mechanism, we study a mathematical model for the influence of the extracellular matrix (ECM) on tumour evolution in terms of a system of partial differential equations, i.e., parabolic equations and a hyperbolic equation for nutrient density, MDE and ECM concentration. Moreover, it contains a non-constant coefficient µ(E) and allows for the movement of ECM fibres. The model is more reasonable and presents a realistic scenario, nonetheless it leads to a more difficult problem to analyse.
In this paper, we propose a model that describes the evolution with tumor microenvironment involving nutrient density, extracellular matrix and matrix degrading enzymes, which satisfy a coupled system of PDEs as follows: and initial conditions: σ 0 ( x), E 0 ( x), m 0 ( x), p 0 ( x). A range of values for the parameters of the problem can be found in [5,6] as well as in [1]. The first diffusion equation in (1) describes the evolution of the nutrient σ within the tumor, c is the ratio of the nutrient diffusion time scale to the tumor growth (e.g., tumor doubling) time scale and νσ describes the rate of consumption by the tumor. For the obtaining of the second equation in (1), it is taken into account that the concentration of the ECM in the system is governed by contributions from three factors: haptotaxis, degrading and production.
E is the extracellular matrix, V represents the velocity of proliferating cells, the term div(E · V) is the movement of ECM owing to the cell proliferation V, −γmE represents the degrading of ECM by MDE, where m is the concentration of MDE and φ(E) is a positive term representing reorganization of ECM. Since the growth rate of ECM is smaller when the ECM is denser, φ(E) is a positive monotone decreasing function of E.
It is well-known that the matrix degrading enzymes (MDE) is produced by the tumor to degrade ECM so that the cells can escape. In the MDE's equation in (1) ∆m represents diffusion, D is the constant diffusion coefficient, α is a constant production rate by the tumor and −βm represents natural decay. The fourth and fifth equations in (1) are deduced by employing the Darcy's law and the conservation of mass: p the pressure within the tumor resulting from this proliferation, the coefficient f (E) depends on the density of the porous medium, representing a mobility that reflects the combined effects of cell-cell and cell-matrix adhesion, i.e., f (E) depends on the amount of ECM present in the tumor. In the last equation S = g(E)(σ − σ) is a linear approximation for the proliferation S, where g(E)σ represents the growth rate and g(E)σ is the death rate from apoptosis, g(E) is a function of the ECM concentration. Hence, it yields div V = S.
Operating with the fourth and fifth equations of (1), the system can be rewritten: (3) In the system under consideration, c, µ, γ, D α, β, σ are constants, as we have previously detailed. The functions f and g are decreasing, strictly positive and bounded. The function φ is positive and decreasing.
The importance of models such as the above, which present a high degree of plausability, lies in their application in sciences such as biology and medicine. Their high nonlinearity, as well as the use of complicated and non-uniform domains, makes necessary the implementation of meshless methods. The Generalized Finite Difference Method (GFDM) is a meshless numerical method based on the Taylor expansion and moving least squares. This method allows us to find the discrete solution of the above system over non-uniform grids. The GFDM provides a simple discretization of the spatial derivatives at some point in terms of the values of the solution at the surrounding ones. These aspects (as well as the small number of points per star) have meant that the method has recently been used widely in [7,8] and is currently expanding [9].
The paper has the following structure: we show the fundamentals of the GFDM for completeness in Section 2. Section 3 is devoted to the proof of the convergence of the discrete solution to the continuous one of the system. We perform in the fourth section three examples with different geometries. We present in the last section some conclusions.

Fundamentals of the GFDM
Let Ω ⊂ R 2 be a bounded domain and For each one of the nodes of the domain Ω, E s -star is defined as a set of selected points E s = {x 0 ; x 1 , . . . , x s } ⊂ M with the central node x 0 ∈ M and x i (i = 1, . . . , s) ∈ M is a set of points located in the neighbourhood of x 0 . In order to select the points, different criteria as four quadrants or distance [10], can be used.
Consider an E s -star with the central node x 0 , where (x 0 , y 0 ) are the coordinates of the central node, (x i , y i ) are the coordinates of the i th node in the E s -star, and is the value of the function at the central node of the star and U i = U(x i ) are the function values at the rest of the nodes, with i = 1, ..., s, then, according to the Taylor series expansion: with i = 1, ..., s. We define: If in (4) the higher than second order terms are ignored, a second order approximation for the function U i is obtained. We denote it by u i . It is then possible to define the function B(u) as where w i = w(h i , k i ) are positive symmetrical weighting functions decreasing in magnitude as the distance to the center increases [10], of decreasing monotonically in magnitude as the distance to the center of the weighting function increases (see also Levin [11] for more details of this election).
Some weighting functions as potential 1 dist n or exponential exp(−n(dist 2 )) can be used [12], where n ∈ N. We minimize the norm given by (7) with respect to the partial derivatives by considering the following linear system where and Positive definiteness of A and the consistency of the approximation is proved in [13]. Thus, spatial derivatives using GFD [12,13] are denoted by with Finally, the time derivative is approximated as follows

GFD Scheme for System (3)
Letσ,m,Ê be aproximated solution of the system (3), that can be discretized as: Using the relation we getÊ

Convergence
Some basic and previous results are needed in order to obtain the proof of the conditional convergence of the methods. For completeness, we reproduce them: Lemma 1 ([14]). Let be a matrix N ∈ M n×n (R). If there exists some matrix norm verifying N < 1, then lim k→∞ N k = 0.
Our numerical convergence result is the following.

Theorem 1.
Let σ, m, E be the exact solution of (3), and µ and p are differentiable functions. If then the GFD explicit scheme (15) is convergent if Proof. We call ρ n j =σ n j − σ n j , η n j =m n j − m n j , ξ n j =Ê n j − E n j . Let us take the difference between the expressions of the numerical solution given by the GFD scheme (15) and the expression for the exact solution of (3). Since the function µ(E) is differentiable, we apply the Mean Value Theorem and by operating, we obtain the following: We now take bounds and callρ n = max i=0,...,s {|ρ n i |} (the same applies forξ n andη n ). Then, we haveρ We rewrite (18) in matrix form, in the following sense In order to guarantee the convergence of the scheme we use Lemmas 1 and 2 and we impose that its norm, for some matrix norm, is strictly less than 1. Take for instance the maximum of the sums by rows (infinity norm), then the condition is equivalent to Thus the inequality (21) is equivalent to Moreover, if we impose C 1 < 1, it yields the inequality (23) is equivalent to For the last term in (20), we have the inequality (25) is equivalent to

Numerical Examples
In this section we solve system (5) for the following parameters and the boundary conditions: We set as initial conditions We test the procedure above in the three clouds of points of Figure 1 (square [0, 1] × [0, 1] with regular distribution of nodes, circle with centre at (0, 0) and radius unity with both regular and irregular distributions) with 575 nodes. The criterion used to choose the nodes in the star is the distance, with eight nodes in each star, and the weighting function We solve numerically, using GFDM, system (5) with boundary conditions (27) and initial conditions (28) in the three discretised domains (clouds of points) represented in Figure 1. For each cloud of points and for time values of 2, 5 and 10 s, the graphs of the functions σ, E, m and p have been obtained. In all the graphs, the discretised domain has been represented in the x − y plane and in the z-axis the values obtained for the previous mentioned functions, i.e., σ, E, m and p.

Example 1: Square
We plot in Figures 2-4 the results obtained from the application of the GFDM to system (5) with the previous parameters for times t = 2, 5 and 10 s.

Example 2: Circle (Regular)
Figures 5-7 shows the numerical solutions given by the GFDM for times t = 2, 5 and 10 s.

Example 3: Circle (Irregular)
Finally, we test the method using the third irregular cloud of Figure 1. We present in Figures 8-10, for times t = 2, 5 and 10 s respectively, the plot of the discrete solution of system (5).

Conclusions
We have propounded a meshless numerical scheme for the recently published model (1) in [8]. Convergence of the explicit formulae is obtained under some restrictions on the time step, ∆t, depending on the distribution of the nodes and the parameters of the problem. Finally, we have applied the method for solving three cases over regular and irregular domains, showing that the GFDM can be a powerful tool for implementing realistic examples of this important model of tumor growth.