COMPARISON OF DYNAMICAL BEHAVIOR BETWEEN FRACTIONAL ORDER DELAYED AND DISCRETE CONFORMABLE FRACTIONAL ORDER TUMOR-IMMUNE SYSTEM

In this paper, we analyze the dynamical behavior of the delayed fractional-order tumor model with Caputo sense and discretized conformable fractional-order tumor model. The model is constituted with the group of nonlinear differential equations having effector and tumor cells. First of all, stability and bifurcation analysis of the delayed fractional-order tumor model in the sense of Caputo fractional derivative is studied, and the existence of Hopf bifurcation depending on the time delay parameter is proved by using center manifold and bifurcation theory. Applying the discretization process based on using the piecewise constant arguments to the conformable version of the model gives a two-dimensional discrete system. Stability and Neimark–Sacker bifurcation analysis of the discrete system are demonstrated using the Schur-Cohn criterion and projection method. This study reveals that the delay parameter τ in the model with Caputo fractional derivative and the discretization parameter h in the discrete-time conformable fractional-order model have similar effects on the dynamical behavior of corresponding systems. Moreover, the effect of the order of fractional derivative on the dynamical behavior of the systems is discussed. Finally, all results obtained are interpreted biologically, and numerical simulations are presented to illustrate and support theoretical results. Mathematics Subject Classification. 34A08, 92D25. Received July 1, 2020. Accepted December 7, 2020.

