Mathematical Model to Understand the Dynamics of Cancer, Prevention Diagnosis and Therapy

Abstract

In the present study, we formulate a mathematical model to understand breast cancer in the population of Saudi Arabia. We consider a mathematical model and study its mathematical results. We show that the breast cancer model possesses a unique system of solutions. The stability results are shown for the model. We consider the reported cases in Saudi Arabia for the period 2004–2016. The data are given for the female population in Saudi Arabia that is suffering from breast cancer. The data are used to obtain the values of the parameters, and then we predict the long-term behavior with the obtained numerical results. The numerical results are obtained using the proposed parameterized approach. We present graphical results for the breast cancer model under effective parameters such as ${\tau }_1$, ${\tau }_2$, and ${\tau }_3$ that cause decreasing future cases in the population of stages 3 and 4, and the disease-free condition. Chemotherapy generally increases the risk of cardiotoxicity, and, hence, our model result shows this fact. The combination of chemotherapy stages 3 and 4 and the parameters ${\tau }_1$ and ${\tau }_2$ together at a low-level rate and also treating the patients before the chemotherapy will decrease the population of cardiotoxicity. The findings of this study are intended to reduce the number of cardiotoxic patients and raise the number of patients who recover following chemotherapy, which will aid in public health decision making.

1. Introduction

Breast cancer is a condition in which the breast cells multiply out of control [1]. According to data from “The Global Burden of Disease Cancer Collaboration”, breast cancer has the greatest incidence of all types of cancer [2,3]. Breast cancer alters the structure of cells and breast tissue, causing it to proliferate uncontrollably [4]. Breast cancer is a possibility for every woman around the globe. Breast cancer was ranked as the second most common cancer by the World Health Organization (WHO) in 2004. According to a WHO poll, 8 to 9 percent of women worldwide develop breast cancer. Experts in global health are still unsure of the exact etiology of breast cancer [5]. According to [5] and the study given in [6], 685,000 people worldwide died in 2020 as a result of breast cancer, which affected 2.3 million women. By the end of 2020, the most common disease in the globe was breast cancer, which had been diagnosed in 7.8 million women in the previous five years. Worldwide, breast cancer causes more lost disability-adjusted life years (DALYs) for women than any other kind of cancer. Every country in the world experiences breast cancer in women after puberty at any age, albeit the incidence rates rise as people age [6].

Over time, there has been a significant advancement in the understanding of the diagnosis, development, progression, metastasis, and treatment of breast cancer. Despite remarkable advances in our knowledge of the condition over the last 50 years, it continues to be a major public health burden and a major worry on a global scale. The invasive type of cancer that affects women the most frequently in the globe is breast cancer [7]. In Saudi Arabia, breast cancer was the type of cancer that was most frequently diagnosed in 2018, and it was also the second-leading cause of death after leukemia [8]. Although Saudi Arabia has a far lower incidence of breast cancer than many Western countries, there is mounting evidence that the figures are rising quickly. Exploring breast cancer characteristics, trends, age, and regional disparities is one of the top-priority study subjects in a developing country such as Saudi Arabia, where more than 70% of women are under the age of 39. From 2001 to 2008, the most recent available full epidemiological statistics on the burden of breast cancer were supplied [9].

Mathematical models may predict how infectious illnesses will progress, indicating the expected result of an epidemic and assisting in the development of public health and health treatments, [10,11]. The voyage of developing mathematical models for cancer began in 1954 [12] to explain the disease, and, subsequently, the researcher analyzed several elements of cancer and tumor progression [13,14]. In [15], the authors constructed a mathematical model of chemotherapy treatment for cancer; the authors illustrated the treatment approach for tumor cancer. Cancer dynamics depict the interplay of tumor cell energy and density, as well as the influence of treatment medications. Recently, a mathematical model was established for low-dose chemotherapy with few parameters [16]; it analyzed angiogenic signals between vasculature and malignancies. Jordao and Tavares established a compartmental model in [17]; they evaluated malignant and healthy cells to analyze the proposed cancer model. S. Khajanchi and Nieto [18] studied the influence of time delay on the dynamics of the tumor system. Mahlbacher et al. [19] established another essential model to comprehend the connections between the immunityand tumor and forecast better cancer recommendations. Many mathematical models have been constructed and formulated in the literature to explore, understand, and depict cancer transmission processes [20,21,22].

In the present work, we investigate breast cancer’s dynamics under real data and consider a novel numerical technique to investigate its numerical results. We consider the modeling of breast cancer growth in detail and then study the essential characteristics of the model. Later, we investigate the equilibrium points and their stability. The equilibria are carried out and their analysis is given. Some interesting results regarding data versus model fitting are given. We predict the long-term behavior of the disease and also provide the importance of the parameters that can possibly reduce the future infected cases of breast cancer in Saudi Arabia.

The work in this manuscript is divided section-wise as follows: In Section 2, we explain, in detail, the modeling of breast cancer. The analysis of the equilibrium point is discussed in Section 3. Section 4 explains the new numerical scheme for the breast cancer model. Numerical simulations and a discussion of the results are shown in Section 5 whereas Section 6 summarizes the results.

2. Model Formulation

This section formulates the model for breast cancer considering five epidemiological compartments. The compartments are classified based on breast cancer epidemiology, which is defined as individuals admitted to the hospital due to the effect of cancer, denoted by X. According to the medical record, the admitted patients are classified into subpopulations of stage 1 to 4. We assume that all individuals with cancer receive chemotherapy treatment. With the severity of cancer in patients over time, each population will experience changes in their respective population. Among these individuals, some have serious conditions, while some recover, and those people experience cardiotoxicity during the chemotherapy process. We briefly explain below the process of modeling cancer patients initially admitted to the hospital and their progress to stages 1 to 4. The individuals initially admitted to the hospital are identified as cancer patients shown by $X\left(t\right)$ who have stage 1, 2X, and 2B. The population of class A generated by ${\mathsf{\Lambda }}_1$ are identified to be cancer patients with stages 1 and 2. The individuals in X with chemotherapy recover with the rate ${\varphi }_1$ (join class R), while those who have the worst conditions are given the rate ${\beta }_2$ (join class B). This discussion leads to the following ordinary differential equations:

