A Fractional Modeling of Tumor–Immune System Interaction Related to Lung Cancer with Real Data

In this study, we investigate a new fractional-order mathematical model which considers population dynamics among tumor cells-macrophage cells-active macrophage cells, and host cells involving the Caputo fractional derivative. Firstly, the stability of the positive steady state of the model is studied. Subsequently, the conditions for existence and uniqueness of the solutions are examined. Then, the least squares curve fitting method (LSCFM) which is one of the prominent methods for parameter estimation is used to fit the parameters of the model. It is aimed to fit the relevant parameters with the help of the tumor tissue samples which were collected from the patient with non-small cell lung cancer who had chemotherapy-naive hospitalized at Kayseri Erciyes University hospital in Turkey. A total of 12 parameters in the model are estimated using the data of lung tumor cells of this patient for 14 days. Moreover, the numerical simulations are given by considering the different fractional orders and different parameters for the model. So, it is achieved how the change in α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} affects the dynamic behavior of the system. In the sequel, to point out the advantages of the fractional-order modeling, the memory trace and hereditary traits are taken into consideration. Finally, the interpretations in terms of biological science are provided in conclusion. We believe that this interdisciplinary study will open new doors for other similar studies and will shed light on the studies to be developed on the use of real data in the mathematical modeling of cancer.


