Numerical Simulation of Nonlinear Dynamics of Breast Cancer Models Using Continuous Block Implicit Hybrid Methods

Abstract

In the search for causes and cures of cancer diseases, many mathematical models developed have resulted in systems of nonlinear stiff ordinary differential equations. With these models, many numerical estimates of biological knowledge of the parameters have been obtained, a number of phenomena interpreted, and predictions were made in order to gain further knowledge of cancer development and possible treatment. In this study, numerical simulations of the models were performed using continuous block implicit hybrid methods and the results obtained support the theoretical and clinical findings. We analyzed the interactions among the various tumor cell populations and present the results graphically. From the graphical representation of results, one can clearly see the effects of all the tumor cell populations involved in the competition, as well as the effects of some treatments by the applications of some therapeutic agents which have been heavily used in the clinical treatments of breast cancer. The treatments in the past were mostly conventional chemotherapies, which were used either singly (alone) or in combination with each other or other therapies, and all played vital roles, except for the side effects that these therapies incur in normal tissues and organs. Thus, from recent research works, it is now clear that in many cases they do not represent a complete cure. Therefore, the need to address not only the preventative measures of breast cancer, but also more successful treatment, is clear, and can be successfully achieved to increase the survival rate of breast cancer patients.

1. Introduction

Flow of blood in the body network systems has a profound influence on the efficient transport of nutrients and oxygen to the cells and the transport of metabolic waste products away from the same cells. We are sometimes faced with pathological conditions in which some specific tissues are prevented (deprived) from getting adequate blood and oxygen supply (hypoxia) [1]. In such a situation where disagreement between oxygen supply and its demand at the cellular level exist, which can be due to many reasons—for example, strangulation, high intake of carbon monoxide, cardiac failure or vasculature arrest, as seen in tumor cells [2,3]—it requires special attention. It is obvious that a tumor environment is always characterized by an increase of specialized cells called neoplastic cells which escape from the primary tumor very early in their development and grow rapidly, creating poor blood supply and a deficiency in oxygen supply. Hence, cells of such regions are considered to be drug-resistant due to lack of adequate blood flow, oxygen supply and delivery of drugs via the circulatory systems. Cancer, which is one of the most dangerous killer diseases worldwide, is a class of diseases resulting from unregulated cell growth which has great ability to spread or invade other parts of the body (see [4]). The spread of cancer is just like that of anomalous diffusion, which is observed in porous media where the process of flow does not follow the Fickian and Darcy law (see Atangana and Jain [5]). They can be classified into carcinoma, sarcoma, lymphoma, leukemia, germ cell tumor, blastema, etc., based on the origin of the tumor cells, and can rapidly metastasize [4]. There are many causative factors of cancer diseases, some of which are known: for example, cigarette smoking, abuse of alcohol and poor diet all cause lung cancer, radiation causes sarcoma and leukemia, while carcinoma is a cancer of the epithelial cells that includes breast cancer, prostate cancer, lung cancer, pancreas and colon cancers, to mention only a few. Beside the few factors mentioned, increased estrogen also helps in the growth of cell populations, especially in women of all ages and men of mostly sixty-seven years of age and above [6]. Estrogen can also destroy DNA, which causes healthy epithelial cells to convert into a malignant tumor, just like the action of a carcinogen.

Breast cancer, which is the main concern of this study, is a malignant tumor which begins in the breast cells; that is, a cancer-cell group which can grow into surrounding tissue [6]. When the DNA is damaged, breast cells become out of control, and then the cancer cell keeps using the same damaged DNA to produce new abnormal cells which are of no use to the body systems. This results into a single modified cell gradually growing into a tumor. Several studies have determined that there are multiple types of breast cancer, some of which are undetectable. For example, ductal carcinoma in situ (DCIS) is non-invasive. This is the most common type of breast cancer, especially in women. DCIS means that the cells are still inside the ducts of the breast, yet to spread to other parts of the breast. This can grow into the fatty tissue of the breast and spread to other parts of the body through the lymphatic system and bloodstream, and then become invasive, which is referred to as invasive ductal carcinoma (IDC). The second type is lobular carcinoma in situ (LCIS), which begins in the milk-producing glands of the breast. This can become dangerous when spread to other parts of the body, which can then be referred to as invasive lobular carcinoma (ILC) and is very difficult to detect even when using mammography (see [6]).

