This article presents mathematical models for investigating the impact of chemotherapy on the state of a cancerous tumour based on various fractional-order derivatives. Finite difference approximations for the fractal operators of Caputo, Atangana-Baleanu, and Caputo-Fabrizio are provided and derived. The Atangana-Toufik method was employed to develop numerical algorithms for models with such fractional derivatives. Using the developed software, the dynamics of stem cell and effector cell concentrations over time in a fractal environment were studied and described, depending on the fractional order of the differential operators. This approach allows for consideration of memory effects and self-organization within the investigated fractal environment, as well as the specific features of the fractional-order operators used.
1. Introduction
Cancerous tumours exhibit complex, nonlinear behaviour, including changes in growth rate, spread, and metastasis. They operate in complex environments with varying tissue densities and blood supply. In the case of cancerous tumours, cellular processes often depend not only on the current state but also on all previous states. Traditional mathematical models may not always effectively simulate these processes. Fractional calculus allows for such features to be taken into account through fractional derivatives, which reflect the dependence of the process on the previous states of the system and enable the modelling of systems with memory effects. Furthermore, fractal analysis methods allow the construction of models that more accurately approximate real observations, such as non-stationary cell growth, the dynamics of interactions between cells, and the tumour micro-environment.
In cancer therapy, fractional calculus can be used to model various treatment approaches (e.g., chemotherapy or radiotherapy), considering how the body responds to treatment over time. This can help optimize therapeutic regimens and predict their effectiveness. Thus, fractal methods enable more precise modelling of the complex and nonlinear processes that occur in the body during tumour development. Their unique mathematical properties allow for a more accurate representation of complex biological processes and cancer cell growth dynamics. This enhances our understanding of how tumours grow, spread, and respond to therapy, which, in turn, can contribute to the development of new treatment methods and predict their efficacy.
Models based on fractional derivatives can account for delays in tumour responses to environmental changes, such as the influence of the immune system or oxygen deficiency. Various researchers have investigated mathematical models of tumour growth and treatment, and the results of these studies are presented in works [1, 2, 3, 4], among others. In particular, in [4], additional review articles are analyzed. These models provide an important role in enhancing the understanding of cancer progression. Additionally, they can be used to improve effective approaches to cancer study, treatment and prevention.
Mathematical models based on fractional differential equations can play an important role in investigating nonlocal processes in tumours and their micro-environments. In studies [5, 6, 7, 8, 9], fractional calculus was used to model diffusion and other tumour-related processes. For instance, in [5], radiobiological factors were explored to understand the complex mechanisms associated with cancer using the apparatus of fractional derivatives. In [6], a mathematical model of breast cancer was introduced using the Caputo-Fabrizio fractional derivative. The findings enable the analysis of how different treatment strategies impact the progression of breast cancer. In [7], the authors proposed a model to examine the interaction between the immune system and malignant tumors. The study [8] is notable for analyzing a fractional-order cancer model that considers the interaction between stem cells and chemotherapy. The model was investigated using the Sumudu transform and the Atangana-Toufik numerical method. In [9], a tumour-immune model was qualitatively and quantitatively studied using the Atangana-Baleanu fractional operator. Numerical results were obtained using Lagrange piecewise interpolation for various fractional-fractional operators.
Due to their importance, mathematical models based on fractional derivatives are useful for investigating dynamic processes in various fields, such as engineering, biology [10, 11, 12, 13, 14, 15, 16], and especially epidemiology [17, 18, 19, 23]. The studies [20, 21] are dedicated to the investigation of deformation-relaxation anisotropic processes and taking into account the fractal structure in heat transfer processes in media. The algorithmic aspects of implementing such mathematical models, particularly the use of parallelization in the computational process, are presented in the research [21, 22]. These scientific publications present different studies of mathematical models for the qualitative and quantitative mathematical description of infectious diseases, taking into account fractional and fractional-fractal operators. Fractional mathematical models are also used to analyze therapy efficacy, particularly when addressing agent resistance and the body's response to treatment over time [24].
In this article, mathematical models based on Atangana-Baleanu derivatives, Caputo-Fabrizio derivatives, Liouville-Caputo derivatives were synthesized using fractal analysis methods, and software-algorithmic tools were developed to assess the impact of chemotherapy on the state of a cancerous tumour, considering the long-term memory effect. The Atangana-Toufik numerical scheme was used to implement the mathematical models numerically. Difference approximations of fractal operators have been constructed within mathematical models. One of the distinguishing aspects of our research was the analysis of the impact of the fractional differentiation parameter for these fractional derivatives within the studied mathematical models. 2. Formulation of the Mathematical Model and Basic Principles of
Fractal Operators We present the mathematical model for the growth dynamics of stem cells [2], which takes into account three populations: tumor cells T(t), immune system effector cells E(t), stem cells S(t), and the concentration of the chemotherapeutic substance M(t).
The dynamics of cancerous tumours during chemotherapy treatment are represented in [25] by the following system of equations: 푑 푑 = 푑 푑 푑 푑 = 1 − 푘 ( − − 휇 + 1 + 1) − 2( +
) = (1 − 푏 ) 푑 푑 =− 2
3 + + 푘 ( ) , (1)
0 = 0 , 0 = 0 , 0 = 0 , 0 = 0 .
In this model, the first equation describes the relationship between stem cells and the concentration of the chemotherapeutic agent. Stem cells can transform into specific cultured cells and lose their concentration over time at a rate of γ_1, k denotes the fraction of stem cells killed by chemotherapy. The second equation pertains to effector cells, which have a constant value α=α_1+α_2. Value α_1 is the natural baseline level of effector cells, and α_2 is the amount of effector cells produced from the transformation of stem cells. The parameter μ is defined as the natural death rate of effector cells. It is optimal for these cells when they are in a state of maturity, as this is when they produce effector cells [26]. The second term represents the mortality rate, which is proportional to the population of effector cells through μE. The third term represents the proliferation term p_1, which describes the process by which effector cells are stimulated by stem cells. In the equation, the terms p_2 ET_ and p_2 EM_ represent the interaction of effector cells with tumour cells and the chemotherapeutic drug at a rate of p_2. The first term, in the third equation (1), is the rate of production of tumour cells, and the second term represents the decay of tumour cells due to interaction with effector cells and the chemotherapeutic drug at rates p_3 and k_T, respectively. In conclusion, the final equation depicts the rate of change in the concentration of the chemotherapeutic drug. 2.1. Caputo Kernel, Riemann-Liouville Operator The Caputo kernel [27] is used to define fractional derivatives of the first-order Caputo derivatives. The Caputo kernel reflects the contribution of previous values of the function and its derivative in computation of the fractional derivative. A Caputo derivative is a modification of a RiemannLiouville derivative, which allows the use of standard initial conditions.
The Caputo kernel is defined as follows for a function f(x) and derivative order α (0 < α ≤ 1): − / (1 − ), 0 ≤
< 1 (τ), = 1 where Γ(·) denotes the gamma-function, and δ(·) is Dirac delta distribution.
The Caputo kernel κ t, α is used in the formula for calculating the Caputo derivative: where ψ τ is the ordinary derivative of the function ψ( τ ), and for T > 0 differentiable on ψ: [0, T] → C.
The Caputo kernel reflects the weight of preliminary values of function and its derivative in process of fractional differentiation. It considers both the prior values of this function and values of its derivative by integrating them with the corresponding weight. The fractional integral operator of Riemann-Liouville [27] is used to define fractional integrals, which are a generalization of ordinary (integer-order) integrals to fractional orders.
The Riemann-Liouville fractional integral operator has the form for T > 0 and an integrable function ψ: [0, T] - С: ∗ τ = 0 − , 푑 , (2) τ, τ = = τ, (τ, ) = (τ) = (τ, ) ∗
(τ) , −1/ , 0 < ≤ 1 where Γ(·) denotes the gamma function, and α is the fractional order of integration (0< α ≤ 1). 2.2. Caputo-Fabrizio Derivative, Liouville-Caputo Derivative, Atangana-Baleanu
Derivative The main formulas for determining fractional derivatives and their fractional integrals are presented here [28]. The derivative of (τ) with respect to Liouville-Caputo operator of order according to the formula: 1 − 0 (τ) = 0τ (τ − ) − −1 ( )푑 , τ > 0 where
For − 1 <
≤ , : ℝ+ → ℝ,
< 1, then the CF (Caputo-Fabrizio) fractional derivative has 0 τ τ = 1 − τ ' 푒 − where
is a normalization function (0) = 1 = 1.
The Caputo-Fabrizio (CF) fractional integral of the order 0 τ = τ 2 1 − 2 − + 2 − 2 τ − 1 −
푑 , of the function
(τ): τ 0 푑 , where ( ) = 2 , 0 < ≤ 1.
2−
Let ∈ ℍ1( , 푏), > 푏 and 0 < < 1, then the fractional derivative of Atangana-Baleanu in the sense of Liouville-Caputo (ABC) is given by the formula: 퐴 0 τ = 퐴 1 − τ − (τ − 1 − ) ' 푑 , − 1 < < , where 퐴 ( ) has the same properties as , that is, 퐴 ( ) = 1 − + represents the function as normalization, and 퐴 function [32].
We proceed to the ABC fractional integral of order and the convolution theorem: (0) = 퐴 (1) = 1.
Γ( ) (. ) is the one-parameter Mittag-Leffler using the inverse Laplace transform [28] AB 0 τ (τ) =
1− 퐴 ( ) (τ) + 퐴 ( )Γ( ) 0τ ( )(τ − ) −1푑 .
(9) 3. Main Material Presentation In order to improve the accuracy of the modelling of the effect of chemotherapy on the state of a cancerous tumour, we replace the integer order model (1) with a fractional order model. This is achieved using three of the most common fractional operators: Caputo, Caputo-Fabrizio and Atangana-Baleanu in the sense of Liouville-Caputo, which allow to analyse memory effects in (3) (4) (6) (7) (8) arbitrary order systems within the fractional calculus framework. For the numerical implementation of the fractional model of the effect of chemotherapy (1), the numerical scheme of TufekciAtangana [29, 30] is used. 3.1. Fractal Mathematical Model with Liouville-Caputo Derivative We write the mathematical model (1) in a form of the Liouville and Caputo fractional derivative as: 0
τ 0 τ = 0
τ − 휇 τ = τ + 1 1 τ = 1 − 푏
τ 0 τ τ − 푘 τ τ + 3 τ + τ + 푘 τ
τ τ τ τ , = 0,1,2, … , equation (11) can be transformed into a fractional integral equation [32]. A two-step Lagrange polynomial interpolation was employed [28, 31], the approximate solution of (11) takes the form: ∗ 퐴 , ,1 − ℎ −1, −1 ∗ 퐴 ,푘,2) , +1 =
0 + 퐴 , ,1 = 1 =0
(ℎ To simplify, we define the following expressions: , − ( − + 1) ( + + 2) − ( − ) ( + 2 − + 2) , 퐴 ,푘,2 = ( − + 1) +1 − ( − ) ( + − + 1)
, (10) (12) (13) (14) (15) (16) 3.2. Discretization of the Fractal Model with the Liouville and Caputo Derivative Using the proposed numerical scheme (12), we obtain the implementation algorithm for model (10).
S +1 = S0 + ∗ 퐴 , ,1− ℎ −1 ) ∗ 퐴 ,푘,2) ∗ 퐴 , ,1 − −1 , −1 ∗ 퐴 ,푘,2) τ + (17) (18) (19)
0 + − ℎ − ℎ Γ 3 Γ 4 ( where 1(τ, (τ), (τ), (τ), 2(τ, (τ), (τ), (τ), 3 (τ, (τ), (τ), (τ), 4(τ, (τ), (τ), (τ), = − 휇E(τ) + m=0 −1 , ℎ
3 −1 , m=0 −1 , (τ)): = (τ)): 0 τ τ = − 휇
τ + 0 τ τ = 1 − 푏
τ 3.3.
Fractal Mathematical Model with the Caputo-Fabrizio Derivative. The Atangana-Toufik Method We rewrite the model given by equation (1) in the form of the Caputo-Fabrizio fractional derivative based on equation (7). We obtain: 0 τ τ = τ − 푘
τ 0 τ 3.4. Finite Difference Approximations of Fractal Operators The initial value problem associated with the Caputo-Fabrizio operator (7) can be expressed as follows: 0 τ (τ) = (τ, (τ)), 0 < τ ≤ , 0 < ≤ 1, (20)
(0) = 0.
The aforementioned equation can be transformed into a fractional integral equation at the specified point τ = τ +1, = 0,1,2, …, by applying the fundamental theorem of fractional calculus: τ − 0 = 1 − τ −1, + 푑τ, (21)
A general numerical algorithm for implementing the Caputo-Fabrizio fractional derivative model is defined herein using the two-step Lagrange polynomial interpolation [28, 31].
− (27) + 퐴 1 ∗ 1(τ −1, + 퐴 1 ∗ 2(τ −1, + 퐴 1 ∗ 1(τ −1,
+ 퐴 1 ∗ 1(τ −1, 1(τ , −1(τ), 2(τ , −1(τ), 1( τ , −1(τ),
1(τ , −1(τ), 1 − +
3 ℎ ( ) 2 ( ) ,
퐴 = (τ), −1(τ), (τ), (τ),
(τ)) − −1(τ),
−1(τ), (τ), (τ), (τ),
( τ)) − −1(τ), −1(τ),
−1(τ), (τ), (τ), (τ),
(τ)) − −1(τ), (τ), −1(τ),
−1( τ), (τ),
(τ), −1( τ), 1 −
+ ( ) 2 −1( τ), −1( τ), ℎ ( )
.
(τ)) − 1 τ
τ τ + 1
τ − τ = − 휇 τ + τ + τ τ = 1 − 푏
τ In order to obtain the numerical solution of (19) in the sense of Caputo-Fabrizio-Caputo, the numerical scheme (22) is applied. The algorithm proposed by [26] is employed to obtain the numerical solution of (19) in the sense of Caputo-Fabrizio-Caputo.
+1 = − 퐴
∗ τ+1 = − 퐴
∗ +1 = − 퐴
∗ +1 = − 퐴
∗ 퐴 1 = 퐴 0 퐴 0 τ
τ where 3.6. Fractal Mathematical Model with Atangana and Baleanu Derivative in the
Sense of Liouville-Caputo The model presented in equation (1) is rewritten in the form of the Atangana–Baleanu fractional derivative in the sense of Liouville-Caputo (ABC), based on equation (9). The result is as follows: 0퐴 τ τ = 1 τ − 푘 τ τ with initial conditions:
0 τ = 0 , 0 τ = 0 , 0 τ = 0 , 0 τ = 0 .
In order to obtain a numerical solution to the fractional model (27), the Atangana-Toufik [29, 30] numerical scheme is employed. 3.7. Finite Difference Approximations of Fractal Operators The initial value problem associated with the Atangana-Baleanu operator in the sense of LiouvilleCaputo (9) can be expressed as follows: 0퐴 τ ( ) = , 0 < ≤ 1, (τ, (τ)), 0 < τ ≤ (28)
(0) = 0.
The aforementioned equation can be transformed into a fractional integral equation [32] by applying the fundamental theorem of fractional calculus [32] at the specified point τ = τ +1, = 0,1,2, … . By employing the two-step Lagrange polynomial interpolation of the function ( , ( )), we arrive at the following equation: yj+1L=et yu0s+deAn1BoC−teα:α Q τj, y τj + ABCα α k=j 0 ℎαQΓ ταk,+y 2τk 퐴 ,푘,1 − ℎαQ Γτk−α1 +,y2τk−1 퐴(2,푘9,)2 . 퐴 ,푘,1 = ( − 푘+1) ( + − 푘 + 2) − ( − 푘) ( + 2 − 푘 + 2) 퐴 ,푘,2 = ( − 푘 + 1) +1 − ( − 푘) ( + − 푘 + 1) 3.8. Discretization of the Fractal Model The proposed numerical scheme (29) allows us to derive the algorithm for the numerical implementation of model (27).
1 − 퐴 ,푘,2 퐴 퐴 (30) (31) (32) ∗ (3,3) ∗ (3, 4) (35) +1 = 0 + ∗ − ∗ − k=0 ℎ { +1 = S0 +
ℎ ∗
푘=0 퐴 (1 − ) ( )
{ ℎ {
− 휇Ek + − 휇Ek−1 + 퐴 − 2 Г(
푘 + + 2) − 휇E(τj) +
Sk1E+kS1k − 2 Γ( + 2) 1 EkSk S푘−1 + 1 − 2
Γ( + 2) 퐴 (1 − ) ( )
{ ℎ {
1 − 푏 k k=0 ℎ { 1 − 푏 k−1 − 1
−+ 1 E(τj)S(τ ) (S(τj) + 1) 퐴 ,푘,1 −
ℎ k + Mk Ek}
퐴 ,푘,1 k−1 + Mk−1 Ek−1} +
퐴 − 2 Г( 푘−1 + + 2)
∗ − 2( (τ ) + M(τj))E(τ )} + +1 = 0 + 1 − 푏 τj τj − 3 E τj + 푘 τj M τj
τj } + 퐴 ,푘,2 .
퐴 ,푘,2 .
k − 3Ek + 푘 kMk Γ( + 2) k−1 − 3Ek−1 + 푘 k−1 Mk−1 Γ( + 2) 1 − 퐴 ,푘,1
k−1} ∗ 푘=0 ℎ
퐴 − 2 푘 + Г( + 2) 퐴 ,푘,1 − + ℎ +
퐴 − 2 푘−1 + Г( + 2) ∗ 퐴 ,푘,2 4. Results and Discussion 4.1. Software Implementation Using the procedure described above, a software implementation of the discrete model with Liouville-Caputo fractional derivatives (10-13) has been carried out. The fundamental concept underlying the algorithm is the iterative calculation of values ( ) , , ( ) , ( ) at specified time intervals. The execution of an algorithm involves the following steps: 1. as outlined in the analysis of the article [26]: S0 = 1, E0 = 1, T0 = 1, ℽ1 = -0.02825, b = 10-9, 푘 = 1, p1 = 0.1245, r = 0.18, p2 = 1, 푘 = 0.9, p3 = 0.9, ℽ2 = 6.4, V(τ) = 1.
Figure 1 presents the UML diagram of the software application. Let's describe each of the classes.
Base Algorithm: This is the base class that represents the general algorithm for computing the mathematical model. When this class is created, the parameter values and the initial conditions of the system are initialized. The equations of the mathematical model are also presented in this class as methods.
GUI: This class is responsible for interaction with the user interface. It receives as input the initial values T(0), E(0), S(0) , and M(0), as well as a number of parameters such as y1, y2, , , p1, p2, p3, r, b, 푘 , 푘 , V(t). Upon clicking the "Run Simulation" button, it initiates the computation of the selected algorithms of the corresponding classes. This class is also responsible for displaying the results in a graphical representation.
Riemann Liouville Algorithm: This class, like the previous one, solves the fractional-order Liouville-Caputo problem based on the Atangana-Toufik method.
Runge Kutta Algorithm: This class is responsible for computing the integer-order results of the mathematical model using the Runge-Kutta method. This algorithm is applied in the software implementation of the mathematical model to compare the integer-order and fractional-order results. Numerical experiments on the effect of chemotherapy on the state of a cancer tumour were carried out, taking into account long-term memory in the context of comparing the use of fractional Caputo, Riemann-Liouville, Caputo-Fabrizio, Atangana-Baleana, using the above algorithms and developed software for their implementation. The functions of tumour cells T(τ), stem cells S(τ), effector cells of the immune system E( τ ), and the concentration of the chemotherapeutic agent M( τ ) were studied as functions of time and the fractional differentiation parameter . The graphical results for the Atangana and Baleanu model in the sense of ABC (Liouville-Caputo) (32 – 35) are presented in Figures 2–5 for various values of = 1, 0.98, 0.96, 0.92 , while the results for the CaputoFabrizio model (23 – 26) are shown in figure 2 for different values of = 1, 0.98, 0.96, 0.92 . Figure 2 presents an enlarged time scale for a more detailed understanding of the differences in the behaviour of functions under different fractional orders.
Figures 7–10 provide a comparative analysis based on the fractional parameter = 0.94 for the studied variables: the concentration of tumour cells T(τ) (yellow line), effector cells of the immune system E( τ ) (red line), stem cells S( τ ) (blue line), and the chemotherapeutic agent M( τ ) (purple line). These results are presented for models using the Caputo, Riemann-Liouville, AtanganaBaleanu, and Caputo-Fabrizio fractional derivatives.
It is noteworthy that the graphical dependencies of the Liouville-Caputo fractional model and the Atangana-Baleanu model closely approximate one another in the sense of Liouville-Caputo. The analysis of the obtained data indicates that the concentration of cancer cells declines over time, approaching zero at a gradual rate. The cell death rate under the specified initial conditions is p3=0.9. Rate of growth of tumor cells was found to be lower than the rate of their interaction with effector cells, with the help of stem cells and chemotherapy drugs. This constitutes a principal outcome of the study, indicating that the immune system undergoes modifications. Therefore, the combination of stem cell therapy with chemotherapy offers promising prospects for optimizing cancer treatment and improving the quality of life for patients. 5. Conclusion In this work, mathematical models of the effects of chemotherapy on the state of a cancerous tumour were synthesized using fractal analysis methods, taking into account the long-term memory effect based on Liouville-Caputo derivatives, Atangana-Baleanu derivatives, Caputo-Fabrizio derivatives. Difference approximations of the aforementioned fractal operators in mathematical models were presented and constructed. The numerical algorithms for implementing fractional mathematical models of chemotherapy effects were developed based on the Atangana-Toufik scheme and an adaptation of two-step Lagrange polynomial interpolation. One of the key stages of the work was the software implementation of the model, the development of the interface, and the visualization of the results. The numerical results are presented in the form of graphical illustrations and were obtained using software.
This study investigated the influence of the fractional differentiation parameter on the time dynamics of tumour cells, immune system effector cells, stem cells, and the chemotherapeutic agent concentration in models based on Liouville-Caputo, Atangana-Baleanu derivatives and CaputoFabrizio.
A comparative analysis was conducted on the dependency of the fractional parameter on the investigated quantities, including tumour cells, immune system effector cells, stem cells, and the concentration of chemotherapeutic agents, for models utilizing Riemann-Liouville, Caputo, Atangana-Baleanu and Caputo-Fabrizio fractional derivatives. Additionally, an analysis was conducted to ascertain the degree of correlation between the results of the fractal models and those of the classical integer-order model.
The results obtained from fractional modelling demonstrate a significant influence of historical factors on the temporal evolution of tumour cell, immune system effector cell and stem cell concentrations. To predict the effect of chemotherapy on the state of a cancerous tumour, the presence of a fractional time derivative in the models as a parameter of the time derivative is important. 6. Declaration on Generative AI
During the preparation of this paper, the authors utilized Grammarly to verify spelling and grammar accuracy. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the publication’s content. 7. References [1] Anderson ARA. A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion. IMA Math. App. Med. Biol. 2005; 22:163–186. [2] Chaplain MAJ, McDougall SR, Anderson ARA. Mathematical modeling of tumor-induced angiogenesis. Annu. Rev. Biomed. Eng. 2006; 8:233–257 [3] Ribba B, Marron K, Agur Z, Alarcón T, Maini PK. A mathematical model of Doxorubicin treatment efficacy for non–Hodgkin’s lymphoma: Investigation of the current protocol through theoretical modelling results. Bull. Math. Biol. 2005; 67:79–99. [4] Friedman A. Lecture Notes in Mathematics. Vol. 1872. Springer Verlag; 2006. Tutorials in Mathematical Biosciences. III. [5] Farayola, M.F.; Shafie, S.; Siam, F.M.; Khan, I. Mathematical modeling of radiotherapy cancer treatment using Caputo fractional derivative. Comput. Methods Programs Biomed. 2020, 188, 105306. [6] Іdrees, M.; Alnahdi, A.S.; Jeelani, M.B. Mathematical Modeling of Breast Cancer Based on the Caputo–Fabrizio Fractal-Fractional Derivative. Fractal Fract. 2023, 7, 805. [7] Jordão, G.; Tavares, J.N. Mathematical models in cancer therapy. Biosystems 2017, 162, 12– 23 [8] Shabir Ahmad, Aman Ullah, Thabet Abdeljawad, Ali Akgül, Nabil Mlaiki,