Introduction
Cancer is the deadliest and complicated disease of our time. This illness is caused by the uncontrolled growth of abnormal or mutated cells in the body, and cancer cells are known for their capability to grow rapidly, divide and proliferate uncontrollably [1,2]. The most important cause of tumor progression is constant proliferation and metastatic potential [3]. Lung cancer, the most common malignancy in the world, is the number one cause of cancerrelated death worldwide and is responsible for approximately 228,820 new cases and 140,730 deaths in 2020 [4,5]. The 5-year relative survival rate is 6 percent for patients diagnosed with advanced disease [6]. There are many causes why lung cancer is frequently diagnosed at an advanced stage, including inadequate early detection biomarkers, disease following methods and physical examination [7]. Lung cancers are conventionally classified as small cell (SCLC) or non-small cell (NSCLC) [8,9]. While SCLCs are malignant tumors that can be identified by neuroendocrine features, accounting for nearly 15 percent of lung cancers [9], NSCLCs account for approximately 85 percent of all lung cancers [8,10]. This discrimination reflects the different clinical recognition, disease course, and therapeutic alternatives of the two subgroups [11]. Adenocarcinomas, accounting for 40 percent of lung cancers, are, for unknown reason, the most common histological subtype seen in NSCLCs over the past 25 years, and NSCLCs include multiple cancer types such as squamous cell, large cell, and mixed histotypes, apart from adenocarcinomas [12]. In addition to this, there is a growing body of experimental and clinical proof bring heterogeneity into focus among NSCLC subtypes. Discovering the molecular mechanism underlying cell proliferation and metastasis in NSCLC could critically aid the development of new and more effective molecular-targeted therapeutic procedures specific to subtypes [12,13].
Although much research has been done on lung cancer, the heterogeneity among its subtypes and the complexity in the mechanism of the disease cause a great challenge in clinical oncology. Because the diagnosis of lung cancer today confides in histopathological analysis, molecular markers and imaging, and therefore early diagnosis is very difficult [14,15]. Mathematical models can be adapted to try to estimate the complex dynamics of disease and by developing models for this, simulate the kinds of treatments that are appropriate and effective for patients in the context of personalized medicine [16][17][18]. Mathematical models that are adaptable for processes important in cancer biology will shed light on unknown points in the field of oncology.
In the world, lung cancer is leading cause of cancer deaths. Although there are several successful treatment alternatives, it continues to exist as a highly dangerous disease that has not been fully cured. The treatment techniques for this illness are still restricted some of these techniques are radiotherapy, hormonal therapy, chemotherapy, gene therapy. Although vaccination studies on cancer have been conducted for years, it still seems to need time to achieve a remarkable result. Cells in cancer tumors also change and evolve like living things in nature. Understanding how this process works will make it easier for us to beat cancer from the very beginning.
Judging by the numbers, victory over cancer still seems far away. A person's life-long cancer risk in the USA is 42 percent in men and 38 percent in women. The Cancer Research Foundation in England gives this rate as 54 and 48 percent, respectively.
As of 2015, the number of cancer patients in England reached 2.5 million. This is an annual increase of 3 percent, in other words, 400,000 extra cancer patients observed in five years.
In destroying cancer cells, immune cells play major role. The most important immune system cells are T cells (such as CD4 + T cells, CD8 + T cells), B cells, macrophages and (NK) natural killer cells. Macrophages have an important role in innate and adaptive immune response [19,20]. They are the most abundant cells in the tumor micro-environment. Research shows that macrophages can break down malignant cells very effectively. Macrophages are also very important as the target of innovative oncological treatments. As a therapeutic goal, tumor shrinkage is achieved by closing macrophages with CSF1 receptors, antibodies and small molecule drugs. However, the events occurring in the periphery of the tumor microenvironment are highly complex and change rapidly. As a result, designing treatment methods to cure or treat cancer is an extremely difficult task. Mathematical models contribute greatly to understanding how immune cells and cancer cells interact, and to define tumor-immune dynamics [21]. So, many researchers have benefited from mathematical models to understand the course of cancer. For example, in [22], Banerjee and Sarkar have taken into account the delay-induced tumor-immune system model and malignant tumor growth. In [23], the authors study the mathematical model of tumor-macrophages cells interaction. In addition, based on tumor radio-biologic mechanisms, the authors propose the tumor growth model considering the radiotherapy effect in [24]. They investigate how the re-oxygenation of hypnotic cells and the radio-sensitivity of radiotherapy influence the effect of tumor radiotherapy. In other study [25], the qualitative analysis of tumor cell-immune cell interaction and chemotherapy with (GWN) Gaussian white noises has been researched. Eftimie and Barelle [26] introduced a model for interactions between lung cancer and macrophages with different phenotypes (M1, M2, M1/M2) with an integer-order differential equation system. Sarmah et al. [27] formulated a seven-dimensional mathematical model connecting p53, DNA damage, and autophagy in lung cancer. They also performed both local and global sensitivity analyses along with parameter recalibration analysis to understand the system dynamics. Feng et al. [28] established a mathematical model for predicting malignancy of lung cancer complicated with Talaromyces Marneffei infection and its chest imaging characteristics, and improve clinicians' understanding of the disease.
Fractional calculus has been used only by mathematicians, physicists and engineers for the three centuries. However, in recent years, fractional operations with its applicability have gained also great importance in unusual areas, such as finance [29,30], medicine [31], physics [32], synchronization of chaotic systems [33], atmospheric ocean problems [34], hydrology, signal processing [35], heat diffusion models [36], competing species [37], and biology [38][39][40][41][42][43][44][45][46]. Among them, in [47], two integer order systems have been investigated with fractional calculus, which explain the tumor-immune system dynamics. In [48], a fractionalorder model which propose the growth of tumor with drug application is studied. In [49], a delay model with fractional order for tumor-immune system with treatments has been proposed. A study on a mathematical model of immune response with cell mediated for tumor phenotype heterogeneity has been presented in [50]. The authors of [51] presented a model using fractional derivative, taking into account cell repair, re-population of the cells and radiotherapy process. A numerical and analytical study of the HIV-1 infection of CD4+ T cells constructed with conformable fractional mathematical model was examined in [52]. Veeresha et al. applied the q-homotopy analysis transform method (q-HATM) to the mathematical model of the cancer chemotherapy effect in the sense of Caputo derivative in [53]. In addition, in [54], the authors analyzed the behavior of solution to the system exemplifying model of tumor invasion and metastasis by the help of q-HATM.
There are only a few studies on fractional-order mathematical modeling of lung cancer in the literature. For example, in [55], the authors proposed two modeling approaches to predict lung tumor dynamics as an effect of radiotherapy. They used the real clinical information of non-small cell lung cancer (NSCLC) patients undergoing stereotactic body radiation therapy (SBRT) as the primary treatment method for numerical simulations. In another study, Ghita et al. [31] concluded that feature extraction from modeling a respiratory function through a specific fractional order impedance model could be transposed to lung tumor dynamics.
Biological phenomena and diseases are one of the most common practice areas of fractional-order mathematical models. Considering the studies until the last decade, while constructing mathematical models, it has generally seen that the authors are used the differential equations of integer order. Mathematical model studies created with these equations have proven to be very valuable in explaining the course of diseases. However, it has been observed that the models made with fractional-order differential equations (FODEs) are more compatible with the truth and provide more advantages when compared with integer-order mathematical models. Because most of the biological systems continue to function using the memory, after-affects and hereditary properties. These effects are neglected in integer-order differential equations models; however, the fractional-order models explain these complex phenomena better than these equations. Before presenting our model, we mention a few mathematical models of tumor-immune systems: In [56], a study on tumor and macrophages cells is modeled as follows: where P(t), Q(t), R(t) are the concentrations of macrophages, activated macrophages and tumor cells. r 1 is the growth rate of macrophages cells and k 1 is the carrying capacity. The activation rate of macrophages cells is b. d 1 and d 1 are the natural death rates of P(t) and Q(t), respectively. After active macrophages destroy tumor cells, they become passive at a rate ε 1 . r 2 is the growth rate of tumor cells and k 2 is the carrying capacity. The conversion rate of normal cells to malignant cells is c. When active macrophages cells destroy the tumor cells, the loss rate of tumor cells is a.
Khajanchi et al. [57] have considered the following tumor-immune competitive system originated from prey-predator model presented by Sarkar and Banerjee [58]: Here, M(t), N (t), S(t) denote the densities of tumor cells, CTLs and resting cells at time t, respectively. The growth rate of M(t) is r and k 1 is the carrying capacity. α is the rate at which tumor cells are cleared by prey cells or CTLs. λ≥ 0, is the discrete time delay, β 1 is conversion rate of resting stage to hunting stage of CTLs population and the natural death of CTLs is d 1 . s is the growth rate of the resting cells and k 2 is the carrying capacity. β 2 is the conversion rate of resting cells to CTLs, and d 2 is the natural death rate. Khajanchi et al. [57] neglected that the term q in [58] is the conversion of normal cells to malignant cells as they assumed that the tumor cells are malignant. Moreover, they consider the conversion of resting cells to hunting cells, and the degradation of resting cells due to hunting cells both must be different, not the same.
In [20], Hu and Jang proposed the integer-order mathematical model of tumor-immune cells interconnection to study the effect of CD4+T on tumor dynamics. Their conclusions show that host cells and CD4+ T cells have great importance in fighting tumor cells and slowing their growth. The ordinary differential equations model is described as follows: where t ≥ 0 and r 1 , b 1 , c 1 , a 1 , δ 1 , β 1 , α 1 , μ 1 , δ 2 , β 2 ,α 2 , μ 2 , r 2 , b 2 , δ 3 are the positive constants. x(t), y(t), z(t) and w(t) are tumor cells, C D4 + T cells, cytokines, and host cells, respectively.
In study [59], the following model was reconstructed by replacing the integer-order derivatives in model (1) with Caputo derivatives: with positive initial conditions. Here, t ≥ 0, α (0 < α ≤ 1) and P(t), Q(t) and R(t) denote the macrophages cells, activated macrophages cells and tumor cells, respectively. a, b, c, k 1 , k 2 , r 1 , r 2 , d 1 , d 2 , ε 1 are the positive constants. The biological meanings of these parameters are the same as the meanings of the parameters in the system (1). In the fractional systems, dimensional consistency is a very important tool, in which the units of measurement from the left-and right-hand sides of the equations are coherent. This consistent can be provided by modifying the parameters involved in the right-hand side of the equations, e.g., raising them to power α. In this context, motivated by [20,59], we have proposed the following the fractional-order system: with the initial conditions:  Table 1: In this study, as in [57], it is assumed that the TCs are malignant. So, we neglected the constant which is the rate of conversion of NTCs to malignant cells in Sarkar [58] and The death rate of active macrophages δ 1 The competition coefficient of NTCs on TCs δ 2 The competition coefficient of TCs on NTCs The growth rate of macrophages cells The growth rate of TCs The growth rate of NTCs μ The rate of destruction of TCs due to attack of active macrophages β 1 The conversion rate of macrophages to active macrophages β 2 The degradation of macrophages cells due to active macrophages Mukhopadhyay [56]. In addition, we consider that the conversion of macrophages cells to active macrophages cells and the degradation of macrophages cells due to active macrophages cells both must be different. Also, we removed the parameter ε 1 (in [56]) is the rate of active macrophages that revert back to passive state after attacking the tumor cells (this is not biologically meaningful). It can be considered that these modifications would make the model more realistic. Motivated by the above discussion, in this study, both fractional modeling has been taken into account and real experimental data from patients have been used. In order to fit the parameters and to minimize the mean absolute relative error between the plotted curve for the tumor cells class and the real data provided by lung cancer patients, we have utilized least-squares curve fitting technique (LSCFT). By doing so, we aim to predict the changing of tumor cells and immune system cells over time by using more accurately generated parameters. Meanwhile, dimensional compatibility has been considered in order to better reveal the effect of fractional-order in the proposed fractional-order lung-cancer system. In addition to this, we have aimed to point out the advantages of the fractional-order modeling, taking into consideration the memory trace and hereditary traits which are capable of integrating all past activities and takes into account the long-term history of the system. In this context, it can be seen that the memory trace dynamics are highly dependent on time. When the fractionalorder α is decreased from unit, the memory trace nonlinearly increases from 0. Hence, the fractional-order system dynamics are quite different from the integer-order dynamics. Best of our knowledge, it is thought that there is not a comprehensive study such this paper in the literature that makes the fractional-order model dimensionally consistent, performs the parameter estimation with lung cancer patient data, and considers memory effect/hereditary properties.
The remaining part of this paper is prepared as follows. In Sect. 2, some definitions of a fractional-order derivative (FOD) and some important theorems for FODs are given. In Sect. 3, the existence and uniqueness conditions of the solutions are given. In Sect. 4, stability theorems for the equilibrium points are examined. In Sect. 5, the parameter estimation method to predict the parameters used in the proposed model has achieved. In Sect. 6, the numerical scheme has been given. In Sect. 7, the effects of the memory trace on the behavior of the system (5) are examined. In Sect. 8, to investigate the effects of different parameter values and different values of α on the dynamic behavior of the proposed model, the numerical solutions have been carried out. In Sect. 9, measurement of memory trace for the proposed fractional-order model is investigated. Finally, the general summary is given in Sect. 10.

