Fractional calculus in mathematical oncology

Even though, nowadays, cancer is one of the leading causes of death, too little is known about the behavior of this disease due to its unpredictability from one patient to another. Classical mathematical models of tumor growth have shaped our understanding of cancer and have broad practical implications for treatment scheduling and dosage. However, improvements are still necessary on these models. The primary objective of the present research is to prove the efficiency of fractional order calculus in mathematical oncology, more specifically in tumor growth modeling. For this, a generalization of the four most used differential equation models in tumor volume measurements fitting is realized, using the corresponding fractional order equivalent. Are established the fractional order Exponential, Logistic, Gompertz, General Bertalanffy-Pütter and Classical Bertalanffy-Pütter models for a treated and untreated dataset. The obtained results are compared by Mean Squared Error (MSE) with the integer order correspondent of each model. The results prove the superiority of the fractional order models. The MSE of fractional order models are reduced at least at half in comparison with the MSE of the integer order equivalent. It is demonstrated in this way that fractional order deterministic models can offer a good starting point in finding a proper mathematical model for tumor evolution prediction. Fractional calculus is a suitable method in this case due to its memory property, aspect that particularly characterizes biological processes.

Cancer is one of the world's deadliest diseases. According to the World Health Organization, there were an estimated 19.3 million new cases of cancer and almost 9.9 million deaths from cancer worldwide in 2020 1 . A lot of research is being done to find proper models and to apply modern control algorithms to discover personalized treatment for patients to, at least, postpone cancer's evolution with minimal impact of treatment side effects (chemotherapy, radiation therapy etc.).
Mathematical modeling is a powerful tool used in finding different structures that can predict the behavior of a system, even when referring to a living organism 2,3 . Using tumor mathematical models, it can be useful in cancer proliferation and survival signaling, tumor immunology, tumor microenvironment, metastasis and anticancer therapeutic research. Exploring models of various cancer types will help understanding cancer and the corresponding therapeutic approach. It can lead to personalized and efficient treatments.
There are several models developed for different type of tumor evolutions and therapy efficiency. In 4 are reviewed several emerging therapeutic strategies and discussed how mathematical models have contributed to the design of such schedules. All researcher in the field agree that mathematical models can be used to describe and forecast the behavior of cancer. This is one of the main objectives of "mathematical oncology" [5][6][7] . The origins of mathematical oncology can be considered the Gompertz population growth model 8 . Another milestone is represented by the Bertalanffy's organism growth model 9 . More recently several models are developed and discussed. An extensive review of recent results proves the ability of such models to simulate and predict the spatiotemporal development of tumors 10 .
As a starting point for the present research, in 11 can be found the fitting of differential equation models to tumor volume measurements of patients undergoing chemotherapy or cancer immunotherapy for solid tumors. There are compared six classical models which are widely used in the field: the Exponential, Logistic, Classic Bertalanffy, General Bertalanffy, Classic Gompertz and General Gompertz model. The study concludes that these models could potentially be effective at predicting treatment outcome.
On the other hand, fractional order calculus becomes more and more used in biological, chemical and medical systems 12 . It generalizes the classical, integer order differential calculus to non-integer orders. Many researchers have already proved the superior performance of fractional order models for describing this phenomenon. Great results are presented in different areas: from biosciences, bioengineering, medicine, economics to physics, www.nature.com/scientificreports/ control, signal processing, neural networks [13][14][15][16] . Many authors agree that fractional operators are useful in capturing and understanding the more relative consequences of physical phenomena with higher nonlinearity and complexity with long-range memory and history-based properties 17 . For example paper 18 analyzes the dynamics of a fractional partial differential equation model of Zika virus. Incorporating the diffusion phenomena using Atangana-Baleanu fractional derivative, the authors explain how the spread of humans and mosquitoes influences the disease's transmission. On the other hand mathematical oncology is both an old and a new field of research. Recently, multidisciplinary approaches became more and more attractive. Predictive mathematical modeling is widely used to fill in gaps between mathematicians, clinicians, biologist and to make predictions of cancer progression and response to therapy on a patient-specific basis. Many papers are published in the field, but all researchers agree that mathematical oncology is still in the early stage. Any contribution to the field is welcomed.
The present work try to combine this two field by extending the classical tumor growth models using the fractional differential operator. The results are compared by analyzing the mean squared errors (MSE) with both real data and the integer order models, highlighting the advantages introduced by fractional calculus. The obtained results show the utility and ability of fractional order models to emphasize certain characteristics that the integer order systems cannot describe.
The paper is structured as follows. After this introductory part, Section "Materials and methods" presents the used methods, while Section "Results and discussions" discusses the obtained results. The work ends with a concluding section.

