Optimal solution of the fractional order breast cancer competition model

In this article, a fractional order breast cancer competition model (F-BCCM) under the Caputo fractional derivative is analyzed. A new set of basis functions, namely the generalized shifted Legendre polynomials, is proposed to deal with the solutions of F-BCCM. The F-BCCM describes the dynamics involving a variety of cancer factors, such as the stem, tumor and healthy cells, as well as the effects of excess estrogen and the body’s natural immune response on the cell populations. After combining the operational matrices with the Lagrange multipliers technique we obtain an optimization method for solving the F-BCCM whose convergence is investigated. Several examples show that a few number of basis functions lead to the satisfactory results. In fact, numerical experiments not only confirm the accuracy but also the practicability and computational efficiency of the devised technique.


Fractional order breast cancer competition model
The mathematical theory and concepts adressing the evolution of diseases and epidemics have been advanced during the last decades. Usually, these formulations consider that all people in a community starts as healthy and that later some of them may be diagnosed with breast (cancer). Abernathy et al. 19 described the dynamic behavior of giving up BCCM. The model was discussed analytically and considering cancer stem cells (C), tumor cells (T), healthy cells (H), immune cells (I), and excess strogen (E). The proposed model is given as In Table 1, we list the description and baseline values of the parameters in the system (2.1). This paper formulates an alternative representation of the BCCM considering the F-CD. The F-CD of order 0 < η ≤ 1 , with respect to t is given by 53,54 : where Ŵ(·) denotes the Gamma function Ŵ(z) = ∞ 0 t z−1 e −t dt, z > 0 . In expression (2.2), the convolution integral represents the memory effect embedded in the fractional derivative. Indeed, the fractional derivative uses the previous values of f (t), and captures the long memory effect of the dynamics. From (2.2), for any r ∈ N , we can write C 0 (t) = C(0), T 0 (t) = T(0), H 0 (t) = H(0), I 0 (t) = I(0), E 0 (t) = E(0). www.nature.com/scientificreports/ Rewriting the model (2.1) in terms of the F-CD, we obtain where C 0 D η i t , i = 1, 2, 3, 4, 5 , represents the F-CD of order 0 < η i ≤ 1. We must note that the dimension of the left-side equations of model (2.4) is (time) −η i , i = 1, 2, 3, 4, 5 . Nonetheless, a close inspection of the right-hand sides shows that the quantities k 1 , k 2 , q, γ 1 , γ 2 , γ 3 , p 1 , p 2 , p 3 , n 1 , n 2 , δ , s, ρ , u, τ , µ , d 1 , d 2 and d 3 have the dimension (time) −1 and, therefore, we need to modify the right-hand sides to match the dimensions. The most straightforward way of doing this gives the following model (2.4) www.nature.com/scientificreports/ This is the system that we will actually use for modeling our problem. Note that in the limit case η i −→ 1 , i = 1, 2, 3, 4, 5 , the system (2.5) reduces to classical one (2.1) 55,56 .

Basis functions
Hereafter, we introduce two classes of the basis functions, namelly the SLP and GSLP, which will be used in approximating solutions of the F-BCCM (2.4).
Approximation by the shifted Legendre polynomials. The Legendre polynomials, defined on the interval [−1, 1] , can be determined with the recurrence formula where P 0 (t) = 1 and P 1 (t) = t . To use Legendre polynomials in the interval [0, 1] we have to define the SLP by means of the change of variable t → 2t − 1 . The SLP, P j (2t − 1) , can be denoted by L j (t) . Therefore, L j (t) follows the relationship: where L 0 (t) = 1 and L 1 (t) = 2t − 1 . The analytical form of the SLP of degree j, L j (t) , is given by: Note that L j (0) = (−1) j and L j (1) = 1. A given function g(t) can be expressed using the SLP as follows where Q n (t) is an (n + 1)-order column vector including the basis functions. The following selections for R T and Q n (t) are considered as and Approximation by the generalized shifted Legendre polynomials. Let us define a set of basis functions based on the GSLP to obtain an efficient solution of (2.4). For m ∈ N , the GSLP, L m (t) , are constructed through a change of variable. Therefore, t i is transformed into t i+α i , i + α i > 0 , in the SLP and are defined by where α i denotes the control parameters. If α i = 0 , then the GSLP coincide with the classical SLP.
By using the GSLP, the functions C(t), T(t), H(t), I(t) and E(t) can be expressed in the following form: www.nature.com/scientificreports/ where and with α i j standing for the control parameters.
(the operational matrix of F-CD of order η i ) is as follows Function approximation. Let us consider X = L 2 [0, 1] . We introduce Let u * (t) ∈ Y m 1 be the best approximation of u(t). Therefore, we have