Recently, some researchers have devoted their academic attention to understand the spread and seek accurate and proper medical help to fight cancer hazards against mankind. Therefore, for proper understanding of the development of cancer diseases by biologists and mathematicians, mathematical modeling has been considered as a tool to design some experimental data and clinical results which provide proper interaction among the cancer population dynamics (see Valle et al. [7]). The majority of the mathematical models considered resulted in stiff nonlinear ordinary differential equations (see, for instance, DeLisi and Rescigno [8], Kuznetsov et al. [9], Adams [10], Kirschner and Panetta [11], Abernathy et al. [12], Solis-Perez et al. [13], Alqudah [14], and Simmons et al. [15]), which can only be handled by methods with a region of absolute stability (RAS). Many mathematical models developed describe the interactions among cancer agents (populations), such as tumor cells, immune cells, excess estrogen, healthy cells, and the growing cancer stem cells. Many of them have proved beyond reasonable doubt that a combination of drugs, such as cytokine interleukin-2 (IL-2), and a ketogenic diet can help to fight cancer tumor cells (see Kirschner and Panetta [11] and Oke et al. [16]).

In this study, in a cooperative effort with our academic colleagues, in particular, clinicians and research oncologists, we devoted our attention to mathematical models of tumor growth with the objective of finding a better method of cancer treatment using mathematical concepts and ideas. This is illustrated in the Figure 1 flow diagram, where the various aspects of breast cancer are provided, which give useful information on its development, inversion and metastasis.

2. Mathematically Modeled Problems

From the various studies in the literature, it has been found both in vivo and in vitro that the growth of a tumor cell population is exponential for small quantities of tumor cells, but growth is slowed at large population sizes. The inhibition of growth may be caused by the competition of cells for metabolites or growth factors. In many instances of non-exponential tumor growth, the kinetics are well described by the logistic equation or Gompertz equation. For example, consider a nonlinear rigid system of ordinary differential equations [9] written as follows:

$$\begin{array}{ll}\frac{d{y}_1(t)}{dt}=\sigma +\frac{\rho y_1(t)y_2(t)}{\eta +y_2(t)}-\mu y_1(t)y_2(t)-\delta y_1(t),\hfill & \hfill y_1(0)=D,\\ \frac{d{y}_2(t)}{dt}=\alpha y_2(t)(1-\beta y_2(t))-y_1(t)y_2(t),\hfill & \hfill y_2(0)=D,\end{array}$$

(1)

where $y_1(t)$ are the effector cells, $y_2(t)$ are the tumor cells and [σ, ρ, η, μ, δ, α, β] are the parameters of the system, with values given as [0.1181, 1.131, 20.19, 0.00311, 0.3743, 1.636, 2.0 × 10⁻³], respectively. The practical solution of interest of the system in (1), that is, of the populations denoted by $y_1(t)$, $y_2(t)$ can readily be obtained, which are assumed to be positive. We also assume that the parameters in the system are positive and the Ds are the initial conditions.

The second system of equations to be considered is another nonlinear stiff system [11], which are of the form:

$$\begin{array}{ll}\frac{d{y}_1(t)}{dt}=c{y}_2(t)-{\mu }_2{y}_1(t)+\frac{p_1{y}_1(t)y_3(t)}{g_1+y_3(t)}+s_1,\hfill & \hfill y_1(0)=D,\\ \frac{d{y}_2(t)}{dt}=r_2{y}_2(t)(1-b{y}_2(t))-\frac{a{y}_1(t)y_2(t)}{g_2+y_2(t)},\hfill & \hfill y_2(0)=D,\\ \frac{d{y}_3(t)}{dt}=\frac{p_2{y}_1(t)y_2(t)}{g_3+y_2(t)}-{\mu }_3{y}_3(t)+s_2,\hfill & \hfill y_3(0)=D.\end{array}$$

(2)

where all the parameter values are given in

Table 1. Parameter values.

$0\le c\le 0.05$	${\mu }_2=0.03$	$p_1=0.1245$	$g_1=2\times {10}^7$
$g_2=1\times {10}^5$	$r_2=0.18$	$b=1\times {10}^{-9}$	$a=1$
$g_2=1\times {10}^5$	$p_2=5$	$g_3=1\times {10}^3$

