MEMORY EFFECTS ON THE PROLIFERATIVE FUNCTION IN THE CYCLE-SPECIFIC OF CHEMOTHERAPY

. A generalized mathematical model of the breast and ovarian cancer is developed by considering the fractional diﬀerential equations with Caputo time-fractional derivatives. The use of the fractional model shows that the time-evolution of the proliferating cell mass, the quiescent cell mass, and the proliferative function are signiﬁcantly inﬂuenced by their history. Even if the classical model, based on the derivative of integer order has been studied in many papers, its analytical solutions are presented in order to make the comparison between the classical model and the fractional model. Using the ﬁnite diﬀerence method, numerical schemes to the Caputo derivative operator and Riemann-Liouville fractional integral operator are obtained. Numerical solutions to the fractional diﬀerential equations of the generalized mathematical model are determined for the chemotherapy scheme based on the function of “on-oﬀ” type. Numerical results, obtained with the Mathcad software, are discussed and presented in graphical illustrations. The presence of the fractional order of the time-derivative as a parameter of solutions gives important information regarding the proliferative function, therefore, could give the possible rules for more eﬃcient chemotherapy.


Introduction
Chemotherapeutic drugs are cycle-specific, namely, they destroy the cells in specific phases of their cycles.Some drugs work in the gap period (phase G1) and in the synthetic period (the replication phase, S).Other drugs are acting in the mitosis phase (the mitotic phase, M) or in the second gap period (G2 phase).Webb [40] has investigated mathematical linear and nonlinear models of cycle-specific chemotherapy.For linear models, the author has studied the advantages of periods of dose with shorter duration.A mathematical model with a constraint equation that describes the effects of the drugs on sensitive normal tissue, based on the pulsed and piecewise continuous chemotherapeutic effects was studied by Panetta and Adam [30].They have identified the optimal period needed for maximal reduction of the tumor.
A mathematical model that takes into account by the cells sensitive to the treatment (the proliferating cells), and the quiescent cells (the phase G0), which are resistant to the treatment has been developed by Panetta [31].The model parameters have been estimated for the breast, ovarian cancers as well as for the bone marrow.
The estimated parameters suggest acceptable treatment strategies regarding the treatment period, time of drug infusion, etc. Liu et al. [24] studied the effects of an M-phase specific drug on the development of cancer by including the resting phase G0 and the immune response.In their model, the authors included a time delay for the passage through the inter-phase and assumed that the immune cells interact with all cancer cells.The authors have found that the dynamics of the G0 phase dictates the dynamics of cancer as a whole.
Adam and Panetta [1] studied the model for cell-specific chemotherapy by reformulating the differential equations of the model as a Schrodinger equation in time.The effects of chemotherapy on cycling tumor cells have been considered as an exponentially decaying function while the potential function is a Morse-type potential.The solutions are expressed in terms of Whittaker functions.
An interesting review in the area of mathematical modeling in biology can be found in references [19,22,34].The readers find mathematical and computational models and methods in the fields of cell physiology, epidemiology, and cancer.A new mathematical model that combines interactions between tumor cells and cells in the immune systems combined with drug delivery has been developed by Unni and Seshaiyer [37].The system of ordinary differential equations that describe the model has been numerically solved.Their results seem to suggest that the model used to study the dynamics of tumor cells; it helps to provide the dynamic interactions between the tumor cells, immune system, and drug response systems.
Weerasinghe et al. [41] have presented a set of mathematical models of cancer cell plasticity, namely the process in which, through genetic and epigenetic changes, cancer cells survive in hostile environments and migrate to more favorable environments, respectively, tumor growth, and invasion.The authors discussed models of plasticity, tumor progression, and metastasis using discrete, continuum, and hybrid mathematical models.Usher [38] has developed a mathematical model to investigate the effect of toxic drugs on tumor growth in order to obtain more efficient chemotherapies.
Dingli et al. [10] have formulated a mathematical model to describing the radiovirotherapy as a combination of virotherapy with radiation, used to eliminate tumors when virotherapy alone fails.This model is based on population dynamics that captures the essential elements of radiovirotherapy.Authors have investigated the existence of equilibrium points related to partial/complete cure, and therapy failure.Mathematical models described by ordinary differential equations, algebraic equations, and partial differential equations used to the characterizing tumor burden dynamics are presented by Yin et al. [44].Also, the authors discussed stochastic and deterministic models for the tumor resistance evolution and emphasized the possibility to develop a novel model that incorporates both tumor dynamics and resistance evolution.
Liu and Yang [25] have studied a periodic mathematical model of cancer treatment by radiotherapy.They have obtained conditions on the existence and globally asymptotic stability of the positive periodic solution and the cancer eradication periodic solution.
Wang and Schattler [39] have formulated cancer chemotherapy as a problem of optimal control in order to find the conditions in which are minimizing the tumor volume and side effects over a specified therapy horizon.A mathematical model to describe the immune response to a vascular tumor under the influence of immunotherapy and chemotherapy and their combinations as well as vaccine treatments has been formulated and studied by Isaeva and Osipov [17].The effect of vaccine therapy has been considered as a perturbative parameter of the model.By numerical simulations, the authors found that the efficiency of vaccine therapy depends on both the tumor size and the condition of the immune system as well as on the response of the organism to vaccination.For a strong immune response, their model predicts the tumor remission under vaccine therapy.
Recently, various fractional-order operators (Riemann-Liouville, Caputo, Caputo-Fabrizio, Hilfer, Prabhakar, Riesz, Atangana-Baleanu, etc.) are used in mathematical modeling of several complex transport processes [5, 14-16, 26, 33].Since most biological systems have memory or after-effects, such as the delay due to the incubation time for vectors to become infectious, the modeling of biological systems by fractional differential equations is more advantageous than the classical modeling, in which the memory effects are neglected.
Ahmed et al. [2] investigated a fractional cancer model with Caputo time-fractional derivative by taking into consideration the cross-reactivity of the immune system.A new fractional-order mathematical model for a tumor-immune surveillance mechanism, based on the fractional derivative with Mittag-Leffler kernel has been studied by Baleanu et al. [4].The authors have developed an efficient numerical procedure to solve the fractional differential equations of the mathematical model.Also, they suggested an optimal control strategy for investigating the effect of chemotherapy treatment.
Iyiola and Zaman [18] have considered a cancer tumor model with Caputo fractional derivative as compared to the model with integer order of time-derivative, for the different net killing rate of the cancer cells.Analytical solutions are determined using the q-homotopy method.Mathematical models of breast cancer that consider population dynamics among the cancer stem cells, tumor cells, healthy cells, the effects of excess estrogen, and the body's natural immune response on the cell populations have been studied by Solis-Perez et al. [32] In the studied models, the authors have used time-fractional derivatives with power-law and exponential law in the Liouville-Caputo sense.
Manimaran et al. [27], by adapting the Faedo-Galerkin method, have proved the uniqueness and the existence of weak solutions of the time-fractional cancer invasion system with non-local diffusion operator.Numerical solutions have been obtained using a finite element numerical scheme.Baba [3] studied a fractional model that considers the effect of intravesical bacillus Calmette-Guerin (BCG) in controlling bladder cancer.Equilibrium points are determined; the existence, uniqueness and stability of solutions of the system have been analyzed.A mathematical model of the type SEIR (susceptible, exposed, infected, and recovered) for the dynamics of infectious diseases has been studied by Khan [20].
The mathematical model of interaction cancer cells and immune system cells has been numerically investigated by Ucar et al. [35] using the Adam-Bashforth-Moulton algorithm.Since the cancer cells have the memory structure, the authors have used fractional models that are closer to the reality and give more meaningful results than the classical models.Using Caputo time-fractional derivative, the authors analyzed the behaviors of system cells by changing of the fractional parameter.Ucar et al. [36] have studied the dynamics of smoking mathematical model with Atangana-Baleanu derivative.They proved the existence and uniqueness of solutions via fixed point theory.A SIR (susceptible, infected and recovered) epidemic mathematical model with exponential decay law has been investigated by Yavuz and Ozdemir [42].The studied model is based by the time-fractional Caputo-Fabrizio derivative.The existence and uniqueness of solutions are demonstrated; the influence of the memory parameter on the infectious diseases evolution also was investigated.Yavuz and Sene [43] developed a fractional predator-prey model with the harvesting rate using the time-fractional Caputo derivative.The existence and uniqueness of the solution, local stability and global stability have been investigated.The authors have introduced a novel discretization depending on the numerical discretization of the Riemann-Liouville integral and the corresponding numerical discretization of the predator-prey fractional model has been obtained.Using the stability criteria, the authors have examined the dynamical behavior of the equilibria.Naik et al. [28] developed a dynamical fractional order HIV-1 model in Caputo sense which involves the interactions between cancer cells, healthy CD4+T lymphocytes, and virus infected CD4+T lymphocytes that leads to the chaotic behavior of the infected cells, while, Evirgen et al. [12] have studied an HIV infection model of CD4+Tcells with the Atangana-Baleanu fractional derivative.The existence and uniqueness of the solutions for fractionalized HIV disease models with the fractional have been investigated.Ozdemir et al. [29] presented a fractional-order SEIR-KS model Caputo time-fractional derivative to analyze the effects of kill signal nodes on the virus propagation.The existence and uniqueness of the model's solutions are studied.Using the Adams-Bashforth-Moulton algorithm, the authors have obtained numerical solutions using the MATLAB program.A dengue epidemic model with Caputo fractional derivative to analyze the effect of temperature on the spread of the vector-host transmitted dengue disease has been formulated by Defterli [8].The stability of the equilibrium points of the considered dengue model was studied and the corresponding basic reproduction number has been derived.By numerical simulations, the influence of the temperature on the dynamics of the vector-host interaction in dengue epidemics was analyzed.Koka [21] developed a mathematical model to describing the rubella spread using the non-local and non-singular Atangana-Baleanu time-fractional derivative for the inclusion of the memory effects.Ertas [11] formulated a particle swarm optimization algorithm and its comparison with the Levenberg-Marquardt and the trust-region optimization algorithms in least squares fitting of the intravoxel incoherent motion model to the breast diffusion MR signals.
The aim of this paper is to develop a generalized fractional model of the mathematical model studied by Panetta [31].Our model describes the effects of cell-cycle-specific drugs on cancer and normal tissues by taking into consideration memory effects.To develop the fractional mathematical model, the time-fractional Caputo derivative is employed.
Even if the classical model, based on the derivative of integer order has been studied in many papers, its analytical solution is presented in order to make the comparison between the classical model and the fractional model.
Using the finite difference method, numerical schemes to Caputo derivative operator and Riemann-Liouville fractional integral operator are obtained.Numerical solutions for the fractional differential equations of the generalized mathematical model are determined for the chemotherapy scheme based on the function of "onoff" type.Numerical results obtained with the Mathcad software are discussed and presented in graphical illustrations.The obtained solutions for the generalized mathematical model studied in the present paper are new in literature.Also, another novelty of this study is the comparison between the proliferative function corresponding to the derivative of integer order and the proliferative function corresponding to the fractional derivative of Caputo type.This comarison highlights that the time-evolution of the proliferating cell mass, the quiescent cell mass, and the proliferative function are significantly influenced by their history.
The presence of the fractional order of the time-derivative as a parameter of solutions gives important information regarding the proliferative function; therefore, the studied model in this paper gives some qualitative ideas on how to better implement cycle-specific chemotherapy.

