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.

Introduction

Cancer is one of the most ubiquitous genetic diseases[1,2]. It is known that the disease arises because of mutations in the cancer susceptibility genes[3]. In order to understand the mechanism of human breast cancer various methods have been advanced. Loeb et al.[4] proposed the existance of a mutator phenotype that acts as a mechanism in tumor processes. Tomlinson and Bordmer[5-7] investigated the mutator phenotype hypothesis and discovered that the selection for clonal expansion of intermediate cells.

Mathematical models describe the fundamental principles of the population genetics and evolutionary mechanisms that govern tumor initiation and progression in cancer biology[8,9]. It was discovered that the dynamics of tumorgenesis is determined by a number of factors including mutation, selection, and tissue types[10-13]. It is also known that an increase in the postmenopausal and estrogen receptor-positive cases may increase the rate of breast cancer among women[14]. Several attempts have been made for the description of breast cancer in a variaty of perspectives. Enderling et al.[15] presented a model of the growth and invasion of a solid tumour in a domain of breast tissue and proposed a scheme for the surgery and radiation treatment of the tumour. Enderling et al.[16] applied surgery, as well as adjuvant external beam and targeted intraoperative radiotherapies, and used a model to identify different sources of local recurrence to analyse their prevention. Simmons et al.[17] presented a brief overview of breast cancer, focussing on its heterogeneity and the role of modelling and simulation in teasing apart the underlying biophysical processes. Nave et al.[18] described the treatment of breast cancer by a model involving nonlinear ordinary differential equations with a hidden hierarchy. Interested readers can follow also some recent works in[19-25].

Fractional calculus is one of the most intensively developing branches of mathematical analysis and deals with derivatives and integrals of arbitrary order[26]. Fractional differential equations have motivated considerable attention in several branches of science[27-32]. In applied sciences, memory properties have been widely found in many complex phenomena. The use of fractional derivatives, instead of integer ones, can potencially lead to better results since one has an extra degree of freedom. Indeed, due to their intrinsic nonlocal property, the fractional differential equations have been successfully used to describe phenomena or processes with memory and hereditary properties in physics, chemistry, biology, and economy. Readers can refer to[33,34]. Fractional calculus is an excellent tool for modeling materials and processes with memory and hereditary properties and, in particular, electrochemical problems[35]. Fractional differentiation and integration operators are also used for generalizing the diffusion and wave equations[36,37] or, more recently, of the temperature field problem in oil strata[38]. In what concerns cancer we can mention Farayola et al.[39], that simulated a radiotherapy cancer treatment process with radiobiological factors, and Valentim et al.[40], that propose a multistep exponential model with a fractional variable-order representing the evolution history of a tumor.

The Legendre polynomials are useful mathematical tools in fractional calculus. Cao et al.[41] proposed a numerical algorithm based on the shifted Legendre polynomials to solve the fractional governing equations of polymethyl methacrylate in the time-domain. Zhijun et al.[42] derived the operational matrix of fractional integration for Legendre polynomials to solve initial value problems of nonlinear fractional differential equations. Wang and Chen[43] used the shifted Legendre polynomials for the dynamic analysis of viscoelastic pipes conveying fluid with variable order fractional model. Xiao et al.[44] presented a finite-time empirical Gramians, constructed from impulse responses by solving block tridiagonal linear systems, to generate an approximate balanced model. Sun et al.[45] designed an algorithm based on improved Legendre orthonormal basis for solving second-order boundary value problems. Hesameddini and Shahbazi[46] approximated the unknown functions based on the two-dimensional shifted Legendre polynomials operational matrix method for the numerical solution of two-dimensional fractional integral equations. Guorong et al.[47] proposed the Legendre orthogonal polynomials to calculate the acoustic reflection and transmission coefficients at liquid/solid interfaces. Singh et al.[48] obtained numerical algorithms by using the Legendre, Galerkin and Legendre wavelet collocation methods for solving one phase moving boundary problem with conduction and convection effects. Rakhshan and Effati[49] developed the generalized Legendre polynomials and derived a general procedure for solving nonlinear autonomous fractional differential equations with time varying delay. Heydari et al.[50] presented a numerical method based on the discrete Legendre polynomials and the collocation scheme for solving nonlinear space-time fractional KdV-Burgers-Kuramoto equation. Kuznetsov[51] used the Legendre polynomials and trigonometric functions for the numerical solution of the Ito stochastic differential equations when approximating multiple Ito and Stratonovich stochastic integrals based on generalized multiple Fourier series. Dehghan[52] considered the shifted Legendre polynomials for class of variable order fractional functional boundary value problems.