and the Ds are also the initial conditions. The $y_1(t)$ is the effector-cells population, and $s_1$ is a treatment parameter, such as lymphokine-activated killed cell (LAK) therapy or tumor-infiltrating lymphocytes (TIL) therapy. The second equation is the tumor cells denoted by $y_2(t)$ and the third equation is the concentration of the interleukin-cells (IL-2) $y_3(t)$. The $s_2$ is a treatment term which represents an external input of IL-2 into the system.

The third system of stiff equations recently considered by [7,14] to improve cancer treatment using chemotherapy and cancer stem cells is presented in (3). This has four systems of ODEs. The dynamic of the system is located in the nonnegative orthant:

$$\begin{array}{l}{R}_{+,0}^4=\left\{y_1\ge 0,y_2\ge 0,y_3\ge 0,y_4\ge 0\right\}.\hfill \\ \frac{d{y}_1(t)}{dt}={\gamma }_1{y}_1(t)-k_s{y}_1(t)y_4(t),\hfill & \hfill y_1(0)=D,\\ \frac{d{y}_2(t)}{dt}=\alpha -\mu y_2(t)+\frac{p_1{y}_1(t)y_2(t)}{(y_2(t)+1)}-p_2(y_3(t)+y_4(t))y_2(t),\hfill & \hfill y_2(0)=D,\\ \frac{d{y}_3(t)}{dt}=r(1-b{y}_3(t))y_3(t)-(p_3{y}_2(t)+k_T{y}_4(t))y_3(t),\hfill & \hfill y_3(0)=D,\\ \frac{d{y}_4(t)}{dt}=-{\gamma }_2{y}_4(t)+v(t),\hfill & \hfill y_4(0)=D\end{array}$$

(3)

In the model, a competition amongst the cancer populations is formulated, which includes concentration of stem cells, effector cells, tumor cells and chemotherapy drug, respectively. We denote $y_1(t)$ as the stem cell concentration, $y_2(t)$ to represent the effector cells, $y_3(t)$ tumor cell concentration and lastly $y_4(t)$ as the chemotherapy drug concentration. All the constant parameters in the system are described in

Table 2. Parameter values in different models of drug concentration.

Parameter	Interpretation	Value
$y_1(0)$	Initial concentration of stem cells	1
$y_2(0)$	Initial concentration of effector cells	1
$y_3(0)$	Density of free tumor	1
${\gamma }_1$	Decay rate of stem cells concentration	−0.02825
$\alpha$	Rate of product of the effector cells	0.17
$\mu$	Natural death rate of the effector cells	0.03
$b$	Carrying capacity of the tumor cells	10−9
$k_s$	Fractional stem cells killed by chemotherapy	1
$p_1$	Maximum proliferation rate of effector cells	0.1245
$r$	Tumor growth rate	0.18
$p_2$	Decay rate of the effector cells killed by tumor cells and chemotherapy	1
$k_T$	Fractional tumor cells killed by chemotherapy	0.9
$p_3$	Decay rate of tumor cells killed by the effector cells	0.9
${\gamma }_2$	Decay rate of chemotherapy drug	6.4
$y_4(0)$	The time dependent external influx of chemotherapy drug	1

.

The last model was studied by [12,13], where $y_1(t)$, $y_2(t)$, $y_3(t)$, $y_4(t)$ and $y_5(t)$ are, respectively, the cancer stem cells, tumor cells, healthy cells, immune cells and excess estrogen:

$$\begin{array}{ll}\frac{d{y}_1(t)}{dt}=k_1{y}_1(t)\left(1-\frac{y_1(t)}{M_1}\right)-{\gamma }_1{y}_1(t)y_4(t)+\frac{p_1{y}_1(t)y_5(t)}{a_1+y_1(t)},\hfill & \hfill y_1(0)=D\\ \frac{d{y}_2(t)}{dt}=k_2{y}_1(t)\left(\frac{y_1(t)}{M_1}\right)\left(1-\frac{y_2(t)}{M_2}\right)-n_1{y}_2(t)-{\gamma }_2{y}_2(t)y_4(t)\hfill \\ +\frac{p_2{y}_2(t)y_5(t)}{({\alpha }_2+y_2(t))},\hfill & \hfill y_2(0)=D\\ \frac{d{y}_3(t)}{dt}=q{y}_3(t)\left(1-\frac{y_3(t)}{M_3}\right)-\delta y_2(t)y_3(t)-\frac{p_3{y}_3(t)y_5(t)}{a_3+y_3(t)}.\hfill & \hfill y_3(0)=D\\ \frac{d{y}_4(t)}{dt}=s+\frac{\rho y_2(t)y_4(t)}{\omega +y_2(t)}-{\gamma }_3{y}_2(t)y_4(t)-n_2{y}_4(t)-\frac{u{y}_4(t)y_5(t)}{v+y_5(t)},& \hfill y_4(0)=D\\ \frac{d{y}_5(t)}{dt}=r-\left(\mu +\frac{d_1{y}_1(t)}{a_1+y_1(t)}+\frac{d_2{y}_2(t)}{a_2+y_2(t)}+\frac{d_3{y}_3(t)}{a_3+y_3(t)}\right)y_5(t),& \hfill y_5(0)=D\end{array}$$