Preliminaries and definitions
Definition 1 [60] The Riemann-Liouville form of fractional derivative of order α > 0 of a function f : Definition 2 [60] The Caputo fractional derivative of order α > 0 of the function that has been given in Definition 1 is presented as: For the convenience, we use the notation of D α f (t) to represent the Caputo fractional integral Theorem 1 [61,62]. Consider the following fractional-order system: with x ∈ R n and α ∈ (0, 1). The equilibrium points of the system (8) are solutions to the equation f (x) = 0. An equilibrium is locally asymptotically stable if all the eigenvalues λ i (i=1,2"…,n) of the Jacobian matrix J = ∂ f ∂ x evaluated at the equilibrium satisfy On the other hand, if |arg (λ i )| < απ 2 , then the equilibrium point is unstable.

Existence and Uniqueness (E&U) of the solution
With the initial conditions as follows: where The definitions required for E&U (existence and uniqueness) are given as follows: If t > σ ≥ 0, one can write C * σ [0, ] and C σ [0, ]. (10). Proof Using the properties of fractional calculus, we write the system (10) as follows: operating with I α we have Then, let F: C * [0, ] → C * [0, ] be given by Then, This gives us that So, the operator F given by (12) has a unique fixed point. That is, (11) has a unique solution X∈ C * [0, ] . From (11), we have and from which we can deduce that and Therefore, (11) is equivalent to the system (10).

Equilibrium points
To find the equilibrium points of system (5), we write Therefore, we have the following equilibrium points: Now, using Theorems 1 and 2, we obtain the stability conditions of these equilibrium points. But, our aim is to examine the conditions under which the patient can survive [20]. So, we only examine the stability of the equilibrium point E 11 . The stability conditions of other equilibrium points have been omitted.

Local stability of the endemic equilibrium
We now examine the local stability of the endemic equilibrium pointE 11 =(T, A, M, W). The Jacobian matrix of model (5) calculated at equilibrium point E 11 =(T, A, M, W) is given by We obtain the characteristic equation such as: By Theorems 1 and 2, one can deduce the following conclusions.

Theorem 4
The equilibrium point E 11 is locally asymptotically stable if one of the following conditions holds.

Parameter estimation
Parameter estimation (PE) is a very important issue in the identification and validation of many cancer models. By performing PE, the parameters of the proposed model are determined within appropriate bounds using real data and thus, it is tried to ensure that the model is the most appropriate model expressing cancer. The most important advantage of parameter estimation is to produce parameter values specific to the model created, instead of using the In this section, it has been aimed to fit the relevant parameters with the help of the data of a lung cancer patient who was treated at Kayseri Erciyes University hospital using the least-squares curve fitting method. A total of 12 parameters in the model have been fitted with the aid of least-square curve fitting method using the data of lung tumor cells of the patient for 14 days from 15 September-28 September. The average absolute relative error between tumor cells real data belonging to lung cancer patient and the solution of model's for the tumor class is tried to be decreased, and the best fitted parameter values have been obtained. In Fig. 1, the red solid circles show the real tumor values, while the best fitted curve of the model is indicated by the blue solid line. The biological parameters in the model along with the best predicted values obtained through LSCFM are given in Table 2.

Providing real data (sample collection and cell culture)
Tumor tissue samples were collected from the patient with non-small cell lung cancer who had chemotherapy-naive. The diagnosis was verified on routine H&E(hematoxylin and eosin)stained slides by a histopathologist. Informed written consent was obtained from the patient before the sample collection. Ethics committee approval for study was obtained from Erciyes   Table 3).
In the literature, cell lines and animal data are also widely used as experimental data. Cell lines are important resources for identifying predictors of unrestricted growth, eligibility for high-throughput screening, response to therapy, and resistance [66,67]. They are also cheaper, easier to use, and experiments on them are more reproducible. However, currently available cancer cell lines have some limitations. Cell lines are selected from tumor subsets grown under in vitro culture conditions [68]. However, it is often unclear from which cancer subtype these cells originate and to what extent they resemble the tissue of origin [69]. This produces cancer cell lines that do not represent the diversity of human tumors [67,68]. Moreover, cell lines grown only in monolayer cultures lack the heterogeneity observed in cells derived from patient tumors [67,70]. Although the development of patient-derived xenograft cancer models provides some improvement over cancer cell line models, it has some limitations. Vaccination and drug validation time in mice often take a long time, and this time delay limits its applicability to patient therapy. In addition, these models are suitable for a limited number of drug combinations and are not amenable to genetic manipulation such as the introduction of transgenes or knockout studies [68]. In primary cell cultures, cells are isolated directly from the tumor site, and detailed pathology information is also available, allowing the characteristics of the culture to be compared with those of the original tumor. In addition, primary cell culture provides a more accurate prediction of patients' responses to chemotherapy treatments [71]. In this study, in order to represent the origin of tumor profile in cancer patients population, cancer cells were isolated from NSCLC patients.