Kuznetsov suggested the following mathematical model: where E and T represent the number of effector cells (ECs) and tumor cells (TCs) respectively, F (E, T ) is the recruitment function which represents the stimulating effect of tumor cells on immune cell proliferation. Biological meanings of other parameters are given in Table 1. Moreover, Kuznetsov et al. suggested that F should be F (E, T ) = pET r + T , and then, they studied the dimensionless form of model (1.1), that is: where x and y denote the dimensionless density of ECs and TCs respectively, x = E/E 0 , y = T /T 0 , σ = s nE0T0 , ω = p nT0 , η = r T0 , µ = m n , δ = d nT0 , γ = a nT0 , β = bT 0 , t = nT 0 t with initial conditions x(0) = x 0 ≥ 0 and y(0) = y 0 ≥ 0.
Besides the above ordinary differential equation model, there is a growing interest in the effects of time delay on modeling real-world phenomena, especially in biological dynamical systems [42,50], as time delay is a suitable tool for modeling dynamical systems that also depend on past data [9]. Time delays incorporated into biological systems are used to reflect some crucial biological facts such as regeneration time, maturation period, reaction time, feeding time, and immunity period. It is well known that the tumor-immune cell interaction is very complicated because competition between them is multi-dimensional, reciprocal, and includes many different processes and dissimilar pathways [48]. Therefore, we think that the delayed version of the model (1.2) can reflect the complex behavior of this interaction better than the undelayed counterpart. Here, the term ωxy η + y is called recruitment function and represents the stimulating effect of the tumor on the proliferation of immune cells. The immune response to the appearance of tumor cells requires the transmission of signals that initiate the proliferation of suitable immune agents and the production of specific cells and proteins. To reflect this phenomenon, the discrete-time lag τ is added to the recruitment function [22,25,48]. In [24,48], this delay is only added to the tumor cells variable y, not in the effector cells variable x, since one should consider that process in the context of per capita (x /x) change. This approach here is the same as in the derivation of the logistic model with delay by Hutchinson [29]. So, we replace the term ωxy η + y in system (1.2) by As tumor cells have a stimulating effect on effector cells, they also have a lethal effect, which is represented by the term −µxy in model (1.2). Nevertheless, this deadly effect is reciprocal. Hence, the term −xy in the second equation represents the defense of effector cells against tumor cells. The dynamics of this mechanism in vivo is complicated and not totally understood. But, we can still say there are two major branches of this mechanism: humoral immunity mediated by antibodies produced by B lymphocytes and cell-mediated immunity mediated by T lymphocytes [21]. When effector cells form tumor-immune complexes, they emanate chemokines that activate more effector cells around the malignant tumor site and destroy them. But, this is not an instantaneous process, as it takes time to detect the expression of some cytokines [48]. So, we replace the term xy in system (1.2) by x(t − τ )y in system (1.3). This delay is only presented in the effector cells x, not in the tumor cells y. Because we want to consider this process in the context of per capita (y /y) change. The method of adding discretetime lag to only one variable in the competition term has been investigated in articles on the Lotka-Volterra prey-predator system [44,54,58] and these discrete-time lags are called pure delays [44].
On the other hand, as mentioned above, tumor cells have the ability to make the immune system ineffective. This mechanism also includes various branches and processes such as local immune suppression, induction of tolerance, and systemic dysfunction in T-cell signaling. Tumor cells can also utilize several distinct pathways to actively prevent destruction by immune cells. This mechanism is also not simultaneously obtained. So, we add the discrete-time delay τ into the system to better reflect this phenomenon. Using the same approach described above, we are adding a pure delay τ only into the tumor cells variable y, and we obtain µxy(t − τ ) in system (1.3).
The fractional derivative operator is non-local, and, consequently, fractional differential equations provide a good tool for the description of memory and hereditary properties of various processes [47], and more accurate results can be obtained when compared with integer-order counterparts [6,11,32,47]. Therefore fractional calculus has been widely applied to various fields of science and engineering such as boundary value problem [32], Vander Pol oscillator [52], harmonic oscillator [7], human liver modeling [8], plasma physics [27], exothermic reactions [39], food web [36], fish farm [53], epidemic diseases [31,37,38], drug diffusion [20,30], tumor-immune interaction [11,12,51], etc. In [11], the authors considered a fractional Gompertz model and showed that the Gompertz model of fractional-order 0.68 produced a better fit to the experimental dataset than the integerorder Gompertz model. In [12], the authors proposed a fractional-order model of brain tumor growth. In [51], the authors provided an extended version of the Kirschner and Panetta model of tumor-immune interaction with time lags and fractional-order. Based on the above discussions, we also consider fractional-order tumor models with respect to two definitions of fractional derivative. A fractional-order version of the Kuznetsov tumor-immune interaction model (1) is given as where D α represents the Caputo fractional derivative operator of order α with α ∈ (0, 1). The Caputo fractional derivative of order α > 0 of a real-valued function f is defined as where m is an integer and m − 1 < α < m.
Then, by replacing the left-hand side derivatives of the delayed model (1.3) with Caputo sense fractional derivatives, we obtain the following model: where α ∈ (0, 1) is the order of the Caputo fractional derivative, τ ≥ 0, The most known fractional derivative definitions are Riemann-Liouville, Caputo, and Grünwald-Letnikov [49]. Recently, new approaches such as conformable fractional derivative, Atangana-Baleanu fractional derivative have also been proposed and applied to various fields [5,13,33]. The left conformable fractional derivative starting from a of the function f : [a, ∞) → ∞ of order 0 < α ≤ 1 defined as follows: provided the limit exist. The right fractional derivative of order 0 < α ≤ 1 terminating at b of f is defined by We also note that if f is differentiable in the ordinary sense, we have the following properties: The above definition has some similarities with the ordinary derivative and has been proposed to eliminate difficulties arising from Caputo fractional derivative in applications [35]. The conformable fractional derivative is perhaps not a fractional derivative, as in the Caputo sense, but has a fractional-order and is a natural generalization of the classical derivative. Some physical and biological applications of the conformable derivative were studied in [3,5,13,34].
equations by conformable fractional derivatives. Secondly, we replace delayed terms in the system (1.5) by PCA, and, consequently, we have the following model: denotes the integer part of t ∈ [0, ∞) and h > 0 is the discretization parameter. Instead of defining a new parameter to the time-delay effect like in [16,28], we use the discretization parameter h to control the magnitude of the time-delay. Here, the biological interpretation of piecewise-constant arguments in those particular terms is the same as the meanings of the time-delay τ in model (1.5). This paper aims to compare the dynamical behavior of the delayed fractional tumor model with Caputo sense and discretized conformable fractional-order tumor model and it is structured as follows. In Section 2, we perform a stability and bifurcation analysis of the tumor model with Caputo fractional derivative for both non-delay and delay systems. Next, in Section 3, a discretization process based on using piecewise constant arguments is applied to the conformable fractional-order form of the model which leads to a two-dimensional discrete dynamical system. We also deal with stability and bifurcation analysis of the discrete model. Our biologically interpreted numerical findings are reported in Section 4, which is followed by some conclusion in Section 5.