(4)

where some of the parameters in the system are described in

Table 3. Parameter values in different models of drug concentration.

Parameter	Interpretation	Value
$k_1$	Rate of $y_1$ cells division	0.75 day−1
$k_2$	Rate of $y_2$ cells division	0.514 day−1
$q$	Rate of $y_3$ cells division	0.70 day−1
$M_1$	Carrying capacity of $y_1$ cells	$2.27\times {10}^6$ cells $y_1$
$M_2$	Carrying capacity of $y_2$ cells	$2.27\times {10}^7$ cells $y_2$
$M_3$	Carrying capacity of $y_3$ cells	$2.5\times {10}^7$ cells $y_3$
${\gamma }_1$	Death rate of stem cell due to immune cell response	$3\times {10}^{-7} cel{l}_{y_4}^{-1}$ day−1
${\gamma }_2$	Death rate of tumor cell due to immune cell response	$3\times {10}^{-6} cel{l}_{y_4}^{-1}$ day−1
${\gamma }_3$	Death rate of tumor cell due to immune cell response	$1\times {10}^{-7} cel{l}_{y_2}^{-1}$ day−1
$r$	Infusion of estrogen	200 (pg/mL)−1 day−1
$u$	Immune suppression by estrogen	0.20 day−1
$\nu$	Estrogen threshold	400 (pg/mL)−1
$\omega$	Immune cell threshold	$3\times {10}^5$ cellT
$s$	Source rate of immune cells	$3\times {10}^5$ cellI day−1
$\mu$	Washout rate of estrogen by the body	0.97 day−1
$n_1$	Rate of death of the tumor cells	0.01 day−1
$n_2$	Rate of death of the immune cells	0.29 day−1
$\delta$	Death rate of healthy cells due to competition with tumor cells	$6\times {10}^{-8}$ day−1${\mathrm{cell}}_{T^{-1}}$

, while others are described as follows: $\left\{\rho \right\}$ represents the immune cell response and $\{a_1,a_2,a_3\}$ represent the number of $y_1(t)$, $y_2(t)$, $y_3(t)$ cells. The values $\{p_1,p_2,p_3\}$ help estrogen to proliferate cancer stem cells, tumor cells and the rate at which healthy cells are lost to DNA mutation due to the presence of estrogen. The absorption of estrogen by cancer stem cells, tumor cells and healthy cells are denoted by the values $\{d_1,d_2,d_3\}$, respectively.

Our major worry about all the mathematical models developed and studied is the determination of the solutions. Almost all the methods used for the solution of the models, except in the first model (1) which is experimental or practical in nature, are of lower orders and tend to have undesirable time-step restrictions and poor stability properties within the domain of solution. It is of interest to have method(s) with high accuracy and the ability to incorporate some function evaluations at some off-grid points which are widely used in continuous applications to solve real-life problems.

3. Numerical Computational Algorithms

Ordinary differential equations play an important role in the physical sciences, engineering, biological sciences, technology, etc. Some of these equations have difficulty (stiff) in obtaining their solutions, even numerically. These type of differential equations have arisen in several areas of interest, such as medicine, anthropology, engineering and many other fields. However, only very few of these differential equations can be solved analytically. In this paper, we develop continuous block implicit hybrid methods (CBHMs) based on Chebyshev polynomial nodes that have better regions of absolute stability (RAS) over a wide parameter range. The numerical computation of the study is called continuous block implicit hybrid methods (CBHMs) obtained from the Chebyshev points, with the general form as follows:

