Investigation of Fractional Order Tumor Cell Concentration Equation Using Finite Difference Method

The fractional differential equation provides a powerful mathematical framework for modeling and understanding a wide range of complex, non-local and memory-dependent phenomena in various scientific, engineering and real-world applications. Liver metastasis is a secondary cancerous tumor that develops in the liver due to the spread of cancer cells from a primary cancer originating in another part of the human body. The primary focus of this study is to understand tumor growth in the human liver, both in the presence and absence of medication therapy. To achieve this, a temporal fractional-order parabolic partial differential equation is utilized, and its analysis is carried out using numerical methods. The Caputo derivative is employed to explore the impact of medication therapy on tumor growth. To numerically solve the mathematical model, the Crank-Nicolson finite difference method has been developed. This method is selected for its remarkable attributes, including unconditional stability and second-order accuracy in both spatial and temporal dimensions. The outcomes and insights derived from this study are effectively communicated through various graphical representations. These visual aids serve as invaluable tools for comprehending the profound impact of medication therapy on the growth of liver tumors. Through the medium of these graphical depictions, one can glean a clearer and more intuitive understanding of the complex dynamics at play. The numeric solution to this intricate problem is achieved through the implementation of algorithms meticulously crafted using versatile and powerful Python programming language. Python’s flexibility, extensive libraries, and robust numerical computing capabilities make it a good choice for handling the complexities of this study.


Introduction
Liver metastases are cancerous tumors that spread from other parts of the body to the liver.The liver performs several activities, including filtering toxic substances from the blood, generating bile, and Published Online First: February, 2024 https://doi.org/10.21123/bsj.2024.9246P-ISSN: 2078-8665 -E-ISSN: 2411-7986 Baghdad Science Journal amount of blood flow from the rest of the body.The cancer cells can travel through the bloodstream and settle in the liver 2 .Breast cancer, lung cancer, colorectal cancer and pancreatic cancer are common names of cancers that spread to liver.The diagnosis of liver metastases is usually made through magnetic resonance imaging (MRI) scans or computed tomography (CT) scans 1,3 .These tests can show the presence of tumors in the liver and help doctors determine the extent of the disease.The type of cancer, the size and number of tumors in the liver, and the overall health of the patient are important factors in the treatment of liver metastases.One of the primary challenges in modeling liver metastases is capturing the temporal dynamics of tumor growth.
Many investigators have developed mathematical models to analyze the growth of tumor cells.In 2003, Swanson et al. 3 studied virtual and real brain tumor for growth and invasion of gliomas.In 2014, Nenad Filipovic et al. 4 have used classical reaction diffusion equation to model liver metastasis disease and compare it to experimental data on liver metastasis disease that were collected by using medication therapy for a specific patient.The literature survey suggests that, numerous models have been developed to study the growth of tumor cells in liver metastases diseases.It is observed that, most of the models are developed by using ordinary and partial derivatives, which means that they are free from the memory effect.In the last two decades, fractional differential equations (FDEs) have been widely used in the modeling of real life problems [5][6][7][8][9] .
In 2014, both Larisse Bolton et al. 10 and Olaniyi Samuel Iyiola and F. D. Zaman 11 conducted significant studies related to fractional order tumor modelling.Bolton et al. 10 utilized a fractional Gompertz model to fit experimental data on tumor growth and they discovered that a fractional order model yeilds superior results compared to an ordinary model.On the other hand, Olaniyi Samuel Iyiola and F. D. Zaman 11 focused on a fractional order tumor cell concentration model and obtained the series solution using the q-Homotopy Analysis Method.In 2020, Rihan and Velmurugan 12 employed a delay differential equation of fractional order to investigate the tumor immune system with drug therapy.They successfully obtained a numerical solution of the model by employing Adams-Bashforth-Moulton predictor-corrector scheme.
Recently, Abaid Ur Rehman M et al. 13 studied a fractional order model of tumor cell in the Caputo sense, describing the rate of killing tumor cells.They obtained the computational solution using the reduced defferential transform method.It is observed that fractional-order modelling in biological systems yields more favorable results compare to integerorder modelling.The latter ignores the memory effects or aftereffects manifesting in these systems 13 .This advantage arises due to the power-law kernal of the Caputo derivative, where the memory effect becomes more pronounced at small time values.Therefore, to foster motivation, the present study utilizes the time fractional-order reaction-diffusion equation in the Caputo sence to analyze growth of tumors in the liver.Unfortunately, finding exact solutions to most FDEs is a challenging task due to their complex nature and non-linearity.To tackle this issue, numerical methods play a crucial role.Among the various numerical methods available, the finite difference method stands out an effective and commonly used approach 14 .The Crank-Nicolson method is one of the finite difference methods renowned for its stability and accuracy, especially when dealing with time-dependent problems 15,16 .
Our goal in this study is to devise and develop the Crank-Nicolson method to better comprehend the growth of tumor cells.