Convergence analysis
The classical Weierstrass theorem (see 57 Remark 1 Let C (I) ⊂ C(I) be the set of all finite linear combinations of the functions defined in (4.2), where 0 < 1 < 2 < · · · and let C β (I) be the set of the functions where {β i } i∈N is a sequence of real numbers, with i + β i ≥ 0 for all i ∈ N . It is obvious that C (I) ⊂ C β (I) and, hence, C (I) ⊂ C β (I) , where C (I) and C β (I) are the closures of the sets C (I) and C β (I) in C(I), respectively. In view of Theorem 1, it is apparent that the set C β (I) is dense in C(I), i.e., C (I) = C(I).
Similarly, if C γ (I) is the set of set of all finite linear combinations of the functions where {γ i } i∈N is a sequence of real numbers, with i + γ i ≥ 0 for all i ∈ N , then, in view of Theorem 1, the set C γ (I) is dense in C(I), i.e., C γ (I) = C(I).
The following two theorems are immediate consequences of Theorem 1 and, therefore, we omit the details. We now investigate the convergence of the proposed method. We first discuss the convergence analysis of the GP expansion by means of the following theorem.
Theorem 4 Let f : [0, 1] → R be a continuous function. Then, for every ǫ > 0 there exists a generalized polynomials 59  www.nature.com/scientificreports/ Proof Let ǫ > 0 be a fixed real number. Chose a sequence {β i } i∈N of real numbers with i + β i ≥ 0 , 0 < i + β i < i + 1 + β i+1 for all i ∈ N and ∞ n=1 1 n+β n = ∞ . Applying Theorem 1, we get the desired result. This completes the proof.
For discussing the convergence analysis of the GSLP expansion we consider the following theorem.
Theorem 5 Let f : [0, 1] → R be a continuous function. Then, for every ǫ > 0 , there exists a GSLP, L m 1 (t) , such that Proof Let ǫ > 0 be a fixed real number. Chose a sequence {γ i } i∈N of real numbers with i + γ i ≥ 0 , 0 < i + γ i < i + 1 + γ i+1 for all i ∈ N and ∞ n=1 1 n+γ n = ∞ . Applying Theorem 1, we get the desired result. This completes the proof.

The proposed strategy
In this section, we design a matrix approach by using the GSLP to solve the problem generated in Eq. (2.5). To carry out this method, we approximate C(t), T(t), H(t), I(t) and E(t) by the GSLP basis as follows  www.nature.com/scientificreports/ where Q i and i , i = 1, 2, . . . , 5 , are the free coefficients and control parameters, respectively. To obtain the optimal values for the free coefficients and control parameters, we consider the following optimization problem subject to Eq. (5.5), where M is the objective function.

(5.5)
To solve this minimization problem, we assume that where denotes unknown Lagrange multipliers. In order to obtain the extremum, the necessary conditions are: Finally, by solving Eq. (5.9) using a software package to calculate the components Q 1 , Q 2 , Q 3 , Q 4 , Q 5 , i and , we obtain the approximate solutions C(t), T(t), H(t), I(t) and E(t) of Eq. (3.5). In this study we used Maple 18 (with 20 digits precision) for the above extracted system and also for all numerical simulations. The step-by-step algorithm of the new technique is summarised as follows:

Remark 2
Abernathy et al. 19 presented a system of five ordinary differential equations which consider the population dynamics among cancer stem, tumor, and healthy cells. They described (i) the effects of excess estrogen and the body's natural immune response on the aforementioned cell populations and (ii) the global dynamics of the F-BCCM (2.5), with η i = 1 , i = 1, 2, 3, 4, 5 , along with various submodels by employing a variety of analytical methods. In this paper we consider the integer order model (2.1) proposed initially in 19 and generalize it to fractional order (2.5). Additionally, we introduced a new basis function (i.e., the GSLP) for solving the F-BCCM (2.5) and obtained optimally the unknoun coefficients and parameters. From the results we verify not only that the approximation by the new method provides good approximate solutions, but also that it is superior to the one proposed by Abernathy et al. 19 .

Epidemiologic and clinical relevance
Tumor cells proliferate abnormally and gradually undergo changes that induce the growth and development of cancer with an high mutation rate and spread, leading to tumor progression. A small number of cancer stem cells can be the source of cancer and cause recurrence, metastasis and resistance to treatment. Presently, stem cells are being targeted for cancer treatment, so that with a lower number of cancer cell we can expect better prognosis. On the other hand, healthy cells must exhibit a normal proliferation and function. In other words, they do not undergo aberrant proliferation and malignant changes. Immune cells play a key role in defending the body against foreign agents and deformed cells such as those of cancer. An high number of immune cells, in particular of T lymphocytes, correlates with a better prognosis for the patient. Epidemiological studies revealed the pivotal role of estrogen in the initiation and progression of breast cancer 60 . The hormone has also influence in the mechanism of some drugs for treatment, since they inhibit estrogen and the duration of exposure to the estrogen increases the risk of breast cancer [61][62][63] . During the progression of the breast cancer there are evidences that an increasing number of immune cell infiltrate in tumor parenchyma, including the cytotoxic CD8 + T, CD4 + T helper, B, macrophages and dendritic cells, the natural killer cells, and cytokines, such as interferons, interleukins, chemokines and growth factors 64,65 . Immune cells contain estrogen receptors and are regulated by estrogens as well. Therefore, strogens could influence immune cells in breast cancer 64 . Parenchymal and stromal cells of breast may be accessible to several immune cells subtypes that lead to decreasing the tumor cells and reducing tumor growth 64 . Therefore, from the medical point of view, the numerical results show that the approximate solution is coherent with real-word experience. From the numerical viewpoint, we must highlight that there is a basic difference between the proposed approach and other spectral methods. In fact, the main idea of spectral methods (based on the Legendre, Chebyshev, Lagrange and Jacobi polynomials) is to express the solution of a differential equation as a sum of the basis functions and then to choose the coefficients in order to minimize the error between the numerical and exact     www.nature.com/scientificreports/ solutions, in some suitable sense. To determine the coefficients, three main techniques, are commonly employed, namely the Galerkin, tau and collocation methods. In the present case, the residual function and its 2-norm are employed for converting the problem to an optimization one, so that the unknown parameters are obtained optimally. As a result, the necessary conditions of optimality are derived in the form of a system of nonlinear algebraic equations with unknown parameters. It is also worth mentioning that approximating any arbitrary smooth function by the eigenfunctions of singular Sturm-Liouville problems, such as the Legendre, Chebyshev, Hermite, Lagrange, Laguerre or Jacobi polynomials, has ′′ spectral accuracy ′′ . This means that the truncation error approaches zero faster than any negative power of the number of the basis functions used in the approximation, as that number tends to infinity. Consequently, these basis functions are not the most adequate for approximating non analytic functions, in contrast which occurs when using the GSLP which prove to be much more efficient.

Conclusion
This paper developed and analyzed the GSLP method for solving the F-BCCM. The proposed approach is based on the operational matrices of the GSLP and the Lagrange multipliers. By adopting the GSLP basis and operational matrices of F-CD, the problem was reduced to the solution of a system of algebraic equations. The convergence analysis for the new algorithm was also carried out. Two numerical examples illustrate the ability and reliability of the algorithm. This method shows that with fewer number of basis functions we can obtain the approximate the solutions. This model and algorithm can be further explored to develop in silico studies of the dynamics and cancer problems. www.nature.com/scientificreports/ 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:// creat iveco mmons. org/ licen ses/ by/4. 0/.