Tumor model with fractional derivative in the Caputo sense
In this section, we will analyze model (1.4). The equilibrium points of model (1.4) are found by solving the equations From D α x(t) = 0, we get x = σ δ+µy− ωy η+y as long as δ + µy = ωy η+y . From the second equation D α y(t) = 0, we get y = 0 or x = γ(1 − βy). By defining the functions f (y) = σ δ+µy− ωy η+y and g(y) = γ(1 − βy), we can classify the equilibrium points of system (1.4) as follows: i) The tumor-free equilibrium point E 0 = (σ/δ, 0), ii) The tumor-infection equilibrium points E = (x, y) where x = g(y), y = y and y is the non-negative solution of the equation By considering equation (2.1), the classification of equilibrium points of the system depending on the parameters can be found in [40].

Stability analysis without delay
For the stability analysis of system (1.4), we use the following theorem.
Theorem 2.1. [4,46] Consider the n-dimensional system of Caputo sense fractional differential equations Let X * be an equilibrium point of the system and J(X * ) be the Jacobian matrix about the equilibrium point X * , then, the equilibrium point X * is locally asymptotically stable if and only if all the eigenvalues Proof The Jacobian matrix of system (1.4) evaluated at the equilibrium point E 0 is given by The eigenvalues of J(E 0 ) are λ 1 = −δ and λ 2 = γ − σ δ . Under the condition R 0 < 1, both eigenvalues are negative real numbers, consequently, E 0 is locally asymptotically stable.
Assume that E = (x, y) is a positive equilibrium point of system (1.4). Linearizing the system about the equilibrium point E = (x, y) gives the following Jacobian matrix and the corresponding characteristic equation is By Theorem 2.1, we deduce the following conclusions: Theorem 2.3. The equilibrium point E = (x, y) of system (1.4) is locally asymptotically stable if one of the following conditions holds.
Although damped oscillations are seen in Figure 1 with respect to the parameter β for the system without delay, these oscillations do not cause Hopf bifurcation or chaotic behavior. We also notice that smaller fractional derivatives damp the oscillation behavior for both effector and tumor cells.

Stability analysis of the delayed model
The equilibrium points of system (1.5) are the same as in system (1.4), and the following theorem explains the stability conditions of the delayed fractional differential equations.
The characteristic equation of system (2.5) is given as det|s α I − A − Be −sτ | = 0. If all the roots of the characteristic equation have negative real parts, then the zero solution of system (2.5) is locally asymptotically stable.
For the numerical simulations of models (1.4) and (1.5), we use the predictor-corrector method introduced in [19] and its generalization to the delayed fractional differential equations [10,56]. Dynamical behaviors of system (1.5) for different values of time delays are given in Figure 2. The system undergoes a Hopf bifurcation at τ ≈ 0.298. Figure 3 displays the dynamical behavior of system (1.5) for different values of the parameter β. For smaller values of β (i.e., for bigger carrying capacity), system (1.5) loses its stability and we observe oscillatory behavior. In Figure 4, we observe stable behavior for smaller fractional derivative. For τ = 0.25 and α = 0.966004, system (1.5) undergoes a Hopf bifurcation and loses its stability.
We rewrite this equation as following: Solving this first-order linear differential equation with respect to t ∈ [nh, t) gives For t → (n + 1)h, we get the difference equation Lastly, substituting E(nh) and T (nh) by E(n) and T (n) gives . (3.1) Applying same procedure to the second equation of system (1.7) we obtain .

(3.2)
By gathering the difference equations (3.1) and (3.2) together, we have the following discrete dynamical system: .