Mathematical model of cancer with derivative of integer order
In this section, we make a brief presentation of the mathematical model of the chemotherapy effects studied by Panetta [31].We present this model together with the analytical solutions to have the basis of the formulation of model based on differential equations with derivatives of fractional order.Also, the comparison of the results offered by the two models is considered in order to determine the conditions in which the effects of chemotherapy lead to the best treatment strategies.
Basic equations of the mathematical model are [30,31] where P (t) is the proliferating cell mass (cells which are sensitive to the treatment), Q(t) is the quiescent cell mass (cells which are resistant to the treatment), α > 0 is the rate from proliferating to resting, β > 0 is the transition from resting to proliferating, γ > 0 is the growth rate of cycling cells, δ > 0 is the natural cells decay, λ > 0 is the quiescent cell loss, s > 0 is the parameter that describes the strength of the treatment and f (t) is a periodic function that describes the effects of the chemotherapeutic treatment.Equations (2.1) and (2.2) can be written in the equivalent forms where (2.5)

Analytical solution for the "on-off " function
In this study, we consider the function f (t) that describes the effects of the chemotherapeutic treatment of the type "on-off", i.e. the drug is either active or inactive.The form of this function is Function (2.6) can be written in the equivalent form where is the Heaviside unit step function.Note that for any t ∈ [(n − 1)T, nT ], n = 1, 2, 3, . . . the values of function f (t) belong the set 0, T t0 , therefore, the system of equations (2.3) and (2.4) is a system of differential equations with constant coefficients.
The graph of function f (t) is given in Figure 1.