$$\begin{array}{cc}\hfill \frac{dX}{dt}& ={\mathsf{\Lambda }}_1-({\varphi }_1+{\beta }_2)X.\end{array}$$

(1)

The subpopulation B contains all those individuals who are identified to be patients with stages 3X and 3B. The population of group B is generated through the rate ${\mathsf{\Lambda }}_2$, which shows the individuals that are suffering from stage 3 and those who have the worst conditions in class X join class B through the rate ${\varphi }_1$. Individuals in class B die due to cancer at the rate ${\mu }_1$. Individuals who are cured after chemotherapy join class R with the rate ${\varphi }_2$, and those have worse conditions join class C with the rate ${\beta }_1$. The individuals in class B with more intensive chemotherapy than class X cause cardiotoxicity with the rate ${\tau }_1$ and join class E. The patients in class R are reinfected with cancer at a rate ${\psi }_1$, so the following evolutionary equations are formulated:

$$\begin{array}{cc}\hfill \frac{dB}{dt}& ={\mathsf{\Lambda }}_2+{\beta }_{2}X+{\psi }_{1}R-({\varphi }_2+{\beta }_1+{\tau }_1+{\mu }_1)B.\end{array}$$

(2)

Subjects who are identified after testing to be at stage 4 are shown by class C. Individuals who are tested and identified as cancer patients enter the population C at a rate ${\mathsf{\Lambda }}_3$. Those from class B with worse conditions are at the rate ${\beta }_1$, and ${\varphi }_3$ is the cure rate due to chemotherapy that joins class R. The rate ${\tau }_2$ towards cardiotoxicity is assumed to be of great value because the patient is experiencing very intense side effects of chemotherapy, the death rate in class C is given by ${\mu }_2$. The individuals in class R again join class C with the rate ${\psi }_2$. The above discussion is shown by the following system:

$$\begin{array}{cc}\hfill \frac{dC}{dt}& ={\mathsf{\Lambda }}_3+{\beta }_{1}B+{\psi }_{2}R-({\varphi }_3+{\tau }_2+{\mu }_2)C.\end{array}$$

(3)

The class R shows individuals with disease-free conditions after using chemotherapy. The individuals in classes X, B, and C showing disease-free conditions after chemotherapy are joining the class R with rates ${\varphi }_1$, ${\varphi }_2$, and ${\varphi }_3$ respectively. The rates ${\psi }_1$ and ${\psi }_2$ show the relapse back to stages 3 and 4, respectively. An extended period of subpopulation R can also experience direct cardiotoxicity with rate ${\tau }_3$. This is shown by the following mathematical expression:

$$\begin{array}{cc}\hfill \frac{dR}{dt}& ={\varphi }_{1}X+{\varphi }_{2}B+{\varphi }_{3}C-({\psi }_1+{\psi }_2+{\tau }_3)R.\end{array}$$

(4)

Subpopulation E represents cancer patients who have cardiotoxicity. This population is generated by individuals from classes R, C, and B with the rates ${\tau }_3$, ${\tau }_2$, and ${\tau }_1$, respectively, while ${\mu }_3$ is the death rate in class E. This is shown by the mathematical expression given by

$$\begin{array}{cc}\hfill \frac{dE}{dt}& ={\tau }_{3}R+{\tau }_{2}C+{\tau }_{1}B-{\mu }_{3}E.\end{array}$$

(5)

The equations given in (1)–(5) can be written as a joint system, as follows:

$$\begin{array}{cc}\hfill \frac{dX}{dt}& ={\mathsf{\Lambda }}_1-({\varphi }_1+{\beta }_2)X,\hfill \\ \hfill \frac{dB}{dt}& ={\mathsf{\Lambda }}_2+{\beta }_{2}X+{\psi }_{1}R-({\varphi }_2+{\beta }_1+{\tau }_1+{\mu }_1)B,\hfill \\ \hfill \frac{dC}{dt}& ={\mathsf{\Lambda }}_3+{\beta }_{1}B+{\psi }_{2}R-({\varphi }_3+{\tau }_2+{\mu }_2)C,\hfill \\ \hfill \frac{dR}{dt}& ={\varphi }_{1}X+{\varphi }_{2}B+{\varphi }_{3}C-({\psi }_1+{\psi }_2+{\tau }_3)R,\hfill \\ \hfill \frac{dE}{dt}& ={\tau }_{3}R+{\tau }_{2}C+{\tau }_{1}B-{\mu }_{3}E,\hfill \end{array}$$

(6)

with the non-negative initial conditions

$$\begin{array}{c}\hfill X\left(0\right)\ge X_0,B\left(0\right)\ge B_0,C\left(0\right)\ge C_0,R\left(0\right)\ge R_0,E\left(0\right)\ge E_0.\end{array}$$

The details of the parameters are shown in

Table 1. Parameter details.