Materials and methods
The concept of integral operator is defined in several ways 19 . The most used is the Riemann-Liouville definition, which states that a fractional order integral of order ℜ(α) > 0 is a natural consequence of Cauchy's formula for repeated integrals, expressed as 19 : where I represents the notation for the integral operator and n is a natural number, the order of integration. As it can be seen, there is a constraint regarding the order n of this operator, due to the term n!. Introducing the Gamma function in the above formula, the notion of integration order can be extended from the set of natural numbers to the set of positive real numbers, so the new form of the integral becomes: where α ∈ R + is the new order of integration and is the Euler's Gamma function which is a generalization of a factorial. The generalization for the whole set of real numbers becomes: which describes the real order integral operator, keeping at the same time the inverse correspondence between differentiation and integration.
This definition became 19 : The fractional-order derivative of order α ∈ R + can be defined using the Riemann-Liouville formula 19 : or by the alternative definition introduced by Caputo 16 : (1) www.nature.com/scientificreports/ The advantage of Caputo fractional derivative 20 is allowing the initial and boundary conditions to be included in the problem. Moreover, this form of the integral overcomes the Riemann-Liouville integral drawback of having the derivative of a constant different from zero.
A useful aspect of the variation of the integration order can be represented intuitively in Fig. 1 where a simple function of the form f (x) = x is considered. It is obvious that the first order derivative is 1, while the derivative of 0 is the function itself. Varying the integration order between 0 and 1 we can obtain a function that oscillates between the two curves as shown in Fig. 1.
Fractional operators are widely used to describe complex behaviors of dynamical systems 21,22 and proved to be a powerful tool in describing and understanding the complex behavior of nonlinear systems.
The mathematical models established and discussed in present work are deterministic models. A deterministic model allows to calculate a future event exactly, without the involvement of randomness. This class of models can predict the outcome with certainty which from the beginning may be a wrong assumption regarding the complex behavior of cancer but can provide a good starting point for finding new models and to extend the research to stochastic models. For all models, the dependent variable is the volume of the tumor as a function of time.
The first model considered is the exponential model. This type of model has the simplest form and is based on the following differential equation 6 : where v is the tumor volume, a represents the growth exponent, the kinetic parameter.
The second analyzed model is the logistic model, with general form 6 : where a is the inherent growth rate, k is the average population size of a species in a habitat (in this case, the volume of cancerous cells in a living organism) and b is an exponent that corrects the expression of tumor growth rate. The third model is the Gompertz model 23 : which is practically a generalization of the logistic model. The constant c introduces the volume minimum carrying capacity. Its advantages are proved mainly in the field of biology, helping to describe the evolution of animals and plants, as well as the number or volume of bacteria or cancer cells in a living organism. The fourth model is the Bertalanffy-Pütter model, having the form 24 : where the first definition represents the generalized form of the Bertalanffy-Pütter model, the second standing for the particular form when the exponents a and b are equal. The parameters in this case are the pair of exponents (a, b) and coefficients (p, q) where p is the intrinsic growth (meaning the number of new cells minus the number of dead cells during a generation) and q is the growth rate declaration factor of antiangiogenic process (referring to the vascularization process). However, these models do not include the effects of memory, which are found in biological systems. Hence, in order to take into account the memory effects by the mathematical formulation, it is introduced a new form for each model by replacing the ordinary derivative with the fractional operator. For each model the fractional order differential correspondent is proposed. The existence and unicity of the above-mentioned fractional order models of tumor growth are proved, starting from the following form of Riemann-Liouville integral: , the following definition has been proposed, by changing the kernel from , thus obtaining the Caputo-Fabrizio definition of the integral operator 25,26 : where M(α) is a constant which depends on α.
Using Losada-Nieto proposition and considering D α f (t) = h(t), it is obtained: In the upcoming demonstration fixed-point theory and Picard-Lindelof technique will be used 27 .
In terms of kernel the following equation holds: where K(v,t) is one of the models mentioned above (e.g. K(v, t) = av(t)).