The optimal active phase t 0 of the treatment
To analyze the efficiency of treatment, the values of functions in the initial moment are compared with the values of functions in the final moment of the treatment interval.More exactly, the values P n1 (τ n ), Q n1 (τ n ), n = 1, 2, . . ., are compared with values P n2 (nT ), Q n2 (nT ),n = 1, 2, . . . .Using equation (2.16) for t = nT , we obtain the following condition: We introduce the matrix A straightforward calculation leads to the following expressions of matrices in equation (2.18): (2.19) 2, 3..., r being the proportionality factor.Equation (2.17) becomes From equation (2.21) it results that the proportionality factor rdenotes the eigenvalue of the matrix (2.18).The eigenvalues of the matrix W (t 0 ) are called by Panetta [31], "the characteristic multipliers".
It is observed that the magnitude of cell mass will decay in a period of treatment if and only if the characteristic multipliers are less than 1.Therefore, we are able to determine the value of the active phase t 0 for which the eigenvalues are less than 1, respectively the value for which after the treatment, the cell mass decreases.
The problem of finding the values of t 0 for which the eigenvalues are less than 1 will be solved graphically using the software Matcad 15.
We denote by Ω 1 (t 0 ) = T r(W (t 0 )), Ω 2 (t 0 ) = det(W (t 0 )).The eigenvalues of the matrix W (t 0 ) are given by Obviously, to analyze the evolutions of cell mass functions P (t) and Q(t) during a treatment period T , we will study the maximum characteristic multiplier r 1 (t 0 ).
Figure 2 shows the dependence of the active phase t 0 on the period of treatment T .Along with a period of T = 20 days and a dose of strength s = 0.12 we observe that, for values of t 0 < 7.647, the maximum characteristic multiplier is greater than 1 and thus there is cancer cell growth, whereas, for values of t 0 ∈ (7.67, 28.25), the maximum characteristic multiplier is less then 1 and thus there is cancer cell decay.
Also, let's note that the acceptable values of the active phase are increasing with period T at the same strength s.Therefore, we can directly see that longer active phases cause cell decay, whereas shorter active phases allow cell growth.
The influence of the dose strength s on the active phase t 0 is emphasized by Figure 3 for the period T = 21.As expected, the values of the active phase are decreasing with parameter s.In the analyzed case, for a period of treatment of T = 21 days, we have obtained the following conditions for the active phase: t 0 > 11.196 if s = 0.11, t 0 > 8.346 if s = 0.12, t 0 > 6.688 if s = 0.13.Another parameter which gives important information about the treatment effects is the proliferative function defined as: . (2.23) The profiles of functions P (t), Q(t) and r(t) are presented in Figure 4, for different values of the parameter α-the transition rate from proliferating to resting.
It is observed in Figure 4 that, along with a period of 21 days and a dose of strength of 0.13, the proliferative function r(t) is decreasing with the rate of proliferating to resting parameter α.The proliferative function profiles shown in Figure 4 are clearly highlighted by the second form of the equation (2.23).It can be seen in Figure 4 that the function P (t) decreases with the parameter α, while the function Q(t) increases with the same parameter; therefore, the ratio Q(t)/P (t) is increasing and the proliferative function is decreasing with the parameter α.

