Mathematical analysis and numerical simulation for fractal-fractional cancer model

1.   Introduction

Cancer is a word used to describe disorders in which aberrant cells divide uncontrollably and can infiltrate neighboring tissues. According to the World Health Organization (2020), cancer is the second leading cause of mortality globally, accounting for approximately one in every six deaths ^[1]. Since the middle of the 1960s, mathematical modeling and nonlinear simulations of the tumor growth process has been researched due to the significant public health issues and the requirement for immediate health measures ^{[2,3,4,5,6,7,8,9]}.

Edward Lorenz, a meteorologist and mathematician, discovered the chaos phenomenon in the unpredictable and irregular behavior of nonlinear dynamical systems in 1963 ^[10]. Chaos can be expressed mathematically via deterministic iterations of nonlinear difference equations or the development of nonlinear ordinary differential equations (ODEs) or partial differential equations (PDEs). The study of chaotic systems has been heralded as one of the most significant scientific accomplishments of the twentieth century. While the field is still in its infancy, there is no doubt that it is becoming increasingly important in various of scientific disciplines. To that end, chaos has been demonstrated to exist in a wide range of systems, including electronics ^[11], chemistry ^[12], economics and finance ^[13,14], biological systems ^{[15,16,17,18]} and so on.

It is worth noting that fractional calculus is a vital branch of mathematics. Because of the memory and genetic peculiarity of fractional-order differential equations, several researchers have modeled biological phenomena using fractional calculus derivatives. As a result, it is a very useful tool for describing genuine natural processes. Many papers on fractional-order dynamical models have recently been published ^[19,20,21].

Atangana ^[22] presented a new advanced type of fractal fractional derivative in 2017, bridging the gap between fractional and fractal calculus. Fractal-fractional operators contain two components: the fractional order and the fractal dimension (order). Differential equations using the fractal-fractional derivative transform the assumed system's order and dimension into a rational order system. The major goal of defining these derivatives is to examine fractal nonlocal boundary and initial value problems in nature. Certain mathematicians developed various results and designed some fractal-fractional models that exhibit improved simulations for representing mathematical structures in this direction ^{[23,24,25,26]}.

The nonstandard finite deferential numerical methods were first introduced by Mickens in 1994 ^[27]. These methods are well known for maintaining the positivity, boundedness and stability of nonlinear systems' equilibrium points ^[27,28].

In the paper ^[15], authors introduced and on studied the following three-dimensional order cancer model

where $x(t)$ stands for the number of tumor cells at time $t \, y(t)$ for the number of healthy host cells at time $t$, and $z(t)$ for the number of effector immune cells present at time $t$ within the single tumor-site compartment, and $x_{0}, \, y_{0}\, \text{ and }z_{0}$ are the associated initial values of system (1.1). Here, the parameters $a, \, b$ and $c$ are positive real numbers that indicate the growth rates of populations of $x(t), y(t)$ and $z({t})$, respectively. If $\alpha_{i} = \alpha$ for every $i = 1, ..., 3,$ then the system (1.3) is called commensurate order; otherwise, it is named incommensurate order ^[29]. The fractional version of system (1.1) was considered in the paper ^[16], and were described by

where ${ }_{0}^{C} D^{\alpha}$ is the ${\alpha}$-order Caputo differential operator.

The three-dimensional fractal-fractional-order cancer model is the main topic of this research:

where ${ }_{0}^{C} D^{\alpha, \beta}$ is the $({\alpha, \beta})$ fractal-fractional-order Caputo differential operator.

The rest of this paper is organized as follows. Section 2 provides some fundamental definitions of generalized Banach spaces in the sense of Perov, its properties and fractal fractional operators in the Caputo sense. Section 3 is devoted to the existence and the uniqueness with Perov fixed point theorem. In Section 4, the suggested model's numerical solution was achieved using the Grünwald-Letnikov nonstandard finite difference scheme of Caputo derivative (in short GL-NSFDM) scheme using MATLAB software. Section 5 presents the Lyapunov's stability of the equilibrium points of the proposed system by varying the fractional order and the set of parameter $(a, b, c)$, and by maintaining the fractal dimension $(\beta_{1}, \, \beta_{2}, \, \beta_{3}) = (1, 1, 1)$. Section 6 shows the Ulam-Hyers stability of the Incommensurate fractal-fractional-order cancer model (1.3). Finally, the discussion and the conclusion are given in the last two sections.

2.   Preliminary