$$y(x)={\displaystyle \sum _{j=0}^r{\nu }_j(x)}{y}_{n+j}+{\displaystyle \sum _{j=0}^s{\upsilon }_j(x)}f({\overline{x}}_{n+j},y({\overline{x}}_{n+j})),j=0,1,2,\dots ,s.$$

(5)

where

$$y_{n+j}=y(x_n+jh),$$

(6)

$${y^{\prime }}_{n+j}=f_{n+j}=f({\overline{x}}_n+jh,y({\overline{x}}_n+jh))$$

We denote the number of interpolation points $x_{n+j},j=0,1,2,\dots ,r-1$, by r taken from $(x_n,x_{n+1})$, and denote the distinct collocation points ${\overline{x}}_{n+j},j=0,1,2,\dots ,s-1$ by s, also chosen from the interval $[x_n,x_{n+1}]$. The continuous coefficients ${\nu }_j(x)$ and ${\upsilon }_j(x)$ in (5) are to be determined. They are of degree $p-1$$(p=r+s)$ given by:

$${\nu }_j(x)={\displaystyle \sum _{i=0}^{p-1}{\nu }_{j,i+1}}{x}^i\mathrm{and}h{\upsilon }_j(x)=h{\displaystyle \sum _{i=0}^{p-1}{\upsilon }_{j,i+1}{x}^i}.$$

(7)

The coefficients ${\nu }_{j,i+1}$ and ${\upsilon }_{j,i+1}$ in (7) are to be determined explicitly.

Theorem 3.1.

Collocation at the Gaussian points together with the two end points lead to order one less than the superconvergence scheme, while collocation at only one additional end point to the Gaussian points resulted in order two less than the upgraded scheme.

Proof .

(see [17,18,19]). □

Case I:

We shall derive one-step collocation methods as continuous block implicit hybrid formulae of orders 3 and 4. We shall consider the zeros of the Chebychev polynomial of degree 3, which can be normalized to give the general finite range in $a\le x\le b$ as:

$$T_n^{*}(x)=\cos (n\arccos (t)),t=\frac{2x-(b+a)}{b-a}$$

Hence, with n = 3, we obtain the zeros of the Chebyshev polynomial as follows:

$$\begin{array}{l}{\overline{x}}_0=x_{n+u},u=(2-\sqrt{3})/4\\ {\overline{x}}_1=x_{n+w},w=\frac{1}{2}\\ {\overline{x}}_2=x_{n+v},v=(2+\sqrt{3})/4\end{array}\}$$

In this family $r=1,s=3,\xi =(x-x_n)$, the continuous scheme is given as:

$$y(x)={\nu }_0(x)y_n+h[{\upsilon }_1(x)f_{n+u}+{\upsilon }_2(x)f_{n+w}+{\upsilon }_3(x)f_{n+v}].$$

(8)

where

$$\begin{array}{l}{\nu }_0(x)=1\hfill \\ {\upsilon }_1(x)=\left[\frac{2{(x-x_n)}^3-3(v+w)h{(x-x_n)}^2+6vw{h}^2(x-x_n)}{6(v-u)(w-u)h^2}\right],\hfill \\ {\upsilon }_2(x)=\left[\frac{2{(x-x_n)}^3-3(v+u)h{(x-x_n)}^2+6vu{h}^2(x-x_n)}{6(w-u)(w-v)h^2}\right],\hfill \\ {\upsilon }_3(x)=\left[\frac{-2{(x-x_n)}^3+3(u+w)h{(x-x_n)}^2+6uw{h}^2(x-x_n)}{6(v-u)(w-v)h^2}\right]\hfill \end{array}$$

Evaluating (8), we obtain the continuous block implicit hybrid method of order 4 with only 3 stages written in the form:

$$y_{n+1}=y_n+\frac{h}{9}\left[2f_{n+u}+5f_{n+w}+2f_{n+v}\right],\mathrm{order}\mathrm{p}=4,C_5=-1.3020\times {10}^{-4}$$

(9)

$$\begin{gathered} y_{n+u}=y_n+\frac{h}{144}\left[\left(16-3\sqrt{3}\right)f_{n+u}+\left(40-24\sqrt{3}\right)f_{n+w}+\left(16-9\sqrt{3}\right)f_{n+v}\right], \\ \mathrm{order}\mathrm{p}=3,C_4=-1.6276\times {10}^{-4} \end{gathered}$$