Mathematical cancer model with Caputo time-fractional derivative
In this part of the paper, we consider a generalization of the model studied in the previous section.The generalized model is described by fractional differential equations with a time-fractional derivative with a powerlaw kernel.By considering the generalized model, the proliferating cell mass is influenced by his history.
First, we will make a brief presentation of the theoretical notions about the Caputo derivative and fractional integral Riemann-Liouville operator necessary in the presentation of the generalized model.

Some preliminaries
Let us remind some well-known mathematical tools regarding to Caputo fractional derivative and Riemann-Liouville fractional integral needed to use in the present study.

Generalized mathematical model of cancer
In this section, we consider the mathematical model that describes the effects of cell-cycle specific drugs on cancer and normal tissue based on the time-fractional Caputo derivative, describes by the following fractional differential equations: along with the initial conditions, The significations of parameter and functions are the same used in previous sections.

Numerical solutions of the fractional model for
Using the expression (2.6) of the "on-off" function for the drug administration and the notations the system (3.14) and (3.15) is written in the equivalent form Applying to the equations (3.18) and (3.19) the Riemann-Liouville integral fractional operator of order ε, and using the property (3.13), we obtain the system A straightforward calculation, based on equations (3.18), (3.20) and (3.21), leads to the following fractional differential equations: The solutions of equations (3.22) and (3.23) will be obtained numerically using approximations of the finite differences type [6,23].For this, we consider a uniform discretization of the interval [0, t 0 ] with the points t n = nh 1 , n = 0, 1, . . ., N 1 , h 1 = t0 N1 .
For the time interval [t 0 , T ] it is considered the uniform discretization by points N2 .The values of functions P (t) and Q(t)in the point t = t 0 = σ 0 , determined in the previous section, are denoted here by P0 = P (σ 0 ) = P (t 0 Using the same method as in the previous section, we obtain the approximation relations respectively, the values of functions P (t) and Q(t) given by

A graphical study of the particular case ε → 1
As it has been shown in Section 3, when the fractional parameter ε → 1, the time-fractional Caputo derivative tends to the derivative of the integer-order.This means that graphs of the solutions of the generalized model must tend to become identically with the graphs of the ordinary model if the fractional parameter tends to 1. Figure 5 is drawn to highlight this property.As seen in Figure 5, when the fractional parameter approaches 1, the corresponding graph for ε < 1, ε → 1 tends to overlap with that corresponding to the classical case ε = 1.For the curves in Figure 5 we used the following values of the system parameters: α = 0.218, β = 0.08, γ = 0.579, δ = 0.477, λ = 0.001, s = 0.13, T = 21.
4.2.Effect of the memory parameter ε and of the time t 0 of the active drug on the proliferative function r(t) Figures 6-8 stand out the influence of fractional parameter (the memory parameter) on the proliferative function.To these figures we have used the parameters values α = 0.218, β = 0.08, γ = 0.579, δ = 0.477, λ = 0.001, s = 0.2 and T = 21.It is observed from these figures that the fractional models can give cases for which the proliferative function decreases and in the interval in which the drug is inactive (t > t 0 ).Such behavior is not found for the ordinary case corresponding to the ε = 1.In the case of the fractional model, a convenient value of the fractional parameter could be chosen so that at the end of the active period or at the end of the total treatment period, the proliferative function has a specified value.Such a choice might suggest certain forms of treatment.Regarding the active time active t 0 , for the values of the parameters considered in this analysis, a lower value of the active time is more efficient regarding the decrease of the proliferative function.As expected, the strength s of the treatment has significant influence of the values of the function P (t) as seen in the Figure 9.The profiles of the proliferating cell mass P (t), have been plotted in Figure 9 for α = 0.218, β = 0.08, γ = 0.579, δ = 0.477, λ = 0.001, t 0 = 12, T = 21 and for different values of the parameters s and ε.
It is observed that for small values of s, the proliferating cell mass decreases if t < t 0 , and increases for t > t 0 .It is important that such a situation must be avoided.
By increasing values of parameter s, acceptable values of the function P (t) could be obtained.Also, pairs of values of the parameters s and ε could be determined to obtain specified values of the function P (t) in t 0 and T.

Conclusion
In this paper the integer order mathematical model of cancer [31] has been generalized by using Caputo fractional derivative.The fractional mathematical models of cancer are appropriate to the real-life problems because the immune system and the cancer cells have memory features.
A generalized fractional model of time-evolution of proliferating cell mass, quiescent cell mass and proliferating function have been formulated using the Caputo derivative of fractional order.
New analytical solutions of the classical model have been determined.Using the obtained solutions, optimal values of the active drug time have been determined.Also, the solutions of classical model are compared with those of the fractional model.
Numerical schemes to Caputo derivative operator, Riemann-Liouville fractional integral operator, and the fractional differential equations of the generalized mathematical model are determined for the chemotherapy scheme based on the function of "on-off" type.The obtained results are new in the literature.
Numerical results obtained with the Mathcad software are discussed and presented in graphical illustrations.The use of the fractional model shows that the time-evolution of the proliferating cell mass, the quiescent cell mass, and the proliferative function are significantly influenced by their history.
The presence of the fractional order of the non-local time-derivative as a parameter of solutions gives important information regarding the proliferative function; therefore, the studied model in this paper gives some qualitative ideas on how to better implement cycle-specific chemotherapy.

Figure 2 .
Figure 2. The characteristic multiplier r 1 (t 0 )versus active phase t 0 for different values of the period T and for s = 0.12.

Figure 3 .
Figure 3.The characteristic multiplier r 1 (t 0 )versus active phase t 0 for different values of the strength s and for T = 21.

Figure 4 .
Figure 4.The profiles of functions P (t), Q(t) and r(t), versus time, for T = 21, s = 0.13 and for different values of the parameter α.

Figure 6 .
Figure 6.The profiles of proliferative function r(t) for t 0 = 8 and different values of the fractional parameter.

Figure 7 .
Figure 7.The profiles of proliferative function r(t) for t 0 = 12 and different values of the fractional parameter.

Figure 8 .
Figure 8.The profiles of proliferative function r(t) for t 0 = 16 and different values of the fractional parameter.

Figure 9 .
Figure 9.The profiles of the proliferating cell mass P (t) for different values of s and ε.

4. 3 .
Combined effects of the fractional parameter ε and strength of the treatment s on the proliferating cell mass P (t)