We present some basic notation, results of generalized Banach spaces in the sense of Perov, matrices converges to zero and Fractal-Fractional calculus in Caputo sense, which will be essential in the next sections. We begin with defining on $\mathcal{M}_{m\times n}(\mathbb{R}_{+})$ the partial order relation as follow: $\text{Let } \varLambda, \; \varUpsilon \in \mathcal{M}_{m\times n}(\mathbb{R}_{+}), m\geq 1\text{ and } n\geq 1. \text{ Put } \varLambda = (\varLambda_{i, j})_{1\leq j\leq m \atop 1\leq i\leq n} \text{ and } \varUpsilon = (\varUpsilon_{i, j})_{1\leq j\leq m \atop 1\leq i\leq n}.$ Then,

and we write $I_{n}$ for the identity $n \times n$ matrix and $O_{n}$ for the zero $n \times n$ matrix.

Definition 2.1. Let $\mathcal{E}$ be a vector space over $\mathbb{K} = \mathbb{R} \text{ or } \mathbb{C}$. A generalized norm on $\mathcal{E}$ is a map

has the next properties

$(i)$ For all $\vartheta\in \mathcal{E}$; if $\| \vartheta\|_{G} = 0_{ \mathbb{R}_{+}^{n}}$, then $\vartheta = 0_{\mathcal{E}}$,

$(ii) \; \|a \vartheta\|_{G} = |a|\| \vartheta\|_{G}$ for all $\vartheta\in \mathcal{E}$ and $a \in \mathbb{K}$, and

$(iii) \; \| \vartheta+ \omega \|_{G}\preccurlyeq \| \vartheta\|_{G} + \|\omega\|_{G}$ for all $\vartheta, \omega \in \mathcal{E}.$

The pair $(\mathcal{E}, \|\cdot\|_{G})$ is called a generalized normed space. Moreover, $(\mathcal{E}, \|\cdot\|_{G})$ is called a generalized Banach space (in short, GBS), if the vector-valued metric space generated by its vector-valued metric $\delta_{G}(x, y) = \| x- y\|_{G}$ is complete.

Let $(\mathcal{E}, \|\cdot\|_{G})$ be a generalized Banach space. In the rest of this article for $r = (r_{1}, \cdots, r_{n}) \in \mathbb{R}^{n}_{+}$, $\vartheta_{0}\in \mathcal{E}$ and $i = 1, \cdots, n$, we denote by:

for the open ball centered at $\vartheta_{0}$ with radius $r$, and by:

for the closed ball centered at $\vartheta_{0}$ with radius $r$. If $\vartheta_{0} = 0$ we simply denote $B_{r} = B({0}, r)$ and $\bar{B_{r}} = \bar{B}({0}, r)$. Finally, we respectively denote by $\overline{\mathcal{K}}$ and ${co}(\mathcal{K})$ for the closure and the convex hull of a subset $\mathcal{K}$ of $\mathcal{E}$.

Definition 2.2. A matrix $\varUpsilon\in \mathcal{M}_{n\times n}(\mathbb{R_{+}})$ is said to be convergent to zero if

Lemma 2.3. ^[30] Let $\varUpsilon\in \mathcal{M}_{n\times n}(\mathbb{R}_{+})$. The following assertions are equivalent:

$(i) \; \varUpsilon^{m} \longrightarrow O_{n}, \quad {\text{as}}\quad m\longrightarrow \infty$.

$(ii)$ The matrix $I_{n} -\varUpsilon$ is invertible, and $(I_{n} -\varUpsilon)^{-1}\in \mathcal{M}_{n\times n}(\mathbb{R}_{+})$.

$(iii)$ The spectral radius of $\varUpsilon$ is strictly less than 1.

Definition 2.4. Let $(\mathcal{E}, \delta_{G})$ be a generalized metric space and $N$ be an operator from $\mathcal{E}$ into itself. $N$ is called $\varUpsilon$-contraction with matrix $\varUpsilon\in \mathcal{M}_{n\times n}(\mathbb{R}_{+})$ that is converges to $O_{n}$, if $\text{for all } \varrho, \; v \in \mathcal{E}$ we have

In the following, an extension of the Banach contraction principle by Perov is given.

Theorem 2.5. ^[31] Let $\mathcal{E}$ be a complete generalized metric space and let $N: \mathcal{E} \longrightarrow \mathcal{E}$ be an $M$-contraction operator. Then, $N$ has a unique fixed point in $\mathcal{E}.$

Next, we give some important concepts from fractal-fractional calculus in Caputo sense. We refer the reader for the reference ^[32] for more details.