Symbol	Definitions
${\mathsf{\Lambda }}_1$	People identified as cancer patients at stage 1 and 2
${\mathsf{\Lambda }}_2$	People suffering from cancer at stage 3
${\mathsf{\Lambda }}_3$	People suffering from cancer at stage 4
${\varphi }_1$	Recovery due to chemotherapy at stages 1 and 2
${\varphi }_2$	Recovery due to chemotherapy at stage 3
${\varphi }_3$	Recovery due to chemotherapy at stage 4
${\beta }_1$	People with worse conditions join stage 4 population
${\beta }_2$	People with worse conditions join class B
${\tau }_1$	People with intensive chemotherapy that causes cardiotoxicity
${\tau }_2$	People at stage 4 cancer chemotherapy who experience cardiotoxicity
${\tau }_3$	People with intensive chemotherapy that causes cardiotoxicity at disease-free stage
${\mu }_1$	Death due to cancer at stage 3
${\mu }_2$	Death due to cancer at stage 4
${\mu }_3$	Death rate of cardiotoxic people
${\psi }_1$	People relapse back to stage 3
${\psi }_2$	People relapse back to stage 4

.

Existence and Uniqueness

We investigate the existence and uniqueness of the model (6). To perform this, we present the model (6) in the form given by

$$\begin{array}{ccc}\hfill X^{\prime }\left(t\right)& =& f_1(t,X,B,C,R,E),\hfill \end{array}$$

$$\begin{array}{ccc}\hfill B^{\prime }\left(t\right)& =& f_2(t,X,B,C,R,E),\hfill \end{array}$$

$$\begin{array}{ccc}\hfill C^{\prime }\left(t\right)& =& f_3(t,X,B,C,R,E),\hfill \end{array}$$

$$\begin{array}{ccc}\hfill R^{\prime }\left(t\right)& =& f_4(t,X,B,C,R,E),\hfill \end{array}$$

$$\begin{array}{ccc}\hfill E^{\prime }\left(t\right)& =& f_5(t,X,B,C,R,E).\hfill \end{array}$$

(7)

We define the norm given by

$$\begin{array}{c}\hfill {\parallel X\parallel }_{\infty }=\underset{t\in U}{sup}\left|X\left(t\right)\right|,\end{array}$$

(8)

where $U\in [0,T]$. Here, it is assumed that $X\left(t\right)$, $B\left(t\right)$, $C\left(t\right)$, $R\left(t\right)$, and $E\left(t\right)$ are bounded in $[0,T]$ for any $t\in [0,T]$, and there exist $S_1$, .... $S_5$ such that $\parallel X\parallel \infty <S_1$, $\parallel B\parallel \infty <S_2$, $\parallel C\parallel \infty <S_3$, $\parallel R\parallel \infty <S_4$ and $\parallel E\parallel \infty <S_5$. We have to prove $f_1,...,f_5$ are bounded.

$$\begin{array}{c}\hfill |f_1(t,X,B,C,R,E)|=|{\mathsf{\Lambda }}_1-({\varphi }_1+{\beta }_2)X|,\end{array}$$

$$\begin{array}{c}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\le {\mathsf{\Lambda }}_1+({\varphi }_1+{\beta }_2)\left|X\right|,\hfill \end{array}$$

$$\begin{array}{c}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\le {\mathsf{\Lambda }}_1+({\varphi }_1+{\beta }_2)\underset{L}{sup}\left|X\right|,\hfill \end{array}$$

$$\begin{array}{c}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\le {\mathsf{\Lambda }}_1+({\varphi }_1+{\beta }_2){\parallel X\parallel }_{\infty },\hfill \end{array}$$

$$\begin{array}{c}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\le {\mathsf{\Lambda }}_1+({\varphi }_1+{\beta }_2)S_1<\infty ,\hfill \end{array}$$

(9)

where $L=t\in [0,T]$. We follow the above method and can easily show the rest:

$$\begin{array}{ccc}\hfill |f_2(t,X,B,C,R,E)|& =& |{\mathsf{\Lambda }}_2+{\beta }_{2}X+{\psi }_{1}R-{\mathsf{\Phi }}_{2}B|,\hfill \end{array}$$

$$\begin{array}{c}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\le {\mathsf{\Lambda }}_2+{\beta }_2{S}_1+{\psi }_1{S}_4+{\mathsf{\Phi }}_2{S}_2<\infty .\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill |f_3(t,X,B,C,R,E)|& =& |{\mathsf{\Lambda }}_3+{\beta }_{1}B+{\psi }_{2}R-{\mathsf{\Phi }}_{3}C|,\hfill \end{array}$$

$$\begin{array}{c}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\le {\mathsf{\Lambda }}_3+{\beta }_1{S}_2+{\psi }_2{S}_4+{\mathsf{\Phi }}_3{S}_3<\infty .\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill |f_4(t,X,B,C,R,E)|& =& |{\varphi }_{1}X+{\varphi }_{2}B+{\varphi }_{3}C-({\psi }_1+{\psi }_2+{\tau }_3)R|,\hfill \end{array}$$

$$\begin{array}{c}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\le {\varphi }_1{S}_1+{\varphi }_2{S}_2+{\varphi }_3{S}_3+{\mathsf{\Phi }}_4{S}_4<\infty .\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill |f_5(t,X,B,C,R,E)|& =& |{\tau }_{3}R+{\tau }_{2}C+{\tau }_{1}B-{\mu }_{3}E|,\hfill \end{array}$$

$$\begin{array}{c}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\le {\tau }_3{S}_4+{\tau }_2{S}_3+{\tau }_1{S}_2+{\mu }_3{S}_5<\infty ,\hfill \end{array}$$

(13)

where ${\mathsf{\Phi }}_2=({\varphi }_2+{\beta }_1+{\tau }_1+{\mu }_1)$, ${\mathsf{\Phi }}_3=({\varphi }_3+{\tau }_2+{\mu }_2)$, ${\mathsf{\Phi }}_4=({\psi }_1+{\psi }_2+{\tau }_3)$. Thus, $X\left(t\right)$, $B\left(t\right)$, $C\left(t\right)$, $R\left(t\right)$ and $E\left(t\right)$ are bounded, and then there exist $J_1,...,J_5$ such that