Stability and bifurcation analysis
In this section, we explore the dynamical behaviors of system (3.3). System (3.3) and system (1.5) have the same equilibrium points. The Jacobian matrix calculated at any equilibrium point (E, T ) is given by where c = −δ − T µ + T ω T +η . The corresponding characteristic equation is ) and To investigate the asymptotic stability of the equilibrium point (E, T ), we use the Jury conditions [45]: Since the expression of the characteristic equation and Jury conditions are too complicated, it is difficult to obtain explicit parametric criteria for the stability of the equilibrium point (E, T ). Thus, we will represent the results using numerical simulations after making the bifurcation analysis of system (1.7). Proof. Firstly, we solve the equation 1 − a 0 = 0 (3.4) to find the eigenvalue assignment condition of Neimark-Sacker (NS) bifurcation. We choose β as the critical bifurcation parameter and we get critical β value from (3.4) as follows: (3.5) For β = β, the eigenvalues of J (E,T ) are and a 0 (β) = 1.
Then, we calculate the norm of the eigenvalues: So, the non-strong resonance is fulfilled. Furthermore, the inequality guarantees the transversality condition thus proving the presence of NS bifurcation.
From this point, we make the calculations numerically after giving the formulas. By using the parameter values (except β) given in Table 1 and by taking the discretization parameter h = 0.25 and the fractional-order parameter α = 0.95, we calculate the critical NS bifurcation value β = 0.00157061. Next, we will apply the projection method [41] to investigate the NS bifurcation around the equilibrium (E, T ) = (1.61493, 8.20002).
Then, we can compute the coefficient k(β) by

Results and discussion
In this work, we propose a modification of tumor-immune interaction model (1.2) introduced by Kuznetsov [40]. First, we add an extra delay parameter to model (1.2), and then we consider the fractional form of this model with both Caputo and conformable sense, which are given as model (1.5) and model (1.7), respectively. We note that we use piecewise constant arguments in model (1.7) instead of the delay parameter in model (1.5). After the discretization process based on the PCA method, we obtain system (3.3).
We perform the numerical simulations for x(t) or E(n) (Effector cells), y(t) or T (n) (tumor cells) at parameter values σ = 0.1181, ω = 1.31, η = 20.19, µ = 0.00311, δ = 0.3743, γ = 1.636 and β = 0.002. We note that these  parameters are the dimensionless form of the model parameters whose biological meanings and numerical values are given in Table 1. Figure 1 reveals the effect of the order of Caputo derivative for system (1.4) without delay at β = 0.002 and β = 0.001. It can be noticed from Figure 1 that although periodic solutions with small oscillations emerge according to the changes in two parameters, that is α and β, these periodic solutions do not cause Hopf bifurcation. Both of the effector cells (x(t)) and tumor cells (y(t)) tend to the positive equilibrium point (x, y) ≈ (1.60920, 8.18970) as a result of small oscillations in populations. Figure 2 displays the effect of time-delay parameter τ upon effector and tumor cells with respect to time for system (1.5). It can be seen from Figure 2 that damped oscillations are observed for small values of τ that is tumor dormancy state (Fig. 2a, b). Tumor dormancy is a situation in which the number of tumor cells remains constant at (x, y) ≈ (1.60920, 8.18970) and neither effector cells nor tumor cells are extinct. When τ exceeds  Table 1.
τ 0 = 0.298105, system (1.5) undergoes a Hopf bifurcation as a result of periodic or quasi-periodic solution (Fig. 2c). This state, in which tumor cells exhibit oscillatory behavior without any treatment, known as Jeff's phenomenon and clinically observed [55]. The positive equilibrium point (x, y) remains oscillatory with higher amplitude mode for τ > τ 0 (Fig. 2d). Figure 3 reveals the effect of parameter β on tumor and effector cells for system (1.5) with delay. The parameter β −1 represents the carrying capacity of tumor cells, and its value depends on nutrients provided by the host together with other growth factors. This β −1 value reflects the biological fact that tumor cells population cannot increase infinitely due to limited resources and a lack of space [17]. For smaller values of β −1 , both effector and tumor cells display damped oscillatory behavior (Fig. 3a, b) and finally reach the stable equilibrium point (x, y) ≈ (1.60920, 8.18970). When β −1 is bigger, effector cells remain insufficient to control tumor cells, and both effector and tumor cells exhibit oscillatory behavior (Fig. 3c, d). Figure 4 shows the effect of fractional-order for system (1.5) with delay for τ = 0.25. We obtain limit cycles in phase diagrams for bigger α values (Fig. 4c, d). For smaller α, the system tends to the coexistence equilibrium point (x, y) ≈ (1.60920, 8.18970) (Fig. 4a, b). This is an expected result since smaller fractional derivatives enlarge the regions of stability of a system.
In Section 3, we discretize system (1.7) using properties of the conformable fractional derivative and piecewise constant approximation method. While constructing system (1.7), we replace the terms involving τ of system (1.5) with piecewise constant arguments [ t h ]h. We also show that the discretization parameter h is actually behaving the same way with τ . Piecewise-constant approximation method is about solving the equation in the interval t ∈ [nh, t) with t → (n + 1)h and then generalizing the solution for t > 0. Hence, taking the discretization parameter h bigger, enlarging the interval [nh, (n + 1)h], and consequently uncovering and increasing the delay effect. After the discretization of system (1.7), we get a two-dimensional discrete system (3.3). We note that we  Table 1. take t (days) as nh for the numerical simulations of system (3.3). As in model (1.5), the parameters β, h, and α play a key role in the dynamical structure of discrete dynamical system (3.3). Figure 5 displays the dynamical behavior of discrete-time conformable fractional-order system (3.3) with varying the parameter β, and Figure 6 shows the bifurcation diagram of the system with changing the parameter β. When the parameter β reaches to β = β = 0.00157061, the system enters Neimark-Sacker bifurcation around the positive equilibrium point (E, T ) ≈ (1.60920, 8.18970) (Figs. 5c and 6). For β < β, the positive equilibrium of the system is asymptotically stable with damped oscillations that is tumor dormancy (Fig. 5a, b). However, if β exceeds β, chaotic behavior occurs for the tumor population that is uncontrolled tumor growth (Fig. 5d). Figure 7 shows maximum Lyapunov exponents varying the parameter β corresponding to Figure 6. This figure demonstrates the existence of the chaotic regions and period orbits in the parametric space. Figure 8 and Figure 9 reveal the effect of the order of conformable derivative and discretization parameter h for system (3.3), respectively. It can be seen from Figures 8 and 9 that decreasing the fractional-order parameter α and increasing the discretization parameter h lead to a chaotic attractor for system (3.3). Figure 10 displays a comparison of dynamical behaviors of delayed fractional-order model (1.5) and discrete conformable fractional-order model (3.3) for tumor cells. It can be noticed from Figure 10 that the discrete system exhibits more complex dynamical behavior than the Caputo fractional-order form such as chaotic dynamics.