Numerical scheme for the provided tumor-immune system model in the Caputo derivative sense
To investigate the dynamics of the proposed fractional-order model (2), we implement the Caputo fractional operator. To provide the numerical simulation of the suggested nonlinear fractional-order system, the Adams type estimator-corrector method [72][73][74][75] is used. The following Cauchy-type ODE is taken into account w.r.t. the Caputo operator of order α: where b = 0, 1, ..., n −1, and n = α . Equation (19) can be turned to the Volterra equation: By considering this proposed predictor-corrector scheme associated with the Adam-Bashforth-Moulton algorithm [73] to have the numerical solutions of the proposed model, we can take h = τ/N , t z = zh, and z = 0, 1, ..., N ∈ Z + , by letting z ≈ (t z ) , it can be discretized as follows, i.e., the corresponding corrector formula [76] where Subsequently, the following step is to construct the coincident predictor formula P F q+1 . One can compute the proposed predictor formula as: where

Memory trace and hereditary traits
To examine the behavior of the tumor-immune-host cells model (2), we use the Caputo operator defined in (2). For α, 0 < α ≤ 1 derivative, let the fractional derivative of variable (t) be [77,78]  Utilizing the one of most common numerical methods, the L1 scheme [78][79][80][81], the numerical approximation of the FOD of (t) is One of the most powerful numerical methods for discretizing the Caputo-FOD in time is L1 scheme. The purpose of implementing the L1 scheme to this research study is its memory term and convergence rate. Memory term is also explicitly present in other numerical methods, but this memory integration term is more clearly defined in the L1 scheme. Considering (24) and (25) together, the numerical solution of (24) is as follows: Therefore, the solution of the FOD (fractional-order derivative) can be defined as the difference between the Markov term and the memory trace. The Markov term weighted by the Gamma function is as follows: The memory trace ( -memory trace since it is related to variable (t)) is The memory trace is capable of integrating all past activities and takes into account the longterm history of the system. For α = 1, the memory trace is 0 for any time t. Memory trace dynamics is highly dependent on time. When the fractional-order α is decreased from unit, the memory trace nonlinearly increases from 0. Hence, the fractional-order system dynamics are quite different from the integer-order dynamics.