$$\begin{array}{ccc}\hfill \underset{L}{sup}\left|f_1(t,X,B,C,R,E)\right|& <& J_1,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill \underset{L}{sup}\left|f_2(t,X,B,C,R,E)\right|& <& J_2,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill \underset{L}{sup}\left|f_3(t,X,B,C,R,E)\right|& <& J_3,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill \underset{L}{sup}\left|f_3(t,X,B,C,R,E)\right|& <& J_4,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill \underset{L}{sup}\left|f_5(t,X,B,C,R,E)\right|& <& J_5.\hfill \end{array}$$

(14)

Further, we have to show

$$\begin{array}{ccc}\hfill |f_1(t,X_1,B,C,R,E)-f_1(t,X_2,B,C,R,E)|& =& |{\mathsf{\Lambda }}_1-({\varphi }_1+{\beta }_2)X_1-{\mathsf{\Lambda }}_1+({\varphi }_1+{\beta }_2)X|,\hfill \end{array}$$

$$\begin{array}{ccc}& \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}<& ({\varphi }_1+{\beta }_2)|X_1-X_2|,\hfill \end{array}$$

$$\begin{array}{ccc}& \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}<& {\overline{u}}_1|X_1-X_2|,\hfill \end{array}$$

(15)

where ${\overline{u}}_1=({\varphi }_1+{\beta }_2)$. With the same approach, we show

$$\begin{array}{ccc}\hfill |f_2(t,X,B_1,C,R,E)-f_2(t,X,B_2,C,R,E)|& <& {\overline{u}}_2|B_1-B_2|,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill |f_3(t,X,B,C_1,R,E)-f_3(t,X,B,C_2,R,E)|& <& {\overline{u}}_3|C_1-C_2|,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill |f_4(t,X,B,C,R_1,E)-f_4(t,X,B,C,R_2,E)|& <& {\overline{u}}_4|R_1-R_2|,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill |f_5(t,X,B,C,R,E_1)-f_5(t,X,B,C,R,E_2)|& <& {\overline{u}}_5|E_1-E_2|,\hfill \end{array}$$

(16)

where ${\overline{u}}_2=({\varphi }_2+{\beta }_1+{\tau }_1+{\mu }_1)$, ${\overline{u}}_3=({\varphi }_3+{\tau }_2+{\mu }_2)$, ${\overline{u}}_4=({\psi }_1+{\psi }_2+{\tau }_3)$, ${\overline{u}}_5={\mu }_3$.

This uniqueness and existence of the solution of the model (6) indicate that the model has a unique system of solutions.

3. Equilibrium Points and Their Analysis

We denote the equilibrium point of the model by $P_1$, and it is obtained as follows:

$$\begin{array}{c}\hfill P_1=(X^{*},B^{*},C^{*},R^{*},E^{*}),\end{array}$$

(17)

where,

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{X}^{*}=\frac{{\mathsf{\Lambda }}_1}{{\mathsf{\Phi }}_1},\hfill \\ \\ B^{*}=\frac{K_1}{K_2},\hfill \\ \\ C^{*}=\frac{K_3}{K_2},\hfill \\ \\ R^{*}=\frac{K_4}{K_2},\hfill \\ \\ E^{*}=\frac{K_5}{\mu K_2},\hfill \end{array}\right.\end{array}$$

where ${\mathsf{\Phi }}_1=({\varphi }_1+{\beta }_2)$, ${\mathsf{\Phi }}_2=({\varphi }_2+{\beta }_1+{\tau }_1+{\mu }_1)$, ${\mathsf{\Phi }}_3=({\varphi }_3+{\tau }_2+{\mu }_2)$, ${\mathsf{\Phi }}_4=({\psi }_1+{\psi }_2+{\tau }_3)$, and