Proof
A. Existence of the solution: Let y(t) = v(t) for the considered model F y, t = y , with y(0) the model's initial condition. Now the Picard iterations are: Repeating the algorithm above we get: For simplicity, we will consider the model (11).
, is a fractional differential equation with a fractional order x.
To prove the stability of the solution v(t) to this fractional differential equation, we can use the concept of Mittag-Leffler stability, which is a generalization of Lyapunov stability for fractional differential equations.
Specifically, we can show that the solutions of the fractional differential equation converge to a steady-state solution as time approaches infinity if the Mittag-Leffler function satisfies certain conditions. The Mittag-Leffler function E α (z) is defined as: where α is a constant exponent that determines the rate of decay or growth of the function for large values of z. The general solution to the fractional differential equation d α v dt = a · v(t) , is given by: where A is a constant of integration.
To show the Mittag-Leffler stability of this solution, we need to demonstrate that E α ((a · t) x ) converges to a finite value as time goes to infinity for certain values of α and a.
For 0 < x < 1 , the Mittag-Leffler function E α (z) is monotonically increasing and bounded for 0 < z < infinity and all α > 0 . Therefore, for 0 < x < 1 and any a, E α ((a · t) x ) converges to a finite value as t goes to infinity. For x ≥ 1 , the Mittag-Leffler function E α (z) is not necessarily bounded, but it has a Laplace transform that is bounded for Re(s) > α , which means that E α ((a · t) x ) decays exponentially for certain values of α and a. Specifically, E α ((a · t) x ) is bounded for 0 < α < x and a ≤ 0 , which implies that the solution v(t) = A · E α ((a · t) x ) is Mittag-Leffler stable for these values of α and a.
Therefore, we can conclude that the solution v(t) = A · E α ((a · t) x ) to the fractional differential equation is Mittag-Leffler stable for certain values of α and a, depending on the order x of the fractional derivative.
Having these proofs, all the above models can be taken into account as being valid to describe the biological process considered. All parameter used in the models are supposed to be positive constants, being a biological system. All parameters have the same definitions as in the integer order case, appearing the additional degree of freedom, α, the fractional order.
The used database is from research 31,32 , where tumor pieces from a mouse tumor were transplanted orthotopically to syngeneic FVB mice. The tumor's volume was measured over a period of 18 days. One group was observed without treatment and a second one was treated based on clinician's protocol. The treatment consisted in drug injection when the tumor volume reached 200 [mm 3 ] and at least 10 days have passed since the last treatment. The injected drug dose used was the maximum tolerable dose of 8 [mg/kg]. If the maximum volume threshold  3 ] was reached, the experiment was stopped. The observations were made regularly, once every two to three days, the system's identification being performed on the interpolated data set using spline curves. There were used two sets of measurements to study the tumor's growth: the first case on which a treatment scheme was applied and the second case where only observations were made on the evolution of the tumor's volume. For all enumerated models the dependent variable is considered the volume of the tumor as a function of time. The goal is to fit each type of models to the entire time series for each dataset. The statistical endpoint of each experiment is the minimum Mean Square Error (MSE): where n represents the number of samples for each model, y i are the real data values and y i the values estimated by the models. For simplicity as initial conditions were used near zero values.
All model fitting procedures were implemented in Matlab ® . For the fractional order calculus FOMCON toolbox is used 33 . The optimal parameters for the integer order models were obtained using the fminsearch function from Matlab ® . The fractional order models were used with the same parameters as already fitted for the integer order ones, establishing the best fit by only changing the fractional order. In this way it can be highlighted the effect of the introduced fractional order on the tumor model dynamics.