Definition 2.6. Let $\varrho$ be differentiable in opened interval $(a, b)$, if $\varrho$ is fractal differentiable on $(a, b)$ with order $\beta$, then the FF-derivative of $\varrho$ of order $\alpha$ in the Caputo sense with power law is given as:

where

Lemma 2.7. The Eq (2.1) can be written as follows:

3.   Existence and uniqueness results

Lemma 3.1. $(x, y, z)$ is a solution of the fractal-fractional-order system (1.3), if and only if it is a solution of the following problem

Theorem 3.2. Suppose that there is a vector with positive entries fulfills

in addition, if the matrix

converges to $O_{3}$, where $\mathcal{H}(\alpha_{i}, \beta_{i})$ denotes the beta function of $\alpha_{i}$ and $\beta_{i}$, then the system (3.1) has a unique solution in the space $\mathcal{C}([0, \, T])\times\mathcal{C}([0, \, T])\times\mathcal{C}([0, \, T]).$

Proof. Let $\mathcal{K}$ be the closed ball $\bar{B}((x_0, y_0, z_0), \varUpsilon)$ on $\mathcal{E} = \mathcal{C}([0, T], \mathbb{R})\times \mathcal{C}([0, T], \mathbb{R})\times \mathcal{C}([0, T], \mathbb{R})$ centered at $(x_0, y_0, z_0)$ of radius $\varUpsilon\succ0_{ \mathbb{R}^{3}_{+}}$ where $r$ satisfies the above inequality in (3.2). We recall that the space ${\mathcal{E}} = \mathcal{C}([0, T], \mathbb{R})\times \mathcal{C}([0, T], \mathbb{R})\times \mathcal{C}([0, T], \mathbb{R})$ is generalized Banach space endowed with the generalized norm

The proof will be broken up into several steps.

Step 1: First, we shall show that the mapping

is G-contraction where $N$ defined by the following formula:

To this end, let $X_1 = (x_1, y_1, z_1), \ X_2 = (x_2, y_2, z_2)\in \mathcal{E}$ and for $t\in [0, \, T]$ we have

And we have for all $A_{1}, \, A_{2}, \, B_{1}, \, B_{2} \in \mathbb{R}$

then

by taking the supermum over $t$ we find

And we have for all positive numbers $\varrho_{1}, \, \varrho_{2}$ and $\gamma\geq1$

then,

It is clear that

then

and

As conclusion

where

Step 2: Our objective here is to prove that the operator $N$ maps $\mathcal{K}$ into itself. To do so, let $X = (x, y, z), \ \in \mathcal{K}$ and for $t\in [0, \, T]$ we have

Then

By the same manner, we find

and

Hence,

By using Perov fixed point Theorem 2.5, we conclude that the system (3.1) has a unique solution in $\mathcal{K}$. □

4.   Numerical method for solving fractal-fractional cancer model

According to Lemma 2.7, the system (1.3) can be written as follows:

The discretization of fractional derivative is given by GL approach ^[33,34]:

where $t_{n} = n \Delta t, \quad \Delta t = \frac{T}{N}$ is the time-step size, $N$ is a natural number, $\mu_{j_i} = (-1)^{i-1}\left(\begin{array}{c}\alpha_{j} \\ i\end{array}\right), \mu_{j_1} = \alpha_{j}, q_{j_i} = \dfrac{i^{\alpha_{j}}}{\Gamma(1-\alpha_{j})}$ and $i = 1, 2, \ldots, n+1, \; j = 1, 2, 3.$ In addition, let us assume that ^[35]:

Using the GL approximation and the NSFD framework ^[27], we discretize the first equation in (1.3) as follows:

where,

From the first equation in (4.1), we have:

hence

Looking that the function $F_{1}$ can be written as next:

By substituting this latter in (4.2), and using the fact that the nonlinear term $g_{1_1}(x_{n}, y_{n}, z_{n})(t_{n})+x_{n}g_{1_2}(x_{n}, y_{n}, z_{n})(t_{n})$ is approximated by $g_{1_1}(x_{n}, y_{n}, z_{n})(t_{n})+x_{n+1}g_{1_2}(x_{n}, y_{n}, z_{n})(t_{n})$ in a nonlocal way, we find that:

Repeating the same procedure to the second and the third equation of the system (1.3), we conclude that the discretization of system (1.3) using GL-NSFDM can be formulated as follows:

5.   Lyapunov's stability of the incommensurate fractional-order cancer model