$$\begin{array}{ccc}\hfill K_1& =& {\mathsf{\Phi }}_3\left({\mathsf{\Lambda }}_1{\psi }_1{\varphi }_1+{\mathsf{\Phi }}_4\left({\beta }_2{\mathsf{\Lambda }}_1+{\mathsf{\Lambda }}_2{\mathsf{\Phi }}_1\right)\right)-{\psi }_2\left({\mathsf{\Phi }}_3+{\varphi }_3\right)\left({\beta }_2{\mathsf{\Lambda }}_1+{\mathsf{\Lambda }}_2{\mathsf{\Phi }}_1\right)+{\mathsf{\Lambda }}_3{\mathsf{\Phi }}_1{\psi }_1{\varphi }_3,\hfill \\ \\ \hfill K_2& =& {\mathsf{\Phi }}_1{\mathsf{\Phi }}_2{\mathsf{\Phi }}_3{\mathsf{\Phi }}_4\left(1-{\mathcal{R}}_0\right),\hfill \\ \\ \hfill K_3& =& {\psi }_1\left({\beta }_1{\mathsf{\Lambda }}_1{\varphi }_1-{\mathsf{\Lambda }}_3{\mathsf{\Phi }}_1{\varphi }_2\right)+{\psi }_2\left({\mathsf{\Lambda }}_1\left({\beta }_2{\varphi }_2+{\mathsf{\Phi }}_2{\varphi }_1\right)-{\beta }_1\left({\beta }_2{\mathsf{\Lambda }}_1+{\mathsf{\Lambda }}_2{\mathsf{\Phi }}_1\right)\right)\hfill \\ \\ & & +{\mathsf{\Phi }}_1\left({\mathsf{\Lambda }}_2{\varphi }_2-{\mathsf{\Lambda }}_3{\mathsf{\Phi }}_2\right)+{\mathsf{\Phi }}_4\left({\beta }_1\left({\beta }_2{\mathsf{\Lambda }}_1+{\mathsf{\Lambda }}_2{\mathsf{\Phi }}_1\right)+{\mathsf{\Lambda }}_3{\mathsf{\Phi }}_1{\mathsf{\Phi }}_2\right),\hfill \\ \\ \hfill K_4& =& \left({\beta }_2{\mathsf{\Lambda }}_1+{\mathsf{\Lambda }}_2{\mathsf{\Phi }}_1\right)\left({\beta }_1{\varphi }_3+{\mathsf{\Phi }}_3{\varphi }_2\right)+{\mathsf{\Phi }}_2\left({\mathsf{\Lambda }}_1{\mathsf{\Phi }}_3{\varphi }_1+{\mathsf{\Lambda }}_3{\mathsf{\Phi }}_1{\varphi }_3\right),\hfill \\ \\ \hfill K_5& =& {\mathsf{\Phi }}_3{\mathsf{\Phi }}_4{\tau }_1\left({\beta }_2{\mathsf{\Lambda }}_1+{\mathsf{\Lambda }}_2{\mathsf{\Phi }}_1\right)+{\mathsf{\Phi }}_3\left(\left({\beta }_2{\mathsf{\Lambda }}_1+{\mathsf{\Lambda }}_2{\mathsf{\Phi }}_1\right)\left({\tau }_3{\varphi }_2-{\tau }_1{\psi }_2\right)+{\mathsf{\Lambda }}_1{\varphi }_1\left({\mathsf{\Phi }}_2{\tau }_3+{\tau }_1{\psi }_1\right)\right)\hfill \\ \\ & & +{\mathsf{\Phi }}_4{\tau }_2\left({\beta }_1\left({\beta }_2{\mathsf{\Lambda }}_1+{\mathsf{\Lambda }}_2{\mathsf{\Phi }}_1\right)+{\mathsf{\Lambda }}_3{\mathsf{\Phi }}_1{\mathsf{\Phi }}_2\right)+{\tau }_3{\varphi }_3\left({\beta }_1\left({\beta }_2{\mathsf{\Lambda }}_1+{\mathsf{\Lambda }}_2{\mathsf{\Phi }}_1\right)+{\mathsf{\Lambda }}_3{\mathsf{\Phi }}_1{\mathsf{\Phi }}_2\right)\hfill \\ \\ & & +{\beta }_1{\mathsf{\Lambda }}_1{\tau }_2\left({\psi }_1{\varphi }_1-{\beta }_2{\psi }_2\right)+{\mathsf{\Lambda }}_3{\mathsf{\Phi }}_1{\psi }_1\left({\tau }_1{\varphi }_3-{\tau }_2{\varphi }_2\right)\hfill \\ \\ & & +{\tau }_2{\psi }_2\left({\mathsf{\Lambda }}_1\left({\beta }_2{\varphi }_2+{\mathsf{\Phi }}_2{\varphi }_1\right)-{\mathsf{\Phi }}_1\left({\beta }_1{\mathsf{\Lambda }}_2+{\mathsf{\Lambda }}_3{\mathsf{\Phi }}_2\right)\right)\hfill \\ \\ & & +{\psi }_2\left({\mathsf{\Lambda }}_2{\mathsf{\Phi }}_1{\tau }_2{\varphi }_2-{\tau }_1{\varphi }_3\left({\beta }_2{\mathsf{\Lambda }}_1+{\mathsf{\Lambda }}_2{\mathsf{\Phi }}_1\right)\right).\hfill \end{array}$$

The real values considered in the fitting model are used to obtain the equilibrium point, which is unique and positive and given by

$$\begin{array}{c}\hfill P_1=(218,750,20,760.7,160,653,71,623.1,1.21144\times {10}^6).\end{array}$$

(18)

We then present at $P_1$ that the breast cancer model is locally asymptotically stable.

 Theorem 1.

The breast cancer model (6) is locally asymptotically stable.

Proof. 

We obtain at $P_1$ the Jacobian matrix:

$$\begin{array}{c}\hfill J=\left(\begin{array}{ccccc}-{\mathsf{\Phi }}_1& 0& 0& 0& 0\\ {\beta }_2& -{\mathsf{\Phi }}_2& 0& {\psi }_1& 0\\ 0& {\beta }_1& -{\mathsf{\Phi }}_3& {\psi }_2& 0\\ {\varphi }_1& {\varphi }_2& {\varphi }_3& {\psi }_2-{\mathsf{\Phi }}_4& 0\\ 0& {\tau }_1& {\tau }_2& {\tau }_3& -{\mu }_3\end{array}\right).\end{array}$$

(19)

The characteristic equation for J has two negative eigenvalues, $-{\mathsf{\Phi }}_1$ and $-{\mu }_3$, while the remaining three eigenvalues with negative real part are calculated from the equation given by

$$\begin{array}{c}\hfill {\mathsf{\Theta }}^3+{\Xi }_1{\mathsf{\Theta }}^2+{\Xi }_2\mathsf{\Theta }+{\Xi }_3=0,\end{array}$$

(20)

where

$$\begin{array}{c}{\Xi }_1={\mathsf{\Phi }}_2+{\mathsf{\Phi }}_3+{\tau }_3+{\psi }_1,\hfill \end{array}$$

$$\begin{array}{ccc}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\Xi }_2\hfill & =& {\mathsf{\Phi }}_4\left[{\mathsf{\Phi }}_2(1-{\mathcal{R}}_2-{\mathcal{R}}_3)+{\mathsf{\Phi }}_3(1-{\mathcal{R}}_2-{\mathcal{R}}_4)\right]+{\mathsf{\Phi }}_2{\mathsf{\Phi }}_3,\hfill \end{array}$$