$$\begin{gathered} y_{n+w}=y_n+\frac{h}{36}\left[\left(4+3\sqrt{3}\right)f_{n+u}+10f_{n+w}+\left(4-3\sqrt{3}\right)f_{n+v}\right], \\ \mathrm{order}\mathrm{p}=3,C_4=1.3020\times {10}^{-4} \end{gathered}$$

$$\begin{gathered} y_{n+v}=y_n+\frac{h}{144}\left[\left(16+9\sqrt{3}\right)f_{n+u}+\left(40+24\sqrt{3}\right)f_{n+w}+\left(16+3\sqrt{3}\right)f_{n+v}\right], \\ \mathrm{order}\mathrm{p}=3,C_4=-1.3020\times {10}^{-4} \end{gathered}$$

Note that the matrix A obtained from the derived method is not a lower triangular matrix hence it is a fully dense matrix, that is, a fully implicit method. In addition, it has higher order than the explicit case, which has fewer stages too. Using the maple package, yields the stability polynomial of the block implicit hybrid method as

$$p(\eta ,z)=\frac{\eta (z^3+\eta z^3+18z^2-18\eta z^2+96\eta z-192\eta +96z+192)}{z^3-18z^2+96z-192}$$

The RAS of the block implicit hybrid method is plotted using the MATLAB package as shown in Figure 2a.

Case II:

Collocating at the three interior points above and one additional end point we obtain the continuous scheme:

$$y(x)={\nu }_0(x)y_n+h[{\upsilon }_1(x)f_{n+u}+{\upsilon }_2(x)f_{n+w}+{\upsilon }_3(x)f_{n+v}+{\upsilon }_4(x)f_{n+1}]$$

(10)

where

$${\nu }_0(x)=1,$$

$${\upsilon }_1(x)=\left[\frac{-{(x-x_n)}^4+4(v+u+w)h{(x-x_n)}^3-6(vu+vw+uw)h^2{(x-x_n)}^2+12vuw{h}^3(x-x_n)}{12(v-1)(u-1)(w-1)h^3}\right],$$

$${\upsilon }_2(x)=\left[\frac{3{(x-x_n)}^4-4(v+w+1)h{(x-x_n)}^3+6(vw+v+w)h^2{(x-x_n)}^2-12vw{h}^3(x-x_n)}{12(v-u)(w-u)(u-1)h^3}\right],$$

$${\upsilon }_3(x)=\left[\frac{-3{(x-x_n)}^4+4(v+u+1)h{(x-x_n)}^3-6(vu+v+u)h^2{(x-x_n)}^2+12vu{h}^3(x-x_n)}{12(v-w)(w-u)(w-1)h^3}\right],$$

$${\upsilon }_3(x)=\left[\frac{-3{(x-x_n)}^4+4(v+u+1)h{(x-x_n)}^3-6(vu+v+u)h^2{(x-x_n)}^2+12vu{h}^3(x-x_n)}{12(v-w)(w-u)(w-1)h^3}\right],$$

$${\upsilon }_4(x)=\left[\frac{3{(x-x_n)}^4-4(u+w+1)h{(x-x_n)}^3+6(uw+u+w)h^2{(x-x_n)}^2-12uw{h}^3(x-x_n)}{12(v-u)(v-w)(v-1)h^3}\right].$$

Evaluating (10) at the off-step points $x_{n+u}$, $x_{n+w}$, $x_{n+v}$ and at the end point $x_{n+1}$, the following continuous block implicit hybrid method is obtained as follows:

$$y_{n+1}=y_n+\frac{h}{9}\left[2f_{n+u}+5f_{n+w}+2f_{n+v}\right],\mathrm{order}\mathrm{p}=4,C_5=-1.3020\times {10}^{-4}$$

(11)

$$y_{n+u}=y_n+\frac{h}{288}\left[\left(38-9\sqrt{3}\right)f_{n+u}+\left(77-48\sqrt{3}\right)f_{n+w}+\left(38-15\sqrt{3}\right)f_{n+v}-9f_{n+1}\right]\mathrm{order}\mathrm{p}=4,C_5=3.9813\times {10}^{-5}$$