Numerical simulations and discussion
For the system (5), the numerical solutions are achieved using the Adams-Bashforth-Moulton Predictor-Corrector method for the parameters in Table 2. With the aid of numerical simulations, the effects of different parameter values and different values of α on the dynamic behavior of the model (5) can be identified. The fitted parameter values used for numerical simulations are given in Table 2.
In this part of the study, the variation of each sub-population over time has been simulated for differing values of the fractional parameter α by using the fitted parameter values given in Table 2. In addition, considering the parameters that significantly change the behavior of the cells, graphics have been obtained for different values of these parameters. The dynamical behavior of each state of the proposed tumor-immune system model has shown in Figs. 2,3,4,5,6,7,8,9 for varying values of the fractional-order parameter α and different parameter values. Figure 2  fractional case. Figure 2 shows decreasing attitude of the tumor cells for increasing values of α. Moreover, in Fig. 2, it has been seen that the total number of tumor cells is more in the case of α = 0.7 than in the case of α = 1. Therefore, the fractional order predicts more tumor cells than the prediction obtained in the integer-order case. In addition, when α = 0.7, it is seen that even around 100 days, tumor cells are still seen in the body. However, it is seen that they disappear from the second day in order of integer-order case. This effect of α = 0.7 indicates that the fractional derivatives are in good agreement with the spread of tumor cells. In Fig. 3, we have changed μ, r 2 and keep the other parameters fixed as in Table 2. One of the main reasons for using FODEs is that fractional differential equations have a wider stability region than integer-order equations, that is, fractional-order equations are at least as stable as integer-order equations. Moreover, solutions in fractional-order equations depend on all previous cases [49,59]. It is concluded that Fig. 3 verifies this theoretical result, that is, the stability region of the proposed model has increased as the α decreases from unit. In Fig. 4, the change in the active macrophages cells with respect to time is examined. One can see from this figure that as the fractional-order α reduces from unit, the memory effect of the active macrophages cells increases. Therefore, these cells take longer to be stable for the non-integer cases. In addition, while the active macrophages cells exhibit an unstable behavior when α is unit, the long-term stability of the system when α is not an integer reveals one of the important properties of the fractional-order derivatives. In Fig. 5, similar behavior like the active macrophages cells has been seen. It has observed that as the α decreases from unit, the memory effect of the macrophages cells increases, which means that the peak of the macrophages cells decreases and the macrophages cells reach the stability after a certain time for the fractional cases.
It can be observed from Fig. 6 that as α reduces from the unit, the host cells take longer time to reach the maximum value due to memory effect. In addition, it is seen that they reach the equilibrium point in a longer time in the fractional-order state. In Fig. 7 for increasing values of μ (the destruction rate of tumor cells due to active macrophages cells), the tumor cells decrease under fractional-order case α = 0.9. In Fig. 8, δ 1 appears to have a significant impact on the regression of the tumor cells. It is seen that the number of cancer cells decreases as δ 1 increases. It is seen that the observed changes in Figs. 2,3,4,5,6,7,8 are compatible with the expected phenomena related to a real cancer patient.
Cancer patients have different tendencies during tumor growth which is difficult to capture by integer order derivative as in seen Fig. 2. In Fig. 3, we have observed that the fractionalorder derivative damps the oscillation behavior about the positive equilibrium point. The existence of periodic solutions is associated with cancer models. It implies that the tumor levels may oscillate around an equilibrium point even in absence of any treatment. Such a phenomenon which is known as "Jeff's Phenomenon" has been observed clinically [82] and has arose in many cancer models [58,59]. It is seen that in Fig. 7, as μ increases, the number of tumor cells decreases. From the biological point of view, the rate of destruction of TCs due to attack of active macrophages (μ) can be increased by the macrophage-based therapy. The macrophage-based therapy is more effective at reducing tumor cell proliferation than standard therapies such as radiotherapy and chemotherapy [83].
It is known that tumor cells and normal tissue cells compete for resources and space [84]. In addition, the signaling interactions between the stromal and neoplastic tissues are important in driving tumor cell proliferation [2]. A certain period of time after the tumor begins to invade the area and source of normal tissues, the normal tissues go on the defensive to protect themselves [2,84]. It takes a certain amount of time for the host cells to take action to attack. In Fig. 6, it is seen that while host cells are not active until the about 20th day, they start to compete with the tumor cells after the about 20th day. In addition, it is seen that in Fig. 10, the memory trace is zero until around the 20th day and it has started to increase after that day. The memory effect seen in Fig. 6 is supported by Fig. 10. The obtained results from Figs. 6 and 10 are consistent with real biological phenomena. Figure 6 shows that fractional-order differential equations play the role of time lag or delay term in ordinary differential model and they are naturally related to systems with memory which exists in tumor-immune interactions [59,85]. In Figs. 10 and 11, we have depicted the memory trace effects on populations for different values of fractional order α. As can be seen from Figs. 10 and 11, when α = 1, it is seen that the memory effect is zero for all populations. As α decreases from 1 to 0.7, the effect of the fractional order and the corresponding memory effect appear. Since the memory effect has a significant importance in the tumor-immune interconnection, it is very important to include the memory effect in the system in order to make the most realistic modeling. It can be seen that the memory effect has a great effect on each compartment in the model we are considering. Considering the values before and after the peak values of the populations, it is seen that there is an important relationship between them and the memory trace. When Figs. 2, 3, 4, 5, 6, 7, 8, 9 and  is noteworthy that memory trace is negative when populations decrease, and memory trace is positive when populations increase. Observing such behaviors clearly demonstrates the importance of including memory trace in the system when modeling cancer. These are the results that have not been addressed in any previous study, revealing the biological reality of our model.