$$\begin{array}{ccc}{\Xi }_3\hfill & =& {\mathsf{\Phi }}_2{\mathsf{\Phi }}_3{\mathsf{\Phi }}_4\left(1-{\mathcal{R}}_0\right)\hfill \end{array}$$

(21)

where

$$\begin{array}{c}\hfill {\mathcal{R}}_0=\underset{{\mathcal{R}}_1}{\underbrace{\frac{{\beta }_1{\psi }_1{\varphi }_3}{{\mathsf{\Phi }}_2{\mathsf{\Phi }}_3{\mathsf{\Phi }}_4}}}+\underset{{\mathcal{R}}_2}{\underbrace{\frac{{\psi }_2}{{\mathsf{\Phi }}_4}}}+\underset{{\mathcal{R}}_3}{\underbrace{\frac{{\psi }_1{\varphi }_2}{{\mathsf{\Phi }}_2{\mathsf{\Phi }}_4}}}+\underset{{\mathcal{R}}_4}{\underbrace{\frac{{\psi }_2{\varphi }_3}{{\mathsf{\Phi }}_3{\mathsf{\Phi }}_4}}}.\end{array}$$

(22)

We can show that ${\Xi }_1{\Xi }_2-{\Xi }_3>0$,

$$\begin{array}{ccc}\hfill {\Xi }_1{\Xi }_2-{\Xi }_3& =& {\mathsf{\Phi }}_2^2\left[{\mathsf{\Phi }}_4(1-{\mathcal{R}}_2-{\mathcal{R}}_3)+{\mathsf{\Phi }}_3\right]+{\mathsf{\Phi }}_3{\mathsf{\Phi }}_4(1-{\mathcal{R}}_2-{\mathcal{R}}_4)\left({\mathsf{\Phi }}_3+{\tau }_3+{\psi }_1\right)\hfill \end{array}$$

$$\begin{array}{ccc}& & \hspace{1em}\hspace{1em}\hspace{1em} +{\mathsf{\Phi }}_2{\mathsf{\Phi }}_3\left({\mathsf{\Phi }}_4(1+{\mathcal{R}}_1-{\mathcal{R}}_2)+{\tau }_3+{\psi }_1\right)+{\mathsf{\Phi }}_2{\mathsf{\Phi }}_4(1-{\mathcal{R}}_2-{\mathcal{R}}_3)\left({\tau }_3+{\psi }_1\right)\hfill \end{array}$$

$$\begin{array}{l}+{\mathsf{\Phi }}_2{\mathsf{\Phi }}_3^2>0.\hfill \end{array}$$

(23)

This ensures that the breast cancer model is locally asymptotically stable at $P_1$. □

4. Numerical Scheme

In this section, we present the numerical algorithm for the numerical solution of the breast cancer model (6). In this part, we explain the history of the parameterized method. When the derivative is classical, we show precisely how to calculate it. It is worth noting that this method of employing the classical derivative was originally presented in [23]. The following generic nonlinear equation is used to derive this approach using the procedure given in [24]:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{y}^{\prime }\left(t\right)=f(t,y\left(t\right))\\ y\left(0\right)=y_0\end{array}.\right.\end{array}$$

(24)

We obtain the following after applying the integral for both sides:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}y\left(t\right)=y\left(0\right)+{\int }_0^{t}f(\theta ,y\left(\theta \right))d\theta \\ y\left(0\right)=y_0\end{array}.\right.\end{array}$$

(25)

Let us consider $t=t_{n+1}$ and $t=t_n$, to obtain after subtracting,

$$\begin{array}{c}\hfill y\left(t_{n+1}\right)=y\left(t_n\right)+{\int }_{t_n}^{t_{n+1}}f(\theta ,y\left(\theta \right))d\theta .\end{array}$$

(26)

Within the interval $\left[t_n,t_{n+1}\right]$, we approximate the function $f(\theta ,y(\theta \left)\right)$ and have

$$\begin{array}{c}\hfill f(\theta ,y\left(\theta \right))\simeq \left(1-\frac{1}{2\varpi }\right)f\left(t_n,y_n\right)+\frac{1}{2\varpi }f\left(t_{n+1},y_{n+1}\right).\end{array}$$

(27)

After the replacement, we have

$$\begin{array}{c}\hfill y_{n+1}=y_n+(1-\frac{1}{2\varpi })hf(t_n,y_n)+\frac{h}{2\varpi }f(t_{n+1},y_{n+1}).\end{array}$$

(28)

As a result of the observation that the obtained results led to an implicit scheme that was occasionally burdensome when undertaking simulation, the $y_{n+1}$ and $t_{n+1}$ on the opposite side of the equation were substituted, and their formulas were specified in [23].

As a result,

$$\begin{array}{c}\hfill y_{n+1}=y_n+(1-\frac{1}{2\varpi })hf(t_n,y_n)+\frac{h}{2\varpi }f({\overline{t}}_{n+1},{\overline{y}}_{n+1}),\end{array}$$

(29)

where ${\overline{t}}_{n+1}=t_n+\varpi h$. We have the following for the function ${\overline{y}}_{n+1}$,

$$\begin{array}{c}\hfill {\overline{y}}_{n+1}=y_n+{\int }_{t_n}^{t_n+\varpi h}f(\theta ,y\left(\theta \right))d\theta .\end{array}$$

(30)

Using the Euler approach, we obtain

$$\begin{array}{c}\hfill {\overline{y}}_{n+1}=y_n+\varpi hf(t_n,y_n).\end{array}$$

(31)

Thus, finally, the scheme for the general case is given by

$$\begin{array}{c}\hfill y_{n+1}=y_n+(1-\frac{1}{2\varpi })hf(t_n,y_n)+\frac{h}{2\varpi }f(t_n+h,y_n+\varpi f(t_n,y_n)).\end{array}$$