$$y_{n+w}=y_n+\frac{h}{36}\left[\left(2-6\sqrt{3}\right)f_{n+u}+13f_{n+w}-\left(2+6\sqrt{3}\right)f_{n+v}+9f_{n+1}\right]\mathrm{order}\mathrm{p}=4,C_5=-2.2786\times {10}^{-4}$$

$$y_{n+v}=y_n+\frac{h}{288}\left[\left(38+15\sqrt{3}\right)f_{n+u}+\left(77+48\sqrt{3}\right)f_{n+w}+\left(38+9\sqrt{3}\right)f_{n+v}-9f_{n+1}\right]\mathrm{order}\mathrm{p}=4,C_5=-3.9813\times {10}^{-5}$$

The continuous block implicit hybrid method of uniform order p = 4 everywhere in the interval of the solution is obtained for the solution of dynamics breast cancer models. Using the MAPLE package yields the stability polynomial of the method as:

$$p(\eta ,z)=\frac{\eta (768\eta +132\eta z^2-19\eta z^3+\eta z^4-480\eta z-768-36z^2-z^4-288z)}{z^4-19z^3+132z^2-480z+768}$$

The RAS of the continuous block implicit hybrid method was plotted using MATLAB package as shown in Figure 2b.

Regions of absolute stability of the continuous block implicit hybrid methods which are A-stable and symmetric about the x-axis, since the enclosed Figures of the regions contain the complex plane outside the domain.

From the two block implicit hybrid methods obtained, it is clear that the Chebyshev polynomial roots are a good approximant of ordinary differential equations. Evidently, Ibiejugba and Onumanyi [20] pointed out that the selection of the elements using the zeros of some appropriate Chebyshev polynomials yield some interesting properties, such as the satisfaction of the usual local support. According to the authors, the idea of selecting the elements based on the zeros of Chebyshev polynomials as described in their paper is in line with the principle of the “orthogonal collocation method” or “method of selected points”.

4. Results and Discussion

Many mathematical models of interactions between the immune system, chemotherapies (anticancer agents) and the growing tumor have been developed in the recent literature. Immunotherapies which are responsible for stimulating antigen-sensitized cells are becoming a very significant component for treating cancer. This can easily be done by increasing the number of effector cells, either by stem cell transformation or by using the combination of cytokine interleukin-2 together with adoptive cellular immunotherapy, which will have the energy to fight the tumor cells by reducing their number and strength. This has been demonstrated in the laboratory and also in clinical experiments (see, for example, Farrar et al. [21], O’Byrne et al. [22], and de Pillis et al. [23]). Apparently, immune system cells usually attacked and killed cancer cells and therefore kept cancer incidence very low (see [24,25]).

From the Figures showing the graphical representation of the results in this section, transients near the steady state show decaying oscillations in the plots. The graphical representations of the results, for example, Figure 3 and Figure 4, are of interest due to the fact that they portrait stable limit cycles, showing that the tumor and the immune system are subject to oscillations. In some cases, cyclic fluctuations in the number of leukocytes were observed. We noticed from the given graphical representation of the results that the system can be stable only if the given immune resistance is greater than the tumor growth. The graphs of Figure 3 for $y_1(0)=5$ and $y_2(0)=50$ in the neighborhood are clearly depicted (presented). The plots of Figure 3c–f for one and two months are in good agreement with certain human leukemia found in the literature (see, for example, [9,11,23]). In the graphs of Figure 3g–j, decaying oscillations to a dormant tumor state are also clearly depicted. However, in the case of higher doses for a longer period of time, the system evolves to an oscillatory state, as shown or depicted in Figure 3g–j. However, in the graphs of Figure 3a,b, as we can see, the effect is small because of the number of days considered. From Figure 3a,b, no oscillation is experienced in the graphs.

The parameter estimates in Table 1 are biologically reasonable because they characterize tumor growth in the absence of an immune response and the data determine the slope and asymptotes of the tumor growth curves. From the efficiency curves generated using Equation (2) or Model II, we can make a number of biological predictions about the interactions of the immune system with the growth tumor, as we did in Equation (1) or Model I. Compared with some published results in scientific articles, cyclic fluctuations in the number of leukocytes were found in several cases (see, for example [8,10]).