In this work, an optimization method based on the generalized shifted Legendre polynomials (GSLP), operational matrix of derivatives and Lagrange multipliers is proposed for solving a fractional order breast cancer competition model (F-BCCM). The algorithm transforms the problem into a system of nonlinear algebraic equations with unknown coefficients, parameters and Lagrange multipliers. The optimal solution of the problem is obtained by solving an algebraic system of nonlinear equations. The results show that we can achieve the approximate solutions by employing only a few number of the basis functions. Moreover, the proposed approach can be adopted for solving other classes of fractional order problems.

The rest of the paper is structured as follows. Section 2 formulates the F-BCCM and the mathematical concept of fractional order Caputo derivative (F-CD). Section 3 discusses the shifted Legendre polynomials (SLP), as well as the GSLP and its operational matrix of derivatives and function approximation. Section 4 addresses the convergence analysis of the proposed method. Section 5 describes the method for finding the solution of the proposed problem. Section 6 discusses several illustrative examples. Section 7 analyses the epidemiologic and clinical relevance of the problem. Finally, Section 8 gives the main conclusions.

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).

Table 1

The parameters of the BCCM model.

Parameters	Description	Value
k_1	Normal rate of division for C cells	0.75 day^-1
k_2	Normal rate of division for T cells	0.514 day^-1
q	Normal rate of division for H cells	0.70 day^-1
M_1	Carrying capacity of C cells	2.27×10^6 cells
M_2	Carrying capacity of T cells	2.27×10^7 cells
M_3	Carrying capacity of H cells	2.5×10^7 cells
γ_1	Cancer stem cell death rate due to immune cell response	3×10^-7 cell^-1 day^-1
γ_2	Tumor cell death rate due to immune cell response	3×10^-6 cell^-1 day^-1
γ_3	Immune cell death rate due to tumor cell response	1×10^-7 cell^-1 day^-1
p_1	Rate at which estrogen helps to proliferate cancer stem cells	600 cell day^-1 (pg/mL)^-1
p_2	Rate at which estrogen helps to proliferate tumor cells	0 cell day^-1 (pg/mL)^-1
p_3	Rate at which healthy cells are lost to DNA mutation by estrogen presence	100 cell day^-1 (pg/mL)^-1
a_1	Number of C cells at which the rate of absorption is at half its maximum	12M_1 cells
a_2	Number of T cells at which the rate of absorption is at half its maximum	12M_2 cells
a_3	Number of H cells at which the rate of absorption is at half its maximum	12M_3 cells
n_1	Normal death rate of tumor cells	0.01 day^-1
n_2	Normal death rate of immune cells	0.29 day^-1
δ	Healthy cell death rate due to competition with tumor cells	6×10^-8 day^-1 cell^-1
s	Source rate of immune cells	1.3×10^4 cell day^-1
ρ	Immune cell response rate	0.20 day^-1
ω	Immune cell threshold	3×10^5 cells
u	Rate of immune suppression by estrogen	0.20 day^-1
v	Estrogen threshold	400 pg mL^-1
τ	Continuous infusion of estrogen	2000 pg mL^-1 day^-1
μ	Washout rate of estrogen by the body	0.97 day^-1
d_1	Absorption rate of estrogen by cancer stem cells	0.01 day^-1
d_2	Absorption rate of estrogen by tumor cells	0.01 day^-1
d_3	Absorption rate of estrogen by healthy cells	0.01 day^-1

This paper formulates an alternative representation of the BCCM considering the F-CD. The F-CD of order , with respect to t is given by[53,54]:where denotes the Gamma function . 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 , we can write

Rewriting the model (2.1) in terms of the F-CD, we obtainwhere , , represents the F-CD of order .