(32)

We then use this scheme for our proposed model (6) and have the following:

$$\begin{array}{c}{X}_{n+1}=X_n+(1-\frac{1}{2\varpi })h{f}_1(t_n,X_n,B_n,C_n,R_n,E_n)+\frac{h}{2\varpi }{f}_1(t_n+h,\hfill \\ X_n+\varpi f_1(t_n,X_n,B_n,C_n,R_n,E_n),B_n+\varpi f_1(t_n,X_n,B_n,C_n,R_n,E_n),\hfill \\ C_n+\varpi f_1(t_n,X_n,B_n,C_n,R_n,E_n),R_n+\varpi f_1(t_n,X_n,B_n,C_n,R_n,E_n),\hfill \\ E_n+\varpi f_1(t_n,X_n,B_n,C_n,R_n,E_n)),\hfill \end{array}$$

$$\begin{array}{c}{B}_{n+1}=B_n+(1-\frac{1}{2\varpi })h{f}_2(t_n,X_n,B_n,C_n,R_n,E_n)+\frac{h}{2\varpi }{f}_2(t_n+h,\hfill \\ X_n+\varpi f_1(t_n,X_n,B_n,C_n,R_n,E_n),B_n+\varpi f_2(t_n,X_n,B_n,C_n,R_n,E_n),\hfill \\ C_n+\varpi f_2(t_n,X_n,B_n,C_n,R_n,E_n),R_n+\varpi f_2(t_n,X_n,B_n,C_n,R_n,E_n),\hfill \\ E_n+\varpi f_2(t_n,X_n,B_n,C_n,R_n,E_n)),\hfill \end{array}$$

$$\begin{array}{c}{C}_{n+1}=C_n+(1-\frac{1}{2\varpi })h{f}_3(t_n,X_n,B_n,C_n,R_n,E_n)+\frac{h}{2\varpi }{f}_3(t_n+h,\hfill \\ X_n+\varpi f_3(t_n,X_n,B_n,C_n,R_n,E_n),B_n+\varpi f_3(t_n,X_n,B_n,C_n,R_n,E_n),\hfill \\ C_n+\varpi f_3(t_n,X_n,B_n,C_n,R_n,E_n),R_n+\varpi f_3(t_n,X_n,B_n,C_n,R_n,E_n),\hfill \\ E_n+\varpi f_3(t_n,X_n,B_n,C_n,R_n,E_n)),\hfill \end{array}$$

$$\begin{array}{c}{R}_{n+1}=R_n+(1-\frac{1}{2\varpi })h{f}_4(t_n,X_n,B_n,C_n,R_n,E_n)+\frac{h}{2\varpi }{f}_4(t_n+h,\hfill \\ X_n+\varpi f_4(t_n,X_n,B_n,C_n,R_n,E_n),B_n+\varpi f_4(t_n,X_n,B_n,C_n,R_n,E_n),\hfill \\ C_n+\varpi f_4(t_n,X_n,B_n,C_n,R_n,E_n),R_n+\varpi f_4(t_n,X_n,B_n,C_n,R_n,E_n),\hfill \\ E_n+\varpi f_4(t_n,X_n,B_n,C_n,R_n,E_n)),\hfill \end{array}$$

$$\begin{array}{c}{E}_{n+1}=E_n+(1-\frac{1}{2\varpi })h{f}_5(t_n,X_n,B_n,C_n,R_n,E_n)+\frac{h}{2\varpi }{f}_5(t_n+h,\hfill \\ X_n+\varpi f_5(t_n,X_n,B_n,C_n,R_n,E_n),B_n+\varpi f_5(t_n,X_n,B_n,C_n,R_n,E_n),\hfill \\ C_n+\varpi f_5(t_n,X_n,B_n,C_n,R_n,E_n),R_n+\varpi f_5(t_n,X_n,B_n,C_n,R_n,E_n),\hfill \end{array}$$

$$\begin{array}{ccc}& & E_n+\varpi f_5(t_n,X_n,B_n,C_n,R_n,E_n)),\hfill \end{array}$$

(33)

where

$$\begin{array}{ccc}\hfill f_1(t_n,X_n,B_n,C_n,R_n,E_n)& =& {\mathsf{\Lambda }}_1-({\varphi }_1+{\beta }_2)X,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill f_2(t_n,X_n,B_n,C_n,R_n,E_n)& =& {\mathsf{\Lambda }}_2+{\beta }_{2}X+{\psi }_{1}R-({\varphi }_2+{\beta }_1+{\tau }_1+{\mu }_1)B,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill f_3(t_n,X_n,B_n,C_n,R_n,E_n)& =& {\mathsf{\Lambda }}_3+{\beta }_{1}B+{\psi }_{2}R-({\varphi }_3+{\tau }_2+{\mu }_2)C,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill f_4(t_n,X_n,B_n,C_n,R_n,E_n)& =& {\varphi }_{1}X+{\varphi }_{2}B+{\varphi }_{3}C-({\psi }_1+{\psi }_2+{\tau }_3)R,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill f_5(t_n,X_n,B_n,C_n,R_n,E_n)& =& {\tau }_{3}R+{\tau }_{2}C+{\tau }_{1}B-{\mu }_{3}E.\hfill \end{array}$$

(34)

5. Numerical Simulation

Here, we consider the data given in [25]. The cancer patients (females) reported in KSA from 2004–2016 are considered in the present study. The time unit considered in data fitting is per year. For the analysis of the data fitting, we used the least-square curve-fitting method. We consider the numerical values of the model variables $X\left(0\right)=$ 30,000, $B\left(0\right)=$ 12,300, $C\left(0\right)=783$, $R\left(0\right)=334$ and $E\left(0\right)=10$. While simulating the model with these values, we obtain the required numerical results shown in Figure 1 and Figure 2. In these simulations, the parameter values that are fitted to the cancer data are shown in