In this section we analyze the dynamics of the incommensurate by taking the initial conditions $(x_{0}, \, y_{0}, \, z_{0}) = (0.4, 0.5, 0.5)$, $(\beta_{1}, \, \beta_{2}, \, \beta_{3}) = (1, 1, 1)$ and by selecting different values of the fractional-orders $\alpha_{1}, \, \alpha_{2}, \; \alpha_{3}$ and varying the set of parameter $(a, b, c).$

Definition 5.1. ^[36] The equilibrium point $E$ is called a saddle point of index one (two) if the Jacobian matrix evaluated at point $E$ has exactly one (two) eigenvalue with non-negative real part. Scrolls are generally created only around the saddle points of index two.

The Following Lemma gives the sufficient condition to exhibit the equilibrium point $E$ a stability nature.

Lemma 5.2. ^[37] The equilibrium point $E$ of the fractional-order system is locally asymptotically stable in the Lyapunov sense if the following condition is satisfied:

where $\Delta(\lambda) = det(J-\operatorname{diag}(\lambda^{\delta \alpha_{1}}, \, ..., \, \lambda^{\delta \alpha_{k}}))$ and $J = (\dfrac{\partial F_{i}}{\partial x_{j}}) _{i, j = 1, ..., k}$ is the Jacobian matrix evaluated at $E.$ The parameter $\delta$ is the least common multiple of the denominators $q_{i}$s of $\alpha _{i}$s, where $\alpha _{i} = \dfrac{p_{i}}{q_{i}}$, $(p_{i}, q_{i}) = 1$, $p_{i}, q_{i}\in \mathbb{Z} ^{+}$.

If the condition (5.1) does not satisfy, we are in the following state.

Lemma 5.3. ^[38] A necessary condition for fractional-order system to exhibit the chaotic attractor is

Furthermore, the number ${\pi }/{2\delta} -\min_{\lambda\in\{\lambda:\, \Delta(\lambda) = 0\}}|\arg(\lambda)|$ is called the instability measure for equilibrium points in fractional order systems (in short IMFOS).

The Jacobian matrix of system (1.3) is

In ^[39], the authors established that the system (1.3) has five real equilibrium points, where four of them are obtained analytically and can be described as follows:

$1) \; E_{0} = (0, 0, 0),$

$2) \; E_{1} = (0, -1, (\frac{b}{b - 1})), \text{ if } b \neq 1$,

$3) \; E_{2} = ((\frac{c}{c - 1}), 0, -1), \text{ if } c \neq 1$,

$4) \; E_{3} = (-1, (\frac{a}{a - 1}), 0), \text{ if } a \neq 1$,

Because they have negative coordinates, the equilibrium points $E_{1}, \, E_{2},$ and $E_{3}$ are irrelevant to the ensuing dynamics (negative populations are not defined and, consequently, the dynamics must take place in the positive octant). The equilibrium point $E_{0}$ relates to a situation in which there is no cell at all. The fifth equilibrium point changes according to the set of parameters $(a, b, c)$. The following Table 1 gives the index of saddle points (ISP), and the IMFOSs of three sets selected parameters and different fractional-orders.

According to Table 1 the IMFOSs for the equilibrium points $E_{*}$ and $E_{**}$ are negative numbers, which implies that $E_{*}$ and $E_{**}$ are stables. Therefore, for the given derivative orders, the systems in case 1 and case 2 do not have the necessary condition to exhibit chaos. Numerical simulation results in Figures 1 and 2, respectively, confirm this conclusion. In the third case, Table 1 shows that the IMFOS is non-negative number, and the equilibrium point $E_{***}$ is a saddle point of index 2. This implies that the system (1.3) in case 3 in The Table 1 satisfies the necessary condition for exhibiting a 1-scroll attractor. As shown in Figure 3, numerical simulation results confirm this conclusion.

6.   Ulam-Hyers stability of the incommensurate fractal-fractional-order cancer model

Here, we are going to demonstrate the stability of Ulam-Hyers sense of the proposed model. We adopt the following definitions from ^[40].

Definition 6.1. Let $(X, d_{G})$ be a generalized metric space and $F: X \rightarrow X$ be an operator. Then, the fixed point equation

is said to be generalized Ulam-Hyers stable if there exists an increasing function $\psi: \mathbb{R}_{+}^{m} \rightarrow \mathbb{R}_{+}^{m}$, continuous in $0_{ \mathbb{R} ^{m}}$ with $\psi(0) = 0$, such that, for any $\varepsilon: = \left(\varepsilon_{1}, \ldots, \varepsilon_{m}\right)$ with $\varepsilon_{i} > 0$ for $i \in\{1, \ldots, m\}$ and any solution $Y^{*} \in X$ of the inequalities