We must note that the dimension of the left-side equations of model (2.4) is (time), . Nonetheless, a close inspection of the right-hand sides shows that the quantities , , q, , , , , , , , , , s, , u, , , , and have the dimension (time) 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

This is the system that we will actually use for modeling our problem. Note that in the limit case , , 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 , can be determined with the recurrence formulawhere and . To use Legendre polynomials in the interval [0, 1] we have to define the SLP by means of the change of variable . The SLP, , can be denoted by . Therefore, follows the relationship:where and . The analytical form of the SLP of degree j, , is given by:

Note that and .

A given function g(t) can be expressed using the SLP as followswhere is an -order column vector including the basis functions. The following selections for and are considered asand

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 , the GSLP, , are constructed through a change of variable. Therefore, is transformed into , , in the SLP and are defined bywhere denotes the control parameters. If , 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:whereandwith standing for the control parameters.

The -th F-CD of is given bywhere (the operational matrix of F-CD of order ) is as follows

Function approximation

Let us consider . We introduce

Let be the best approximation of u(t). Therefore, we have

Since , there exist the unique coefficients , such thatwhere and are defined in Eqs. (3.9) and (3.12).

Convergence analysis

The classical Weierstrass theorem (see[57], Theorem 7.26) states that polynomials are dense in C(I), the space of all continuous complex functions on the closed interval with the supremum norm. In other words, the set of all finite linear combinations of the functionsis dense in C(I). This can be expressed by saying that the functions (4.1) span C(I). The following question arises naturally, namely, if , under what conditions, then is it true that the functionsspan C(I)? The neat answer is that the functions (4.2) span C(I) if and only if .

Theorem 1

(The Müntz-Szasz Theorem):[58] Suppose that and let X be the closure in C(I) of the set of all finite linear combinations of the functions

If , then .

If and if , , then X does not contain the function .

Remark 1

Let be the set of all finite linear combinations of the functions defined in (4.2), where and let be the set of the functionswhere is a sequence of real numbers, with for all . It is obvious that and, hence, , where and are the closures of the sets and in C(I), respectively. In view of Theorem 1, it is apparent that the set is dense in C(I), i.e., .

Similarly, if is the set of set of all finite linear combinations of the functionswhere is a sequence of real numbers, with for all , then, in view of Theorem 1, the set is dense in C(I), i.e., .

The following two theorems are immediate consequences of Theorem 1 and, therefore, we omit the details.

Theorem 2

Suppose that and let be the set of all finite linear combinations of the functions defined in (4.4).

If , then , where is the closure of the set in C(I).

If , and if , , then Y does not contain the function .

Theorem 3

Suppose that and let be the set of all finite linear combinations of the functions defined in (4.5). 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.

If , then , where is the closure of the set in C(I).

If , and if , , then Z does not contain the function .

Theorem 4

Let be a continuous function. Then, for every there exists a generalized polynomials[59], , such that

Proof

Let be a fixed real number. Chose a sequence of real numbers with , for all and . 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 be a continuous function. Then, for every , there exists a GSLP, , such that

Let be a fixed real number. Chose a sequence of real numbers with , for all and . 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 followswhere , , , , and , , are undetermined vectors including the free coefficients and control parameters, and and , , are defined in Eqs. (3.9–3.14). Regarding (3.16), one has

Moreover, we approximate the initial conditions given in Eq. (2.5) via the GSLP as follows

Now, we define the residual function by using Eq. (2.5) and Eqs. (5.1–5.2), so that

Meanwhile, from Eqs. (2.5) and (5.3), we have

We can generate the 2-norm of the residual functions aswhere and , , 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 problemsubject to Eq. (5.5), where is the objective function.

To solve this minimization problem, we assume thatwhere 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 , , , , , 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:

Numerical simulations

In this section, two examples are considered to illustrate the applicability of the proposed technique when solving the F-BCCM (2.5). The numerical simulations are based on the packages Maple 18 (with accuracy 40 decimal digits) and Matlab.

Example 1