Measurement of memory trace for the proposed fractional-order model
As given in detail in Sect. 7, we have numerically integrated the fractional derivatives using the L1 scheme and the solutions of the FODEs for T (t) , A (t) , M (t) , W (t) , have described as the similar way as in (25). Thus, the numerical approximation of the fractional-order By combining (29) and the first equation of system (2), the numerical solution of active macrophages cells A (t) is given by where and Memory trace = By following the same steps, the numerical approximations of the fractional-order derivative of T (t), M (t), W (t) have been achieved. By considering the corresponding results above, numerical analyses have been provided to visually see the effect of memory trace on each sub-population of system (5). Figure 10 and 11 show the effects of memory trace on the dynamics of populations for varied values of fractional order α. As can be seen from Figs. 10 and 11, when α = 1, there is no memory effect in the system. As α decreases from 1 to 0.7, the impacts of the fractional order and the existence of a memory effect become clear. Actual results and estimates of the tumor-immune system interaction can be obtained in the presence of a memory effect. Therefore, the presence of a memory effect is very important for the biological models. When Figs. 2, 3, 4, 5, 6, 7, 8, 9 and Figs. 10 and 11 are examined together, it is seen that there is a significant relationship between the memory effect and the peak values of the behavior of the system. For example in Figs. 2 and 10, it is seen that the memory trace starts to become negative where tumor cells for different values of α begin to decrease after the peak value.
In biological events, when the immune system first encounters with the tumor cells, the memory response occurs and it is effective during a certain period of time [86]. When Figs. 10 and 11 are examined, it is observed that the memory effect approaches zero after a certain time. That is, these results are expected in real cancer phenomena. It can be understood from the results obtained from the graphs that the FODEs reveal the memory effect of the system quite successfully without the need for any other factors.