Table 2. Details of the parameters.

Symbol	Value	Ref.	Symbol	Value	Ref.
${\mathsf{\Lambda }}_1$	14,000	Estimated	${\mathsf{\Lambda }}_2$	80	Estimated
${\mathsf{\Lambda }}_3$	90	Fitted	${\beta }_1$	0.01	Fitted
${\beta }_2$	0.034	Fitted	${\psi }_1$	0.03	Fitted
${\psi }_2$	0.3	Fitted	${\tau }_1$	0.03	Fitted
${\tau }_2$	0.1	Fitted	${\tau }_3$	0.2	Fitted
${\mu }_1$	0.0256	[26]	${\mu }_2$	0.0256	[26]
${\mu }_3$	0.0256	[26]	${\varphi }_1$	0.03	Fitted
${\varphi }_2$	0.4	Fitted	${\varphi }_3$	0.01	Fitted

. Figure 1 shows the data fitting versus the model simulations for the patients with cancer at stage 4. The data fit well with the model and can be useful to predict future cases in the country. A prediction of the cases versus model for a long time is shown in Figure 2. The result of the long-time behavior shows a better fitting result to the data and also predicts the disease to remain in the population for a long time until the government takes steps regarding the minimization of cancer cases.

We consider the simulation of the model compartments and obtain their graphical results. Figure 3 shows the dynamics of the model compartments for a long time ($t=120$). It can be seen that at the present stage if there is no prevention or possible treatment of cancer patients then the number of breast cancer patients will be increased, which will be an alarming situation for the country.

If there is an increase the in the value of ${\tau }_1$, then the number of patients over a long time will be decreased in the population of classes B, C, and R, while it will increase a little in class E (see Figure 4). For example, the parameter ${\tau }_1$ shows intensive chemotherapy and its good impact on the populations of stages 3 and 4 and those with disease-free conditions, but it does not have a good impact on the patients with cardiotoxicity (see Figure 4). Cancer therapy has advanced dramatically in recent years, and it has been shown to raise the rate of cure and prevent recurrences in breast cancer. The danger of cardiotoxicity, however, limits the application of these medications [27]. Cardiotoxicity is one of the most significant side effects of chemotherapy, significantly raising mortality as well as morbidity rates [28,29], which is one of the reasons that a small increase in Figure 4d is observed.

Figure 5 describes the impact of intensive chemotherapy on the individuals at stage 4, which can cause an increase in morbidity and mortality in the population of those who have cardiotoxicity. Increasing the chemotherapy treatment of individuals in stages 3 and 4, individuals with disease-free condition, and those who have cardiotoxicity, we can observe that the stage 3 patients have little improvement (see Figure 5a), and stage 4 and disease-free condition have good improvement in disease reduction (see Figure 5b,c), while for the case of cardiotoxicity, a good improvement is not observed. This is due to chemotherapy’s bad impact on people that increases the chances of affecting their heart-related complications [27,28,29,30].

Figure 6 describes the simulation of the model compartment when there is intensive chemotherapy in the disease-free case. The results indicate a small increase in the cardiotoxicity population, while some decreases can be seen in the populations of B, C, and R. For example, we can see the result in Figure 6 showing that there is little improvement in the decrease of the cases in stage 3 patients; similarly, in Figure 6b of the stage 4 patients, a good decrease in the cases is observed. The population with disease-free condition has also a good decrease in the number of patients (see Figure 6c), while in Figure 6d, a minor increase in the number of cases of cardiotoxicity is observed.

Figure 7 suggests the simulation of the cardiotoxicity individuals with the variations in the chemotherapy parameters of stages 3 and 4. It can be observed from Figure 7 that the number of cardiotoxicity individuals is decreasing with the decrease in the chemotherapy.

Usually, chemotherapy increases the risk of cardiotoxicity; however, according to certain studies, patients without a known cardiovascular history may experience symptoms of heart failure in direct relation to the total amount administered. This finding has led to the adoption of lower chemotherapy doses, which has decreased their effectiveness [31]. Cardiologists and oncologists should assess the patient for cardiovascular risk before the start of cancer therapy and decide on the best course of action based on the findings. This will help prevent cardiotoxicity [32].

6. Conclusions

In the present work, we considered breast cancer and its model formulation. We studied the model and obtained its dynamic results. The existence and uniqueness of the model were proven. We initially discussed the formulation of the model in detail and later carried out the mathematical results. We explored the equilibrium points and discussed their stability analysis. The data of Saudi Arabian female breast cancer patients was considered. In addition, with the fitted parameters, we obtained a unique positive endemic equilibrium. The equilibria of the model were examined and found to be locally asymptotically stable. Further, we utilized breast cancer data in Saudi Arabia and presented the data fitting with the model. The data fitting to the model is reliable and useful for future prediction of the disease in the country. We considered a new numerical approach for the ordinary nonlinear differential equations using the concept of the parameterized approach and presented a numerical scheme for its solutions. The novel numerical scheme was utilized to obtain various graphical results.

Some recommendations were given for the elimination of the disease in the country, i.e., chemotherapy treatments should be started as soon as possible whenever a person is identified as a cancer patient. This model shows that the risk of cardiotoxicity cases can be decreased if the cancer patients, before the chemotherapy, are treated according to their disease history. This study focused on the real data of female breast cancer patients in KSA and obtained the results. The parameters and the initial values of the variables were limited to the dataset used. In the future, this model shall be extended to add more information regarding cardiotoxicity and how cardiotoxicity can be minimized during chemotherapy at all stages of the disease.