Consider the F-BCCM (2.5) and the parameters listed in Table 1. The algorithm is used to solve the F-BCCM (2.5) with , , , and , for different values of and , . The approximate solutions of the five state variables , with , , , for , , , , and also , , , , are shown in Figs. 1 and 2, respectively. The plots of the approximate solutions with , , for , , , , and also , , , , are illustrated in Figs. 3 and 4, respectively. The runtime of the proposed method and the optimal values of the residual function with different choices of , are reported in Tables 2 and 3, respectively. According to Figs. 1, 2, 3, 4, the number of cancer stem and tumor cells and estrogen, C, T and E, decreases, while the healthy and immune cells, H and I, increase. This effect is associated with better prognosis and longer patient survival. The obtained results with the proposed method provide a meaningful solution by applying even a small number of basis functions. From Table 2 we verify that if we select large values of , , then we pose an higher computational load.

Figure 1

Evolution of vs time for the F-BCCM with , , , , , , and for Example 1.

Figure 2

Evolution of vs time for the F-BCCM with , , , , , , and for Example 1.

Figure 3

Evolution of vs time for the F-BCCM with , and , , , , and for Example 1.

Figure 4

Evolution of vs time for the F-BCCM with , and , , , , and for Example 1.

Table 2

The runtime (in seconds) of the proposed method with different choices of , , for Example 1.

Case	m_1	m_2	m_3	m_4	m_5	CPU time	CPU time
						η_1=0.10,η_2=0.15,η_3=0.20,η_4=0.25,η_5=0.30	η_1=0.98,η_2=0.99,η_3=0.97,η_4=0.98,η_5=0.99
1	2	3	3	3	4	15.23	15.29
2	3	4	4	4	5	19.74	19.76

Table 3

The optimal values of the residual function with different values of , , for Example 1.

Case	m_1	m_2	m_3	m_4	m_5	Residual function	Residual function
						η_1=0.10,η_2=0.15,η_3=0.20,η_4=0.25,η_5=0.30	η_1=0.98,η_2=0.99,η_3=0.97,η_4=0.98,η_5=0.99
1	2	3	3	3	4	7.4167E-07	2.9968E-07
2	3	4	4	4	5	8.2858E-08	7.4587E-09

Example 2

Let us consider the F-BCCM (2.5) and the parameters values given in Table 1. We consider , , , and for different values of and , . The new method is applied to obtain the numerical solution of the F-BCCM (2.5). The approximate solutions of the five state variables with , and for , , , , and also , , , , are shown in Figs. 5 and 6, respectively. The plots of the approximate solutions with , for , , , , and also , , , , are illustrated in Figs. 7 and 8, respectively. The runtime of the proposed method and the optimal values of the residual function with different choices of , , are reported in Tables 4 and 5, respectively. According to Figs. 5, 6, 7, 8, the number of cancer stem and tumor cells and estrogen, C, T and E, decreases, while the healthy and immune cells, H and I, increase. This effect is associated with better prognosis and longer patient survival. Also, from Table 4 we conclude that if we choose large values of , , then leads to an higher computational load.

Figure 5

Evolution of vs time for the F-BCCM with , , , , , , and for Example 2.

Figure 6

Evolution of vs time for the F-BCCM with , , , , , , and for Example 2.

Figure 7

Evolution of vs time for the F-BCCM with , , , , , and for Example 2.

Figure 8

Evolution of vs time for the F-BCCM with , , , , , and for Example 2.

Table 4

The runtime (in seconds) of the proposed method with different choices of , , for Example 2.

Case	m_1	m_2	m_3	m_4	m_5	CPU time	CPU time
						η_1=0.08,η_2=0.17,η_3=0.13,η_4=0.11,η_5=0.23	η_1=0.96,η_2=0.98,η_3=0.99,η_4=0.96,η_5=0.95
1	4	4	5	7	7	26.41	26.43
2	6	6	6	8	8	35.08	35.18

Table 5

The optimal values of the residual function with different choices of , , for Example 2.

Case	m_1	m_2	m_3	m_4	m_5	Residual function	Residual function
						η_1=0.08,η_2=0.17,η_3=0.13,η_4=0.11,η_5=0.23	η_1=0.96,η_2=0.98,η_3=0.99,η_4=0.96,η_5=0.95
1	4	4	5	7	7	3.9316E-11	4.8842E-12
2	6	6	6	8	8	8.0023E-14	1.1164E-14

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 , , 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-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 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.