Conclusions
In this study, a new fractional-order differential equation model for tumor-immune system interaction related to lung cancer has been studied by taking into account the models given by [20], and [59]. We have constructed a new model clarifying the interactions between tumor cells, macrophages cells, active macrophages cells and normal tissue cells to evaluate the role of macrophages and normal cells on tumor growth and regression. We have also taken into consideration the Caputo type fractional derivative instead of the integer-order derivative. By doing so, we have guaranteed that the fractional-order system (5) is dimensionally consistent: The units of measurement from the left-and right-hand sides of the equations agree. It has been achieved by modifying the parameters involved in the right-hand side of the equations, e.g., raising them to power α.
We have showed that the compartmental system has a solution by benefiting from the fixed-point theorem, and we have provided the equilibrium points of the system. We have also investigated the local stability of the endemic equilibrium point. By using the real data of a lung cancer patient from Erciyes University hospital in Turkey, sum of the squares of the difference between the numerical solution to the tumor cells and real data has been minimized and the best fitted curve has been obtained (see Fig. 1). In addition, we have performed mathematical investigation along with parameter recalibration analysis to understand the system dynamics and the best fitted parameter values have been obtained, accordingly. In order to achieve numerical simulations of model (5), population dynamics of tumor-immune system have been visualized with the help of different fractional orders and the fitted parameters using Adams-Bashforth-Moulton method. Numerical simulations have been used to investigate the effects of the parameters on the growth and regression of the tumor cells.
Although there have been many studies that discuss the tumor-immune interaction in the literature, our model differs from them both biologically and mathematically in terms of investigating the relationship between tumor, macrophages, active macrophages and normal tissue cells. Normal tissue cells (NTCs) and macrophages cells play an significant role in the prevention of tumor growth along with the production mechanism of macrophage cells. According to the results, if the host cells are more competitive, then the tumor cells undergo a considerable loss. In addition, if tumor growth rate (r 2 ) is small, then cancer cells can be eradicated in a shorter time. When the simulation results have been examined, it has been observed that as α changes, the patient's immune system cells and the number of tumor cells also changes significantly. Moreover, it is seen that the oscillations increase as the fractionalorder changes from zero to one. Cancer has a different mode of progression when the tumor starts to appear in the body. With the help of fractional-order equations, as can be seen from the obtained results, this different structure of cancer can be better explained.
It is a very useful method to use experimental data to test the accuracy of the established mathematical models. But, finding the required data is often difficult. Among the studies on the mathematical model of lung cancer, there are studies that benefit from experimental data. These data are mostly cell lines and animal data. However, it is seen that there are not enough studies in the literature on fractional modeling using experimental data on lung cancer. In this context, both fractional modeling has been taken into account and real experimental data from patients have been used in this study. In order to fit the parameters and to minimize the mean absolute relative error between the plotted curve for the tumor cells class and the real data provided by lung cancer patients, we have utilized least squares curve fitting technique. By doing so, we aim to predict the changing of tumor cells and immune system cells over time by using more accurately generated parameters. In addition to this, we aim to point out the advantages of the fractional order modeling, taking into consideration the memory trace and hereditary traits.
If cancer specialists can change the parameters of the cancer-immune system relationship by applying different treatment methods (such as macrophage-based therapy, BMT therapy, hormone therapy, and immunotherapy) according to the results obtained from the mathematical models, it may be possible to stabilize the tumor cells at a certain value or to terminate the cancer. Therefore, such examinations provide other scientists as well as cancer specialists studying the tumor progression with insight into the control of the cancer and may help develop further treatment strategies. In this context, this study will shed light on future studies.
In future works, one can consider that cancer stem cells can be incorporated into the model in which will be developed.
of Turkey).The authors are very much thankful to the reviewers and editor for their valuable suggestions that have improved the quality of the manuscript.

Conflicts of interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Consent for publication Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.