Practically, Equation (2) described the interaction amongst the effector cells $y_1(t)$, tumor cells $y_2(t)$ and the cytokine interleukin-2(IL-2) $y_3(t)$. As shown in the graphical plots of Figure 4, the initial points, even before introducing the treatment, were almost the same: these are $y_1(0)=1\times {10}^7$, $y_2(0)=3\times {10}^5$, $y_3(0)=0.25\times {10}^{10}$ or $y_1(0)$ = $y_2(0)$ = $y_3(0)$ = 1 or $y_1(0)$ = $y_2(0)$ = $y_3(0)$ = 0. The idea, then, is to improve the natural immune system by using (IL-2), which is usually used together with adoptive cellular immunotherapy (ACI). Cytokines, which are stimulating protein hormones, enhance both natural and specific immunity (effector cells). Here, we indicate the treatment $s_1$ in the equation as the use of adoptive cellular immunity (ACI). Administering a large amount of the interleukin-2(IL-2), which is denoted by $s_2$ in the equation, yields an interesting result. It completely cleared the tumor cells but also has the side effect of capillary leakage syndrome (or vascular leakage syndrome) according to some eminent authors (see, for example, [26,27]). As observed in Figure 3i,j, which are very similar, the difference in oscillations is connected or linked with the values of the initial dose.

As shown in the graphs of Figure 5, the tumor cell population reduced by using only a small constant initial dose of chemotherapy for a short period of time. The idea of using the continuous block implicit hybrid methods is to accurately determine solutions of the modeled differential equations so that even a minimum initial dose of treatment, that is, the last term of Equation (3), $V_M$, will be able to reduce (Figure 5a,b,g,h) or eliminate (Figure 5c,d,e,f) the tumor population described by Equation (3), that is, the term $V(t)$ in the equation. In the recent literature, some authors, for example, de Pillis et al. [28], suggested the use of a high dose for a few days to eliminate the tumor. Though, if one uses high dose for few days, we observed that the case of increase in immune cells will cause a serious side effect of the treatment, as discussed in the case of the graphs of Figure 4 above. As can be seen in Figure 5, our solution graphs differ from that of Figure 4 in terms of the initial dose and the period of applications, which guaranteed complete tumor clearance without much side effect.

In Figure 6, we solved and analyzed a model of cancer disease which considered five populations involving $y_1(t)$, $y_2(t)$, $y_3(t)$, $y_4(t)$ and $y_5(t)$, which are cancer stem cells, tumor cells, healthy cells, immune cells and excess estrogen, respectively [6,13]. As seen in Figure 6a,b, which depicts the behavior of the graphs of the sub-model without estrogen over time, with cancer stem cells, tumor cells, healthy cells and immune cells all beginning at their carrying capacities, cancer still persists at a very high level. If observed carefully, the natural immune cell is not sufficient to even reduce the cancer, much less completely eradicate the tumor, as observed in [6,13,29].

Further, from the graphs of the Figure, it is clearly shown that the presence of excess estrogen, tumor cells increase, but immune cells and healthy cells decrease. In such a situation the normal cells are very much being affected, which makes the system unstable. As can be seen in Figure 6c–h, with the increase in estrogen levels, the immune levels are completely reduced, thereby weakening the body’s immune system [13]. This suggests that, in the presence of excess estrogen, all breast tissue eventually becomes infected with tumor cells. Here, we advise that taking extra levels of estrogen, either as hormonal birth control or to enhance beauty, ends up practically increasing the risk of developing breast cancer and should therefore be avoided.

5. Conclusions

Biological systems are very close to clinical medicine. This is because experiments are performed to test samples. This is where mathematical models are significant, especially in biomedical science where some complex biological models are tested. In this study, we have illustrated how mathematical techniques could be used in biology and medicine, which are difficult to achieve through experimental work alone. These results could prove important and useful in our efforts to effectively deliver a variety of novel and exciting immunotherapies, particularly the cytokine interleukin-2, adoptive cellular immunotherapy, etc. Overall, the study shows that early treatment is the best way to eliminate the tumor. Hence, our results presented should be useful, especially for choosing the optimal treatment strategy which does not cause great negative impact on the cancer patient as demonstrated above graphically. It is interesting to note that in this article we provide a graphical representation of the results obtained from the solutions of modeled equations where the singular control is optimal with respect to the solutions. Therefore, the simulated results presented here show that the CBHMs are good candidates for generating results that represent various asymptomatic behaviors that track the real data more closely than most of the methods used in the literature. We hope that the methods used here can be applied to other biologically complex models.