Results and discussions
The simulation results for two dataset for both treated and untreated case are plotted in Fig. 2. These graphics reveal that the fractional order derivative has a significant impact on the dynamics of the tumor evolution.
For each considered fractional order, MSE is computed for the resulting models. The impact of the fractional order on the MSE is plotted in Fig. 3. It can be observed that there is an optimum value in each case.
The resulted best fitting integer order models and the corresponding fractional order equivalent, having the same coefficients, are presented in Table 1, while Table 2 highlights the obtained mean squared errors: The evolutions of the MSE for each considered model for both integer and fractional order version for the treated and untreated case data are plotted in Fig. 4. With blue are plotted the integer order model results and with read the fractional order correspondents.
All obtained results for the integer case are in accordance with the results published in the field 34 . For the fractional order models no results were found in the literature. From the MSE evolution plots for different fractional order values can be concluded that in each case there is an optimum value. Decreasing the fractional order leads to the decrease of MSE until this minimum MSE is reached.
It can be noted that the MSE of fractional order models are reduced at least at half in comparison with the MSE of the integer order equivalent. This good model fitting with fractional order models is even more obvious for the untreated case, where even a 26.21% of MSE reduction can be obtained.
The best integer order model for the treated case is the Gompertz model. However, with the fractional order equivalent a MSE reduction from 0.0218 to 0.0146 is obtained. It is interesting to see, that a MSE not very far from this value (0.021) is obtained with the more simple fractional order logistic model. The extra degree of freedom offered by the fractional order makes a useful tool even from a relatively simple model, like the logistic models.
The best integer order model for the data obtained from untreated tumors is the generalized Bertalanffy-Pütter model, with a MSE of 0.8557. The fractional order equivalent leads to a MSE of 0.2243, which means a 26.21% MSE reduction. As in the untreated case, it can be observed that even the simplest model, the exponential model, can lead to a small MSE if fractional order is used.
After all these tests, performed to answer the question how well classical differential equation models can fit tumor volume trajectories both in treated and untreated case, it is found that in all cases fractional order versions of each model outperforms the corresponding integer order model. Fractional order models are more flexible in fitting empirical data, captures features of compliance data and offer improved model predictions, proved by the reduced MSE.

Conclusions
Classical mathematical models are in principle useful to model cancer growth for both treated and untreated case. Moreover, although several studies analyze such models on different data, to the best of our knowledge, no fractional order models are developed for tumor models. Deterministic structures can offer a good starting point in finding a proper mathematical model for tumor evolution prediction, but they can be improved by using fractional differential calculus in order to improve the approximation obtained from the integer order differential equation. Fractional calculus is a suitable method in this case due to its memory property, aspect that particularly characterizes biological processes.
The present study discusses the fractional order generalization of four generally recommended mathematical models: exponential, logistic, Gompertz and the generalized and particular form of the Bertalanffy-Pütter model.
As presented in the above results, the best results for the treated tumor case were obtained using the fractional order logistic model, which offers MSE similar to the integer order Gompertz model, recommended in all published studies. This model offers good results for the untreated tumor case as well. For this case the MSE is 0.3002. However, the smallest MSE in this case is offered by the generalized Bertalanffy-Pütter model, having a MSE of 0.2243. It is explicable, being the more complex model, with the most parameters.
Although those models offer a good understanding of the tumor's growth process and cancer development, they have a major disadvantage by not being able to provide the necessary structure of a model (unlike a statespace model) so that it can be used in developing a control strategy. In addition to this, those models do not www.nature.com/scientificreports/ www.nature.com/scientificreports/ include the presence of an input (explicit treatment used on mice), disturbances (changes in the metabolism of the studied specimen or the appearance of certain diseases or abnormalities) or even a resistance to the applied treatment.
To sum up, further research is to be done regarding the use of fractional calculus in control theory, but the preliminary results show their utility and ability to emphasize certain characteristics that the integer order systems cannot observe. Also, stochastic models are proven to be a powerful tool in finding better predictions for processes that are too complex to be described in a deterministic way. Future research include a mix between stochastic and fractional calculus, which can be a solution in finding a proper model to describe any biological process inside a living organism, including cancer. www.nature.com/scientificreports/ www.nature.com/scientificreports/