there exists a solution $X^{*}$ of (6.1) such that

Consider a small perturbation $\Phi: = (\Phi_{1}, \Phi_{2}, \Phi_{3}) \in \mathcal{C}([0, {T}])\times\mathcal{C}([0, {T}])\times\mathcal{C}([0, {T}])$ such that $\Phi(0_{ \mathbb{R}^{3}}) = 0_{ \mathbb{R}^{3}}$. Let

● $|\Phi_{i}(t)| \leq \varepsilon_{{i}}$, for $\varepsilon_{{i}} > 0 \; i = 1, ..3.$

●

Lemma 6.2. The solution of the perturbed model

fulfills the relation given below

Proof. The solution of (6.3) is given by

Then, we have

Repeating the same procedure to the second and the third equations of the system (6.3), we have

Hence, the proof is completed. □

Theorem 6.3. If the matrix $\varTheta$ (3.3) converges to $O_{3}$, then (1.3) is generalized Ulam-Hyers stable.

Proof. Let $X = (x, y, z)$ be any solution of the inequality (6.4), and let $X^{*} = (x^{*}, y^{*}, z^{*})$ be the unique solution of (1.3), then

By the same manner, we find

and

Consequently, one can write

Since the matrix $\varTheta$ converges to zero, then we have

Hence, the solution of the proposed problem is generalized Ulam-Hyers stable. □

Example 6.4. Consider the following fractal-fractional-order cancer model

Note that for $\beta = (1, 1, 1)$ the system (6.6) was stated in the third case in Table 1 and exhibit a chaotic behavior, as it have shown in Figure 3. By taking $\varUpsilon = (1, 1, 1)$ we find

Furthermore,

Then, the eigenvalues of matrix $\varTheta$ are as follows:

hence, the spectral radius of $\varTheta$ is $\rho = 1.6069\times 10^{-2} < 1.$ As a result, the system (6.6) has a unique solution that is generalized Ulam-Hyers stable according Theorem 3.2 and Theorem 6.3. Figure 4 illustrates the conclusion.

7.   Discussion

This section is devoted to numerical simulations of the proposed model under investigation in the present paper. As described in Table 1 and Example 6.4, the approximate solutions of the fractal-fractional system (1.3) are given in Figures 1–4 with varying values of fractional-order parameters $(\alpha _{1}, \alpha _{2}, \alpha _{3})$, the fractal dimension $(\beta _{1}, \beta _{2}, \beta _{3})$ and the set parameters $(a, b, c)$. We briefly presented the simulation of this model using the GL-NSFDM numerical method as given in (4.3). The time interval is $[0, 3000]$ and $\phi(t) = exp (t) -1$. MATLAB computer language was used to accomplish all computations in this work.

By maintaining the fractal dimension $(\beta_{1}, \, \beta_{2}, \, \beta_{3}) = (1, 1, 1)$, the phase plots in Figures 1 and 2 presented above indicated that the system exhibits a non-chaotic behavior, counter to Figure 3, where the phase portrait shows a chaotic behavior. By varying the fractal dimension of the previous case to $\beta = (0.88, 0.8, 0.9)$, the fractal-fractional cancer model widths a stable behavior in Figure 4. This demonstrates how the fractal dimension -which is absent in both the fractional and classical models- can turn the behavior of the solutions from chaotic into a stable state and vice versa. By returning to the proposed model, we observe that the parameters $(a, b, c)$ are constants. However, from the biological point of view, these parameters may show a randomness behavior; that is a limitation of our study. In future research, we are focusing on replacing these parameters with Ornstein-Uhlenbeck process to make the model more realistic.

8.   Conclusions

A mathematical study of the growth of tumor has been discussed in this paper. The contribution is based on describing the cancer process by a novel fractal-fractional order model. This is inspired by population dynamics and contains terms that refer to tumor cells, effector immune cells and healthy tissue cells. The study in ^[16] is a special case from the present paper. Perov's fixed point theorem showed the existence and the uniqueness result. Besides, the numerical simulations were received with the Grünwald-Letnikovv nonstandard finite difference scheme. The dynamics of the proposed Incommensurate fractal-fractional cancer model were analyzed by varying the value of the fractional order, the fractal dimension and the values of the system parameters. The obtained results are also compatible with theoretical analysis. The proposed model could describe a wide range of biologically observed tumor states, including stable and chaotic states.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

All sources of funding of the study must be disclosed

Conflict of interest

The authors declare there is no conflict of interest.