Basic Definitions
This section provides some basic definitions about fractional derivative operators.

Materials and Methods
The time fractional tumor growth equation (TFTGE) is derived from the classical tumor cell growth model proposed by Mubarak et al. 1 It can be represented as follows: The initial condition for the equation is given by: Additionally, the boundary conditions are as follows: In the equation, the term is defined in the Caputo sense.Here,  = (, ) represents the concentration of tumor cells at a specific point in space  and time , D denotes the coefficient of diffusion,  represents the rate of cell proliferation,   is the concentration of normal cells at time  = 0, ℎ() is the starting shape of tumor and (, ) represent a function showing an effect of applied meditation at time .The term  −   represents the relationship between normal cells and tumor cells.This term describes the flow of tumor cells in the liver.
where 𝛼 and In matrix form, the representation of the Crank-Nicolson method given in Eqs.13-17 is as follows: 1 =  0 +  0 , for  = 0, 18 Here,  is an identity matrix of order ,  0 = ( 1 0 ,  2 0 , ⋯ ,   0 )  ,   = ( otherwise. And Using the numpy, scipy and matplotlib modules of the Python programming language, the matrix system in Eqs.18-20 is easily solved to obtain the numerical solution of TFTGE in Eqs.3-5.

Stability
Here, the stability of the solution obtained by the Crank-Nicolson method Eqs.13-17 developed for the TFTGE Eqs.3-5 is proven.

Lemma 3:
The Matrix A is characterized by the following properties: Proof: Let  be an eigenvalue of a square matrix [  ].Then, using Gerschgorin's disc theorem 18 , each eigenvalue  of a square matrix [  ] satisfies at least one of the following inequalites: Now, inequality Eq. 22 is used to establish the condition (i) for the matrix A as follows: Therefore, based on condition (i), it can be deduced that ∥  ∥ 2 ≥ 1.Consequently, ∥  −1 ∥ 2 ≤ 1.

Lemma 4:
The Crank-Nicolson method, given in Eqs. 13 Consequently, for each row, the inequality 1 + 2 > 2 holds, thus confirming that matrix A is strictly diagonally dominant.As a result, matrix A is invertible, thereby ensuring the solvability of the finite difference scheme.Proof: To demonstrate the unconditional stability of the developed finite difference scheme, it is necessary to show that Where  is a positive number that does not depend on  and .
For =1, Eq. 18 can be written as This proves the result for  = 1.
Let us assume that, ∥   ∥ 2 ≤  ∥  0 ∥ 2 , for all  ≤ .Now, for  =  + 1, Eq. 20 can be written as Therefore, by induction, for all positive integer , is the vector of exact solution, the Eqs.18-20 can be written as Thereom 2: The Crank-Nicolson method described by Eqs.13-17 developed for the time fractional tumor growth equation Eqs.3-5 is unconditionally convergent.
Proof: Let,   = ( 1  ,  2  , … ,    )  be the vector of error appear in numerical solution at time level   .Furthermore, assume that Then, from Eq. 13, Now, from Eq. 14, Further, from Eq. 15, This is true for every .Therefore, it follows that Hence, by induction, for every positive integer ,

Results and Discussion
In this section, the following time fractional tumor growth equation is studied:  This phenomenon can be attributed to the application of the medication, which is intended to impede tumor growth.At the initial stage of medication, the concentration of the drug within the tumor is higher, leading to a more significant inhibitory effect on tumor growth.Consequently, the tumor size exhibits a notable reduction during this period.As time progresses and the medication continues to take effect, the impact of the medication gradually decreases.This decrease in the medication's efficacy is often attributed to various factors, such as tumor resistance or the drug eliminations from the body.As a consequence, the rate of tumor size reduction diminishes, resulting in a slower decrease in tumor size as compared to the initial stages of medication.
Fig. 2(a-d) displays the prediction of tumor size response over time in the absence of medication for distinct tumor thicknesses.It is observed that when no medication is used to stop the tumor growth, the tumor exhibits uncontrolled and unrestricted growth.The absence of any inhibitory effect from medication allows the tumor cells to proliferate rapidly, leading to a significant increase in tumor size over time.This unregulated growth can have severe implications for the patient's health and may lead to further complications and metastases.It is observed that despite the differences in tumor thickness among various cases, tumors with different thicknesses still follow the same growth pattern.This finding suggests that the tumor's intrinsic biological properties play a more dominant role in determining its growth behavior rather than its initial thickness.The results presented in Fig. 2 underscore the importance of effective medication in controlling tumor growth.
In the model developed by Filipovic et al., the behaviour of normal cells against the tumor cells is not consider.In 2017, Khanday and Nazir 2 conducted a study on the bioheat model to investigate the thermal distribution in cancerous tissue, providing valuable insights into the understanding of heat transfer in tumors.In 2020, Saqib Mubarak et al. 1 modified the model of Filipovic et al. by including the relationship between tumor cells and normal cells.They also considered the medication term to be both time and space dependent.

Theorem 1 :
The Crank-Nicolson method Eqs.13-17 used to solve the time fractional tumor growth equation Eqs.3-5 yields an unconditionally stable solution.

𝜕 2 𝑐𝜕𝑢 2 +
( −   ) − (, ),  ∈ [0, 6],  ∈ [0, 40], 0 <  ≤ 1, 27 with the initial condition: (, 0) =  0 ,  ∈ [0, 6], 28 and the boundary conditions: (0, ) =  0 , (6,)  = 0,  ∈ [0, 40].29 The literature-based physiological indicators given by Mubarak et al. 1 are considered to assess the tumor's size and its decreasing rate after the application of medication.The starting size of tumor is considered to be constant, that is,  0 = 205.67449 3 , the effect of applied medication at time  = 0 is (, ) = 0 and (, ) =   , if  > 0, where  = 1.075 and  = −0.05are simulation parameters.The concentration of normal cells is,   = 191.91449 3 .The proliferation rate is,  = 0.0281/day.The value of diffusion coefficient is,  = 0.005 cm 2 / .The system in Eqs.18-20 is solved for the time fractional tumor growth equation Eqs.27-29 using the Python program, and the results are presented graphically.The numerical solution of time fractional tumor growth equations Eqs.27-29 is obtained by setting Δ = 0.1, Δ = 0.1, and considering various values of .The model, based on the parabolic partial differential equation for predicting tumor growth in liver metastases disease with a time fractional derivative in the Caputo sense, is analyzed using the proposed Crank-Nicolson method.Graphical visualizations for different tumor thicknesses () and parameter values () are presented at various time intervals to estimate the effect of the applied medication on the tumor size.The solution of the time fractional tumor growth equation Eqs.27-29 is obtained by solving the system of equations Eqs.18-20 which provides information about the tumor size at any given time .In Fig. 1 and Fig. 2, the dynamics of tumor growth are plotted in two different situations: one in presence of medication term (i.e., (, ) ≠ 0), and the other in the absence of medication term (i.e., (, ) = 0) respectively, both for distinct thicknesses () of the tumor.Fig. 1(a-d) depicts the profiles of decreasing tumor size after medication for various tumor thickness values and different values of .From the observations made based on Fig. 1, it is evident that, during the medication period, the size of all tumors decreases rapidly.

Figure 2 .
Figure 2. Tumor size in absence of medication with respect to time for different tumor thickness  = ., . , . , . and for  = ., ., ..By comparing the dynamics in the presence and absence of medication, it becomes evident that the applied treatment plays a crucial role in restraining tumor growth and can be a critical factor in managing , and both  = [  ] and  = [  ] are square matrix of order , where