Conclusions
The aim of this paper is to compare the dynamical behaviors of fractional-order delayed model (1.5) and discrete conformable fractional-order tumor model (3.3). Theoretical and numerical results show that both models have same equilibrium points and exhibit similar dynamical behaviors for the tumor populations such as stable steady state, quasiperiodic solutions with damped and higher amplitude oscillations, and Hopf (Neimark-Sacker) bifurcation (Fig. 10a-c). These dynamic structures emerge due to the change in parameters β, τ (h), α which play a key role in the dynamical behaviors of both models. Hence, it can be concluded that both of  Table 1.
the models of fractional-order are a manageable, easy to employ, and powerful computational mechanism for analyzing the tumor-immune system interaction. We can also say that fractional-order differential equations provide some advantages in modeling tumor-immune interactions. In fractional-order systems, the solution continuously depends on all the previous states. So, using the fractional-order can help to reduce the errors arising from the neglected parameters. In addition, we also observe that the discretization parameter h is actually behaving the same way with delay parameter τ. While Hopf bifurcation occurs in the Caputo sense model when the delay parameter reaches τ = 0.298105, Neimark-Sacker bifurcation exists at the critical value h = 0.267877 in the discrete model (Fig. 10c). However, the discrete system exhibits more complex dynamical behavior than the Caputo fractional-order model such as chaotic dynamics (Fig. 10d). This is an expected result of working on difference equations. Although the number of tumor cells in two graphs overlap better for the initial subinterval of t, the graphs begin to differ as t increases because the amplitude of the periodic solutions is larger and the period of cycles is shorter for the discrete model.
Compared to the Caputo fractional-order model, the discrete conformable fractional-order model has certain advantages. With their simpler form, they allow a relatively simple analysis of stability, bifurcations, and it is easy to calculate population density by an iteration without the need for any numerical methods. Since difference equations exhibit richer dynamic behaviors than continuous models, these equations reflect the chaotic behaviors that already exist in tumor dynamics. Although the conformable derivative is not a fractional-order derivative, it shows a local fractional-order effect on sub-intervals in our dynamical systems. This is not a result to be ignored while considering the difficulty of working with the Caputo fractional-order model. Finally, we can say that the results of the present research work are very important and useful in understanding tumor-immune system interaction.