Next Article in Journal
Determinants of Longitudinal Change of Lung Function in Different Gender in a Large Taiwanese Population Follow-Up Study Categories: Original Investigation
Previous Article in Journal
Clinical Significance of Peritumoral Adipose Tissue PET/CT Imaging Features for Predicting Axillary Lymph Node Metastasis in Patients with Breast Cancer
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Mathematical Model of Breast Tumor Progression Based on Immune Infiltration

1
Department of Mathematics and Statistics, University of Massachusetts Amherst, Amherst, MA 01003, USA
2
Department of Ophthalmology, Ohio State University Comprehensive Cancer Center, Columbus, OH 43210, USA
3
Department of Civil and Environmental Engineering, University of Massachusetts, Dartmouth, MA 02747, USA
4
Department of Biomedical Engineering and OHSU Center for Spatial Systems Biomedicine (OCSSB), Oregon Health and Science University, Portland, OR 97239, USA
5
Department of Mathematics, The Pennsylvania State University, University Park, PA 16802, USA
6
Mathematical NeuroOncology Lab, Precision Neurotherapeutics Innovation Program, Mayo Clinic Arizona, Phoenix, AZ 85054, USA
7
Department of Pharmacology and Chemical Biology, University of Pittsburgh, Pittsburgh, PA 15219, USA
8
Department of Bioengineering, UPMC Hillman Cancer Center, University of Pittsburgh, Pittsburgh, PA 15219, USA
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
J. Pers. Med. 2021, 11(10), 1031; https://doi.org/10.3390/jpm11101031
Submission received: 8 September 2021 / Accepted: 12 October 2021 / Published: 15 October 2021

Abstract

:
Breast cancer is the most prominent type of cancer among women. Understanding the microenvironment of breast cancer and the interactions between cells and cytokines will lead to better treatment approaches for patients. In this study, we developed a data-driven mathematical model to investigate the dynamics of key cells and cytokines involved in breast cancer development. We used gene expression profiles of tumors to estimate the relative abundance of each immune cell and group patients based on their immune patterns. Dynamical results show the complex interplay between cells and molecules, and sensitivity analysis emphasizes the direct effects of macrophages and adipocytes on cancer cell growth. In addition, we observed the dual effect of IFN- γ on cancer proliferation, either through direct inhibition of cancer cells or by increasing the cytotoxicity of CD8+ T-cells.

1. Introduction

Breast cancer is one of the most common cancers in women, and it is estimated that about 43,600 women in the United States will die because of breast cancer in 2021 [1]. There are different subtypes of breast cancer based on molecular-level analysis of gene expression patterns, such as luminal A (LumA), luminal B (LumB), and triple-negative/basal-like [2]. LumA is the most frequently seen subtype of breast cancer with the lowest mortality rate among other subtypes [3], and triple-negative breast cancer (TNBC) is the most aggressive subtype with a poor overall clinical outcome [4]. Several treatment options are available for breast cancer, including surgery (lumpectomy and mastectomy), radiotherapy, chemotherapy, hormone therapy, and other novel therapies that are determined by clinicopathological aspects of each patient [5]. Although targeted therapies have proven effective for some types of breast cancer [6], chemotherapy remains the standard form of treatment for subtypes of TNBC [7], and the development of resistance to chemotherapy is the main challenge for TNBCs [8,9].
It has been shown in many studies that the tumor’s microenvironment has a significant role in breast cancer development, progression, and also response to therapies [10,11,12,13]. Immune cells interact with tumor cells directly or indirectly via chemokine and cytokine signaling in the tumor microenvironment, and they are the major players in the behavior of the tumor and the efficiency of the treatments [14]. When cancer cells die, they release a protein called high mobility group box-1 (HMGB-1) that causes dendritic cells to become activated [15]. Antigens presented by dendritic cells cause T-cells to become activated and kill cancer directly [16,17,18]. IFN- γ is produced by helper T-cells, cytotoxic T-cells [19,20,21], and dendritic cells [19] to inhibit tumor growth [21]. On the other hand, certain immune cells can help cancer progress by either promoting it or having a dual impact. For example, regulatory T-cells control the development and activity of helper and cytotoxic T-cells, hence reducing the immune response and indirectly encouraging malignancy [22,23].
The relationship between clinical outcome and immune cells in breast tumor has been found in many studies. Breast cancer is distinguished by a significant population of tumor associated macrophages [24], and it has been shown that higher macrophage density is related to a poor outcome [25]. Moreover, the presence of CD8+ T-cells has been linked to considerable reductions in the relative risks of death from different subtypes of breast cancer [26]; and in ER-negative cancers, CD4+ and CD8+ T lymphocytes are more closely related to better outcomes than in ER-positive tumors [27]. In addition, it has been found that advanced stage breast tumors have more T-reg cells and a lower ratio of T-helper/T-reg cells [28]. All this evidence suggests that the relative numbers of distinct immune cells, and their interaction network, play key roles in the initiation and progression of breast cancer. Thus, to effectively simulate cancer progression, it is important to divide patients into similar cohorts based on their tumor-infiltrating immune cells and investigate the tumor progression of each group independently. However, adding all immune cells to the model increases the complexity and uncertainty of it. We therefore only considered the above-mentioned key immune cells (macrophages, CD4+ T-cells, T-reg cells, cytotoxic cells, and dendritic cells, which have a huge role in activating T-cells) in the progression of breast cancer.
Many mathematical models have been developed to study the relationships among tumors’ initiation, dynamics, and therapeutic responses to discover the best therapeutic combinations and overcome drug resistance in diverse cancer types [29,30,31,32,33,34,35,36,37], including breast cancer [38,39,40,41,42]. Some of the mathematical models for breast cancer investigate the relationships among immune cells, tumor cells, and some treatments [36,43,44]. For example, the effects of trastuzumab on HER2 overexpressing breast cancer in a mouse model system have been studied by integrating mathematical and experimental models [43]. In addition, a mathematical model has been developed to investigate the interactions between the MCF-7 breast cancer cell line and some immune cells. These models [43,45] include only a few immune components, such as NK cells and CD8+ T-cells, similarly to the mathematical model developed for the treatment of the murine 4T1 TNBC cell line with a high-dose chemotherapy drug [36].
Some of the outstanding challenges in the mathematical modeling of cancers are the existence of many unknown parameters and the limited number of datasets. For this reason, it is a routine practice to assume some values for parameters (see, e.g., [46,47,48]), or use estimated parameters from other diseases or models in the mathematical modeling of cancer. For example, the parameter values obtained for sarcoidosis were used to estimate the parameters of a mathematical model for pancreatic cancer, and the values estimated for pancreatic cancer in the mathematical modeling of breast cancer [23,49,50]. In addition, in a mathematical model of breast cancer [51], the values of some of parameters were chosen in line with a mathematical model of tumor invasion not validated for breast cancer [52]. Hence, the available mathematical models for breast cancer only consider a small subset of immune cells, and they assume all breast tumors behave similarly as they use the same parameter values for all tumors. However, since the evolution of a breast tumor depends on its specific immune profile, it is better to first find such tumors’ immune variations and understand the mechanism of growth in the absence of treatment for each of these immune variations. For this reason, we present a mathematical model for the interaction network given in Figure 1. We used a system of ordinary differential equations (ODEs) to investigate the differences in tumor progression of patients with different immune profiles. We clustered breast tumors based on their estimated immune cell frequencies using their gene expression data. We then estimated the parameters of the mathematical model for each tumor group separately. We found parameter values by reviewing the available literature and estimating the rest using what data were available. Importantly, we then performed global sensitivity analysis on the non-dimensionalized system to find the most sensitive parameters and investigate their impacts. Lastly, we analyzed the dynamics of the tumors in each group and compared them with patients’ clinical information to explore the connections between the tumor microenvironment and the progression of breast cancer.

2. Materials and Methods

There are many players in the progression of breast tumors. However, to avoid too much complexity and to reduce the uncertainty of the model, we only considered very important players. The main cell types that we modeled were cancer cells, necrotic cells, T-cells, macrophages, dendritic cell, and adipocytes (Figure 1).

2.1. Interaction Network of Cells and Molecules—ODE Model

2.1.1. T-Cells

In this model, T-cells are divided into four subgroups of naive, helper, cytotoxic, and regulatory T-cells. Naive T-cells, denoted by T N , are not necessarily part of the tumor microenvironment, as they are usually activated within lymph nodes. We excluded them from the total number of cells in the microenvironment. Even though there are other methods, such as introducing non-linear terms in the ODEs to avoid unlimited exponential growth, making activation rates for other types of T-cells proportional to the number of naive T-cells was the most convenient way to create a controlled system given the complexity of our model. Thus, we summarize the equation for the dynamics of naive T-cells after deriving the equations of other types of T-cells. The variables T h , T C , and T r , respectively, denote the numbers of activated T-helper cells, cytotoxic T-cells, and T-reg cells.

CD4+ Helper T-Cells ( T h )

CD4+ T-cells are activated by dendritic cells [16], HMGB1 [53], and IL-12 [23,54]. CD4+ T-cells’ phenotype expression is also promoted by estrogen [19]. On the other hand, regulatory T-cells [22] and IL-10 [55] inhibit CD4+ T-cells. Therefore, we modeled the dynamics of T-cells using the following equation:
d [ T h ] d t = λ T h H [ H ] + λ T h D [ D ] + λ T h I L 12 [ I L 12 ] + λ T h E [ E ] [ T N ] δ T h T r [ T r ] + δ T h I L 10 [ I L 10 ] + δ T h [ T h ] .

Cytotoxic T-Cells ( T c )

Estrogen promotes the CD8+ T-cell phenotype’s expression [19,56]. Dendritic cells [18,57] and IL-12 activate naive CD8+ T-cells [23,54]. On the other hand, CD8+ T-cells’ function is suppressed by regulatory T-cells [23] and IL-10 [55]. Since natural killer (NK) cells also directly kill cancer cells, we assume this group includes both CD8+ T-cells and NK cells. Therefore, we modeled the dynamics of cytotoxic T-cells in the following way.
d [ T c ] d t = λ T c E [ E ] + λ T c D [ D ] + λ T c I L 12 [ I L 12 ] [ T N ] δ T c T r [ T r ] + δ T c I L 10 [ I L 10 ] + δ T c [ T c ] .

Regulatory T-Cells ( T r )

Dendritic cells stimulate the formation [17] and activation of regulatory T-cells [57]. Furthermore, estrogen enhances the actions of regulatory T-cells [19,58]. Hence, we used the following equation for the dynamics of T-reg cells.
d [ T r ] d t = λ T r D [ D ] + λ T r E [ E ] [ T N ] δ T r [ T r ] .

Naive T-Cells ( T N )

By combining the above-mentioned equations for the activation of naive T-cells and introducing independent naive T-cell production rate A T N , one can get the following equation for dynamics of naive T-cells:
d [ T N ] d t = A T N λ T h H [ H ] + λ T h D [ D ] + λ T h I L 12 [ I L 12 ] + λ T h E [ E ] [ T N ] λ T c E [ E ] + λ T c D [ D ] + λ T c I L 12 [ I L 12 ] [ T N ] λ T r D [ D ] + λ T r E [ E ] + δ T N [ T N ] .

2.1.2. Dendritic Cells (D)

Dendritic cells are activated by cancer cells [23] and HMGB1 [59,60,61]. Moreover, estrogen enhances the metabolism, proliferation, differentiation, and functionality of dendritic cells [19], and multiple factors induced by cancer cells may promote natural dendritic cell death [57]. Denoting A D N as the production rate of naive dendritic cells, we made the following system of equations for dynamics of naive dendritic cells ( D N ) and activated dendritic cells (D):
d [ D N ] d t = A D N ( λ D C [ C ] + λ D H [ H ] + λ D E [ E ] ) [ D N ] δ D N [ D N ] ,
d [ D ] d t = ( λ D C [ C ] + λ D H [ H ] + λ D E [ E ] ) [ D N ] ( δ D C [ C ] + δ D ) [ D ] .

2.1.3. Macrophages (M)

Since macrophages have many phenotypes and frequently change their phenotype, the breakdown of them into M1, M2, and other subsets would have tremendously complicated the model. Therefore, we modeled all activated macrophages as a single variable denoted by M. Tumor associated macrophages (TAMs) are activated by IL-10 [62,63]. IL-12 and IFN- γ activate M1 macrophages [24,62,64,65], and M2 macrophages are activated by IL-4 and IL-13, which are secreted by helper T-cells [62]. Moreover, estrogen exposure leads to alternative macrophage activation [19,58]. Denoting naive macrophages by M N , activated macrophages by M, and the production rate of naive macrophages by A M , we made the following system of equations for the dynamics of naive and activated macrophages.
d [ M N ] d t = A M λ M I L 10 [ I L 10 ] + λ M I γ [ I γ ] + λ M I L 12 [ I L 12 ] + λ M T h [ T h ] + λ M E [ E ] [ M N ] δ M N [ M N ] ,
d [ M ] d t = λ M I L 10 [ I L 10 ] + λ M I γ [ I γ ] + λ M I L 12 [ I L 12 ] + λ M T h [ T h ] + λ M E [ E ] [ M N ] δ M [ M ] .

2.1.4. Cancer Cells (C)

Since cancer cells proliferate at an abnormal rate, the proliferation of cancer cells is traditionally modeled by [ C ] ( 1 [ C ] / C 0 ) , where C 0 is the maximum capacity. In addition, IL-6 promotes the proliferation of cancer cells [66,67]. Additionally, adipocytes, releasing metabolic substrates, promote the proliferation of breast cancer cells [68]. On the other hand, activated CD8+ T-cells kill cancer cells [23,69], and IFN- γ initiates the elimination of cancer cells by inducing cell cycle arrest, apoptosis, and necroptosis [21]. The dynamics of cancer cells was modeled by the following equation.
d [ C ] d t = λ C + λ C I L 6 [ I L 6 ] + λ C A [ A ] ( 1 [ C ] C 0 ) [ C ] ( δ C T c [ T c ] + δ C I γ [ I γ ] + δ C ) [ C ] .

2.1.5. Cancer Associated Adipocytes (A)

The direct crosstalk of cancer cells with tumor-surrounding stromal components, such as tumor-surrounding adipocytes, promotes tumor progression [70]. We modeled the proliferation of adipocytes similarly to the cancer cells’ proliferation.
d [ A ] d t = λ A [ A ] 1 [ A ] A 0 δ A [ A ] .

2.1.6. Necrotic Cells (N)

Necrosis occurs when cells are under metabolic stress as their resources are being consumed by cancer cells [60]. Cells that go through the process of necrosis are denoted as necrotic cells. Since resources are limited in the cancer microenvironment, some cells will undergo necrotic cell death instead of other types of cell death [60,71]. Activated CD8+ T-cells kill cancer cells [23], and CD8+ cytotoxic T-cells produce IFN- γ , which then eliminates cancer cells [21]. Since a fraction of cancer cells can go through first becoming necrotic cells, the production rate of necrotic cells was modeled by the fraction ( α N C ) of dying cancer cells.
d [ N ] d t = α N C δ C I γ [ I γ ] + δ C T c [ T c ] + δ C [ C ] δ N [ N ] .

2.1.7. Molecules

The dynamics of above mentioned molecules were modeled in the following way.

HMGB1 (H)

Damage-associated molecular pattern (DAMP) molecules are danger signals that promote inflammation and immune responses once released from dead or stressed cells [72]. The DAMP molecule, high-mobility group box 1 (HMGB1), exerts immune promoting activity by inducing angiogenesis, proliferation, and invasiveness of cancer cells via recruiting immune inflammatory cells [59]. HMGB1 is secreted by mature dendritic cells [59,60,73], necrotic cells [23,73,74], macrophages, [73,75,76], natural killer (NK) cells (which behave like cytotoxic T-cells) [73,77,78], and cancer cells [23,59,60].
For this reason, the dynamics of HMGB1 was modeled by the following equation.
d [ H ] d t = λ H D [ D ] + λ H N [ N ] + λ H M [ M ] + λ H T c [ T c ] + λ H C [ C ] δ H [ H ] .

IL-12 ( I L 12 )

IL-12, which stimulates the growth and functions of T-cells, is involved in the differentiation of naive T-cells into helper T-cells. IL-12 is secreted by macrophages and dendritic cells [23,57,64]. Helper and cytotoxic T-cells also produce IL-12 [19]. We modeled the dynamics of IL-12 using the following equation.
d [ I L 12 ] d t = λ I L 12 M [ M ] + λ I L 12 D [ D ] + λ I L 12 T h [ T h ] + λ I L 12 T c [ T c ] δ I L 12 [ I L 12 ] .

IL-10 ( I L 10 )

IL-10, which inhibits protective immune response (helper and cytotoxic T-cells), is produced by macrophages [62,79], dendritic cells [57,80,81], T-reg cells [55,82], CD4+ helper T-Cells [19,83], CD8+ cytotoxic T-cells [19,82], and cancer cells [17]. Therefore, the dynamics of IL-10 was modeled in the following way.
d [ I L 10 ] d t = λ I L 10 M [ M ] + λ I L 10 D [ D ] + λ I L 10 T r [ T r ] + λ I L 10 T h [ T h ] + λ I L 10 T c [ T c ] + λ I L 10 C [ C ] δ I L 10 [ I L 10 ] .

Estrogen (E)

Adipocytes are the primary producers of estrogen [84,85]. In breast tumors of postmenopausal women, estrogen can reach levels orders of magnitude greater than the low levels circulated in the body [85]. In general, adipose tissues in the breasts, brain tissues, osteoblasts, and other tissues locally produce estrogen, which circulates throughout the body. This amount of reproduced estrogen in the body depends on the pre-exiting existing amount. For this reason, we modeled the production rate of estrogen throughout the body using λ E [ E ] ( 1 [ E ] E 0 ) . This gave us the following equation for the dynamics of estrogen.
d [ E ] d t = λ E A [ A ] + λ E [ E ] ( 1 [ E ] E 0 ) δ E [ E ] .

IFN- γ ( I γ )

CD8+ T-cells and CD4+ T-cells release IFN- γ [19,20,21]. Dendritic cells also secrete IFN- γ , but when exposed to estrogen, this production is increased [86]. Therefore, we modeled the dynamics of IFN- γ in the following way.
d [ I γ ] d t = λ I γ T c [ T c ] + λ I γ T h [ T h ] + λ I γ D [ E ] [ D ] δ I γ [ I γ ] .

IL-6 ( I L 6 )

The important cytokine that leads to proliferation of cancer cells is IL-6 and is secreted by cancer associated adipocytes [64,66,67], macrophages [23,62,64,65,87], and dendritic cells [19,23,57].
d [ I L 6 ] d t = λ I L 6 A [ A ] + λ I L 6 M [ M ] + λ I L 6 D [ D ] δ I L 6 [ I L 6 ] .

2.2. Data of the Model

2.2.1. Breast Cancer Patients’ Data

There are some popular tumor deconvolution methods used to estimate the percentages of different immune cell types from the gene expression profiles of tumors. Among these methods, CIBERSORTx [88] has shown great performance in several studies [71,89,90,91]. In this study, we applied CIBERSORTx B-mode to gene expression profiles of primary tumors of breast cancer patients obtained from the Cancer Genome Atlas (TCGA) project of breast cancer (BRCA) [92] and the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) cohort [93]. There are 1904 primary breast tumor samples in the METABRIC microarray data downloaded from cBioPortal [94] and 1218 primary breast tumors with RSEM normalized RNA-seq data in log 2 scale in the TCGA data downloaded from the University of California Santa Cruz (UCSC) Xena web portal [95]. Before performing CIBERSORTx B-mode on both datasets separately, we transformed TCGA data to the linear space and dropped samples of normal breast tissues. After estimation of cell proportions, we only considered cases with CIBERSORTx B-mode p-values < 0.05 and continued our study with 2993 patients. We used expression values of genes encoding the molecules in the dynamical model and combined some immune cells frequencies to estimate the values of model’s variables as described in Table 1. Note that the genes’ expression levels in TCGA data were scaled depending on METABRIC data to eliminate the effect of different ranges in the two datasets on the results.

2.2.2. Patient Data Analysis

We grouped patients based on their estimated immune cell fractions using K-means clustering and determined the value of K using the elbow method. As a result, there were five distinct immune patterns of breast tumors. Figure 2 shows the average cell fractions in tumors of each cluster for the most variant immune cells among clusters, plus the standard deviation of each cluster as a bar.
Since the deconvolution method only provides us the percentage of each immune cell type in primary tumors, we used tumor weight for METABRIC data and tumor size for TCGA data, as described below, to estimate the numbers of immune cells, cancer cells, and necrotic cells in each tumor. For numerical stability, if the percentage of an immune cell that was used in the model was zero, it was substituted with 10% of the smallest positive cell fraction in the deconvolution data. We also excluded patients if their necrosis percentage or their tumor size in TCGA data, or the weight of their tumor in METABRIC data, were not available.
First, for simplicity, we assumed that the average number of cancer cells was twice the average number of total immune cells [37]. Therefore, using the necrosis percentage given in the TCGA data, the average ratio of immune cells to cancer cells to necrotic cells is approximately 0.3:0.6:0.1 in breast tumors. Additionally, the epithelial cell density has been reported as 4.5 × 10 4 cells/cm 3 in breast cancer [96]. We therefore choose the scaling factor α = 4.5 × 10 4 so that the average density of cancer cells across all patients was close to that value.
For each patient P t c g a in TCGA data, the total cell number ( TCN t c g a ) was assumed to be proportional to the weight of the tumor:
TCN t c g a = α tumor weight ( P t c g a ) 1 K all P t c g a tumor weight ( P t c g a )
where K is the number of patients in TCGA data. For TCGA data, numbers of necrotic cells ( N tcga ), cancer cells ( C tcga ), and total immune cells ( TIC tcga ) were calculated using the given necrotic percentage ( N p ) for each patient:
N tcga = TCN t c g a N p , C tcga = 2 3 TCN t c g a ( 1 N p ) and TIC tcga = 0.5 C tcga .
For METABRIC data, tumor size was used to calculate total cell number ( TCN m e t a ) similarly to the TCGA data:
TCN m e t a = α tumor size ( P m e t a ) 1 K all P m e t a tumor size ( P m e t a ) .
where K is the number of patients in the METABRIC data. Since the necrotic percentage was not available for METABRIC data, we used immune cells proportion from deconvolution data to calculate the total number of immune cells ( TIC meta ) for each patient ( P m e t a ):
TIC meta = 0.3 α Immune Cells Ratio ( P m e t a ) 1 K all P m e t a Immune Cells Ratio ( P m e t a ) so that
C meta = 6 7 ( TCN m e t a TIC meta ) and N meta = C meta 6 .
Once we have the amount of all the cells and molecules we can extract useful information such as initial conditions (Table 2) and steady state values (Table 3) for each cluster.

2.2.3. Parameter Estimation

We used 17 equations and 75 parameters in our model, which needed to be determined. We found the values of δ T N ,   δ T c ,   δ T h ,   δ T r ,   δ D ,   δ M ,   δ A ,   δ H ,   δ E ,   δ I L 6 ,   δ I L 10 ,   δ I L 12 and δ I γ by finding their corresponding cell and molecule half-lives in the biological literature. For the rest of the parameters, we derived them so that the dynamics of all cells and molecules reached their steady state within our simulation time. In other words, for an ODE of the type d X ¯ d t = F ( X ¯ , θ , t ) , we solved:
F ( X , θ , T ) = 0
for the parameter vector θ = θ 1 , , θ N , where T is the maximum simulation time and X is the vector of steady-state values for the state variables given in Table 3. However, even with all of the known death rate parameters, we still had 62 more to determine and only 17 steady state equations of type (18). To circumvent this issue, we added some assumptions which were basically relationships between the parameters. These relationships were based on assuming some production rates or death rates within the same ODE would be more effective than the rest.The relationship formulations were created by the authors to make sure we got a positive set of parameters. For more technical details of the parameter estimation process, please refer to Appendix A.1.

2.2.4. Sensitivity Analysis

As we mentioned, because of the lack of biological information about many of parameters’ values, we made several assumptions about these parameters given in Appendix A.1. To remedy this important limitation of the model, we perform a global gradient-based sensitivity analysis by changing each of these assumptions 5000 times, and obtaining a new set of parameter values with each new assumption. Since the number of parameter values (number of assumptions × number of variations = 38 × 5000 = 190,000 ) still is a finite number, this limitation of the model must be considered when the results of the model are used. For better numerical stability, we performed the sensitivity analysis using the non-dimensionalized system (see Appendix A.2 for more details).
Generally, the level of sensitivity of a variable X, to its vector of parameters θ = θ 1 , , θ N for an ODE of the type d X ¯ d t = F ( X ¯ , θ , t ) , is calculated by:
s i = d X ¯ d θ i , for i = 1 , , N .
We calculated the sensitivity of cancer cells and the total number of cells for each parameter θ i at the steady state. In other words, for the steady state values X ¯ * the sensitivity vector s = s 1 , , s N was obtained by differentiating F ( X ¯ * , θ ) = 0 with respect to θ . We got this analytical formula:
s = d X ¯ * d θ = F ( X ¯ * , θ ) 1 F ( X ¯ * , θ ) θ ,
where F ( X ¯ * , θ ) 1 is the numerical inverse of Jacobian of F with respect to X ¯ .
To calculate the global sensitivity, we followed the following steps:
  • First, we define a local sensitivity measure for each parameter θ i in the neighborhood Ω k ( θ i ) as
    S i k = Ω k ( θ i ) s i ( θ ) d θ .
    These neighborhoods were acquired by varying assumptions (A1)–(A22) by scaling factors 0.01 to 100. The integration was carried out numerically using sparse grid points [97,98].
  • Second, we found weights for the aforementioned neighborhoods. Scaling each assumption provides us with a new set of parameters. The weights were then determined by calculating the distance of each resulting parameter set to a fixed base parameter set. We assigned higher weights to the parameters that were closer to the base values. We denote each weight by w i k for i = 1 , , N and k = 1 , , K corresponding to the parameter and its neighborhood, respectively.
  • Finally, we obtained the global sensitivity level S i to each parameter θ i by
    S i = k = 1 K w i k S i k .

3. Results

3.1. Data Analysis of the Clusters

We also use TCGA and METABRIC clinical data to analyze the clinical features of each cluster and their associations with tumor immune microenvironment and dynamics. We see that cluster 2 and cluster 3, respectively, include a significant high number of the youngest and oldest patients compared to the other clusters (Figure 3A). We also observe that cluster 5 has the highest percentage of tumors with estrogen receptor (ER+) and human epidermal growth factor receptor 2 negative (HER2-) tumors, and ultimately a high number of LumA, LumB and normal-like tumors. On the other hand, clusters 1 and 4 have a higher percentage of aggressive subtypes such as HER2-enriched and Basal than other clusters, which explains having the highest percentage of tumors without estrogen receptor (ER-) in these clusters (Figure 4A,C,D). In addition, survival probabilities of cluster 1 and cluster 4 are lower than other clusters while cluster 5, which has more non-aggressive tumor, has the best survival results (Figure 3C). Note that survival status of the patients is given as alive or dead in TCGA data while in METABRIC data dead patients are separated into two categories ‘died of other causes’ and ‘died of disease’ so we translate ‘died of other causes’ into alive status. For this reason, we see survival probability of cluster 3 is lower in Figure 3C while the percentage of alive patients in cluster 3 is higher (Figure 4B).

3.2. Dynamics of the Breast Cancer Microenvironment

We find the dynamics of each variable involved in tumor microenvironment by solving the ODEs (1)–(17) with parameters acquired from the steady state assumptions for each cluster (Appendix A.1). We derive the dynamics of all of the cells and molecules based on the non-dimensional parameters in Table A4, half-lives and estimated death rates from Table A3, and initial conditions acquired from the tumors with the smallest size in each cluster (Table 2) and their steady state values in Table 3. We also plot the changes in total cells. We define total cells as:
Total Cells = T h + T c + T r + D N + D + 0.2 M N + M + C + N + A
We exclude T N since they are mainly found in the circulation and lymphatic system [99]. The rest of the T-cells are known as tumor infiltrating T-cells and can be found abundantly in tumors’ microenvironment. Dendritic cells primarily get activated inside of the tumor and cancer cells, necrotic cells and adipocytes are the other main components of the breast tumor, [12]. Finally, since most of naive macrophages polarize into M1 (anti-tumor) and M2 (pro-tumor) phenotypes inside of the tumor we consider a 20 % factor for M N ,100].
The dynamics of cells given in Figure 5 shows that for all clusters except cluster 2, the naive T-cells increase from their initial values and reach the steady state quickly, while the naive T-cells in cluster 2 increase and then decrease in a slow manner to reach the steady state. Helper T-cells increase from their initial condition and then decrease to reach the steady state quickly. Cytotoxic T-cells also follow the same pattern. T-reg cells reach the steady state quickly and remain constant at the steady state. While the naive dendritic cells start with a higher population initially and decrease eventually to reach the steady state, mature dendritic cells start with a low initial population and increase eventually to reach their steady states. This is because mature dendritic cells are derived from naive dendritic cells. On the other hand, both naive macrophages and macrophages increases within a short amount of time before reaching their steady states.
Cancer cells dynamics show an exponential growth until they reach the steady state. Since cancer cells go through necrosis, the necrotic cells’ growth behaves similarly. Adipocytes activation increases over time until they reach the steady state for all clusters. The amount of estrogen also increases from their initial conditions in each cluster as it is produced by adipocytes.
The initial concentration of HMGB1 in clusters 4 and 5 is higher than other clusters and stay higher in their steady state. The general trend for HMGB1 is increasing for all clusters before it reaches the steady state.
Although the concentration of IL-12 increases at the beginning, it lowers from the maximum concentration to reach the steady state. IL-10 and IL-6 increase from their initial value to reach the steady state. IFN- γ increases first then decreases over time as the cancer cells increase.
Due to the shear abundance of cancer cells and adipocytes compared to other cells, the total number of cells is more significantly affected by these cells. We can see that cluster 5 has one of the lowest cancer cells and the lowest adipocytes population, so the total cells population is the lowest in cluster 5. Additionally, cluster 2 has the highest cancer cells and one of the highest adipocytes in its steady state. As a result, total cells dynamics is the highest for cluster 2. Since total cells population is proportional to tumor sizes, it can be inferred that cluster 2 has larger tumors than the other clusters.
While cluster 3 has the fastest cancer growth, cluster 2 has the highest cancer population at the steady state among the other clusters. Additionally, cluster 2 has the highest population of dendritic cells that inhibit cancer cells and being activated by cancer cells. As dendritic cells are also activated by estrogen and HMGB1, their concentrations are the lowest in cluster 2. Thus, the higher population of dendritic cells may be mostly due to cancer cells’ activation. Note, cluster 2 includes the primary tumors of the youngest patients (Figure 3B) and aggressive tumors (Figure 4A), and it is known that aggressive breast tumor in younger patients grow faster [101].
On the other hand, cluster 4 and 5 have the slowest cancer growth, while cluster 1 has the lowest cancer population at the steady state compared to other clusters. Interestingly, the amount of immune cells that are known to be correlated with good prognoses such as cytotoxic T-cells and helper T-cells are the lowest in the cluster 1 which might cause making a wrong prediction at first glance. Thus, it is more important to consider interactions between the immune cells and cancer cells than just looking at their quantities to understand cancer prognosis.
Cluster 5 has the best survival probability according to Figure 3C. The population of adipocytes is the lowest in cluster 5 compared to the other clusters. Dendritic cells population, and the concentration of HMGB1 and estrogen are very low. Furthermore, cluster 5 has one of the lowest cancer cells population at the steady state. Furthermore, Figure 4C,D illustrate that cluster 5 has the lowest HER2-positive and ER-positive patients that are indicators of a better prognoses.
Although dynamical results show a lower cancer growth for cluster 1 and 4, clinical data analysis shows that cluster 1 and 4 have the lowest survival probability among the other clusters as illustrated in Figure 3C. This can be due to the fact that cluster 1 and 4 have one of the highest population of adipocytes at the steady state and thus providing more resources for cancer cells. Furthermore, based on Figure 4A these clusters have the highest percentage of aggressive tumors.

3.3. Sensitivity Analysis

After finding the base parameters from the steady state assumptions in Appendix A.1 and the steady state values given in Table 3, we carry out the global sensitivity analysis mentioned in Section 2.2.4.
We investigate the global sensitivity of cancer cells and the total number of cells (Equation (23)) to each parameter. The results are given in Figure 6 and Figure 7. As expected, we see the cancer population is the most sensitive to cancer related parameters (Figure 6), and since total cells includes high number of cancer cells (Figure 5), they are also sensitive to these parameters (Figure 6 and Figure 7). In addition, total cells and cancer populations in all clusters are significantly sensitive to δ A and λ A . This is reasonable, since all clusters have high numbers of adipocytes that directly contribute to cancer proliferation Equation (9). All of the sensitivity plots imply that less decay and more production of adipocytes lead to larger number of cancer and total cells.
Finally, cancer population in all clusters are directly sensitive to macrophage promoting parameters and reversely sensitive to the macrophage decay rates. We have modeled both macrophage phenotypes (M1 and M2) together. However, in all the clusters the frequency of M2 (pro-tumor) phenotype is more dominant (Figure 2). Therefore, more macrophages should result in more cancer cells and consequently higher number of total cells. The presence of δ M N as a sensitive parameter in Figure 7 can be justified in the same way as a significant number of naive macrophages turn into activated macrophages. Moreover, in clusters 1, 2, and 4, we see that clearance of necrotic cells has an opposite effect on the total number of cells which is reasonable, given the presence of the necrotic cells in the tumor microenvironment. In clusters 3 and 5, δ N is replaced by λ M I γ , because cluster 3 has a very high number of macrophages and cluster 5 has a high level of I γ , Figure 5.

3.4. Dynamics with Varying Assumptions

We investigate the effect of changing parameters’ assumptions on cancer dynamics in clusters 1–5. To avoid cluttering, we only target assumptions that incorporate sensitive parameters from the previous section, and we scale them by factors of 0.2 and 5. These scaled assumptions will produce new sets of parameters for each cluster leading to new dynamics. Furthermore, we locally perturb all the sensitive parameters from Section 3.3 to acquire a 10% variation region for cancer dynamics.
Notice that some sensitive parameters are either solely obtained from ODEs’ steady states or taken from the literature and are not explicitly included in these assumptions. Therefore, we only scale the following assumptions by 1.0, 0.2 and 5.0:
Scale × δ C T c [ T c max ] = 6 δ C I γ [ I γ max ]
Scale × δ C = 6 δ C I γ [ I γ max ]
Scale × 2 λ C A [ A max ] = λ C I L 6 [ I L 6 max ]
Scale × λ M I γ [ I γ max ] = λ M I L 10 [ I L 10 max ]
Scale × λ M E [ E max ] = λ M I L 10 [ I L 10 max ]
Scale × λ M T h [ T h max ] = λ M I L 10 [ I L 10 max ]
The results are shown in Figure 8. We can see that scaling the assumptions (24) and (25) had the largest effect on the dynamics of cancer cells across the five clusters. On the other hand, assumptions (26)–(29) negligibly affect the time of reaching steady state in clusters 4 and 5 when they are up-scaled. The assumption (24) delayed the time of reaching steady states for all clusters while it was scaled down and sped it up when it was scaled up. Clusters 2, 4, and 5 were the most affected in the case shown in Figure 8A; and clusters 1, 2, 4, and 5 were the most affected in the case shown in Figure 8B. Below is a comparison between the values of three crucial parameters for Figure 8A.
According to Table 4, downscaling the assumption (24) increased the value of δ C T c and decreased the values of δ C I γ and δ C . This is because both δ C and δ C T c are related to each other via δ C I γ . Since clusters 2, 4, and 5 had the highest levels of cytotoxic cells, increasing cytotoxic cells’ cancer inhibition effect ( δ C T c ) remarkably decreased the number of cancer cells. Similarly, we got Table 5 when scaling the assumption (25). In this case, upscaling the assumption (25) decreased the value of δ C and increased the values of δ C I γ and δ C T c . For the same reason, cancer saturation in clusters 2, 4, and 5 underwent a delay. Moreover, cluster 1 which has a low level of cytotoxic cells but a reasonable level of IFN- γ can more effectively remove cancer cells with this combination of parameter values.
None of the cases above completely neutralized the tumor. Even in the extreme cases, such as clusters 4 and 5 in Figure 8A,B, the cancer population growth was just delayed and would reach its steady state in later times. However, the mentioned parameters can cause significant delays in cancer growth, especially in clusters 1 and 4, whose patients had low survival probabilities. Hence, targeting cells and molecules corresponding to these parameters can be crucial in cancer treatments.
To even further investigate the parameters involved in cancer ODE and have more flexibility in finding scenarios in which cancer completely vanishes, we carried out a bifurcation analysis (see Appendix A.6). We assumed that [ I L 6 ] , [ A ] , [ T c ] and [ I γ ] were at their steady state values for each cluster. Next, we chose one parameter as the bifurcation variable and let the rest of them take their estimated values from Table A4 for each cluster.
Our results show that three parameters from Table 4 and Table 5 gave non-zero equilibria for cancer for very small values, and they caused cancer to vanish for a large range. This is in-line with how sensitive the dynamics results were to these parameters. On the other hand, parameters such as λ C , λ C I L 6 , and λ C A contribute to larger steady state values as they get larger. Saliently, cancer in clusters 1–5 vanished for very small values of λ C , while small values of λ C A could significantly decrease the steady state population of cancer only in cluster 5. This result is very interesting, because cluster 5 had the highest survival rate (Figure 3), the highest number of LumA subtypes (Figure 4), and the lowest adipocyte population (Figure 5). For more details, see Appendix A.6.

3.5. Dynamics with Different Initial Conditions

We investigated the effects of different initial conditions on the dynamics of each cluster. We varied the initial conditions by extracting individual patient data from each cluster. Due to the large number of patients, we plotted the dynamics of tumors with the number of cancer cells below the 40th percentile as the initial conditions. We observe that the dynamics did not change regardless of different initial conditions. All the patients from one cluster converged to the same steady state, no matter whether they started from higher or lower initial quantities (Figure ). This is reasonable, since the parameters were derived based on the steady state assumptions for each cluster.
We also looked at the dynamics of small tumors in one cluster converging on a large tumor at the steady state in another cluster. We observe that even when the initial conditions were not from the same cluster as the parameters’ values, the dynamics of tumors quickly converged on the dynamics of tumors in the cluster of the steady state. This testifies that the impact of parameters on dynamics is more significant than the initial conditions (Appendix A.5)s.

4. Discussion

Understanding cancer mechanisms and components in in vitro and in vivo studies is time consuming and does not provide comprehensive results to explain cancer complexity, since each cancer component is studied one at a time [102]. With the advancements in tumor deconvolution methods and the availability of more data [103], the data-driven mathematical models have become popular for exploring the complexity of the system more effectively. In this study, we developed a mathematical model that explored the characteristic differences of breast tumors based on their tumor microenvironments. This model allowed us to understand tumor progression in more detail.
The mathematical results showed that tumor growth in each cluster demonstrates unique interactions with its tumor microenvironment. For example, as the number of cancer cells increases, the numbers of cytotoxic cells, helper T-cells, and IFN- γ increase at the beginning and then decrease to reach a steady state. Furthermore, while regulatory T-cells in clusters 1, 2, and 3 reached the steady state very early and remained constant throughout the time, in clusters 4 and 5, they increased very early and then decreased to reach the steady state. Importantly, it has been found that interactions among immune cells can impact the immunological response in various cancer types [104,105]. Our results also show that it is crucial to investigate the interactions among immune cells rather than simply looking at the quantity of a certain immune cell type to make appropriate prognostic predictions for breast cancer patients. The results also show that some of the variables, such as T c , T h , D, E, and I γ , stayed approximately constant over time. Hence, one might simplify the model by removing the ODEs of these variables and assuming their steady state values in the remaining equations.
When we compare dynamical results with clinical information of the clusters, we can see that cluster 5, which has the best survival probability and one of the most non-aggressive tumor types (Figure 3C), had some of the lowest tumor growth and had a smaller cancer population at the steady state (Figure 5). Furthermore, we want to clarify that while cluster 3 contained mostly LumA and LumB breast tumors that tend to have better prognosis and better survival times compared to other subtypes of breast tumor, the survival probability of this cluster is low. This is because cluster 3 significantly consisted of the oldest patients compared to the other clusters (Figure 3A). In addition, while cluster 1 included more aggressive tumors than cluster 2, tumor growth in cluster 2 was the highest at the steady state. This might be due to the fact that cluster 2 had the youngest patients (Figure 3A), and aggressive breast tumors in young patients tend to grow more rapidly than ones in older women [101].
Having a large number of unknown parameters and difficulties in deriving their values are the main challenges in mathematical modeling of cancer. In this study, we derived some parameter values from experimental studies published in the literature, and we calculated others using the steady state assumptions and maximum values of the variables, using a similar method given in recent studies [71,106]. We estimated parameters uniquely for each group of patients with different immune profiles, because many studies have demonstrated that patients with different immune compositions show different prognoses and respond differently to therapy [107,108,109]. Thus, having separate parameters for each cluster allowed us to see the effects of immune cells on tumor growth more accurately.
In addition, we performed a global gradient-based sensitivity analysis to investigate the effects of each parameter on the dynamical system. The results of sensitivity and bifurcation analyses provided interesting insights about the biological mechanism, especially the importance of IFN- γ in controlling cancer growth (Figure 8 and Figure A3). We saw that the cancer population is primarily sensitive to parameters that directly promote or inhibit cancer cells, and secondarily, to parameters which promote or inhibit macrophages. In other words, cancer dynamics positively react to an increase in macrophages. It is reported that high levels of tumor associated macrophages are correlated with worse patient prognosis [100]. Moreover, to acquire patient specific parameters, we had to devise some mathematical assumptions about the parameters of our system. These assumptions were mostly mathematical artifacts ensuring nonzero parameter values. We chose the assumptions which involved the sensitive parameters in our sensitivity analysis and investigated the effect of scaling these assumptions to see how the cancer dynamics would change. We found out that we could delay cancer recovery time through certain assumptions. More specifically, we saw that the parameter δ C I γ was involved in the most interesting cases. In these cases, we saw a significant delay in cancer saturation time. Furthermore, according to our bifurcation analysis (see Appendix A.6), these parameters have a large range, which led to a vanishing cancer population that could translate into promising treatment options. The therapeutic potential of IFN- γ in later stages of cancer has been observed by researchers [110,111]. It is known that IFN- γ can directly inhibit breast cancer by restoring ICI-induced apoptosis in breast cancer cells that have acquired resistance to this antiestrogen [112], and the combination of anti-growth factor receptor MAb and cytokines such as IFN- γ may be useful in the treatment of breast cancer [113]. On the other hand, IFN- γ can increase the cytotoxicity of T-cells, leading to more effective inhibition of cancer cells [114,115]. Our bifurcation results also showed that controlling the production of cancer via adipocytes can significantly decrease the cancer population at the steady state for LumA subtypes. It has been observed that an excessive amount of adipocytes is associated with advanced T stages of the triple-negative and LumA subtypes, and the modulation of adipocytes can provide novel therapy targets for breast cancer [116]. Besides these, we also investigated how a dynamical system would be affected if we used different initial conditions for each cluster, and we saw that the parameter values that were used in each cluster were more important than the initial conditions for determining the dynamics.
The results of this study should be used by considering the limitations of this study. While the use of time course data would be ideal to obtain parameters for each cluster, the availability of such data is limited. To combat this challenge, we used large amounts of tumor data in each cluster as the steady state values, to estimate parameters. Despite this limitation due to the lack of time course data, our model still presented essential information about the progression of breast tumors, and our study can be used to develop new models using only gene expression data of the patients. Moreover, one might utilize different approaches, such as applying partial differential equations, to investigate tumor growth in space and time, or utilize other parameter fitting algorithms [117,118,119,120] to obtain more accurate dynamical systems. Lastly, it has been shown that treatment options in breast cancer have significant effects on the immune system. This model can be used for inferring treatment strategies in breast cancer by adding the interactions between the effects of the treatments on the tumor microenvironment.

Author Contributions

Conceptualization, M.H.A.-R., A.A., C.M.C., Y.H.C., W.H., P.R.J., A.V.L., D.G.S., Z.T., I.K.Z. and L.S.; methodology, N.M.M., S.S., D.S. and L.S.; software, N.M.M. and S.S.; validation, N.M.M., S.S. and D.S.; formal analysis, N.M.M., S.S., D.S. and M.H.; investigation, N.M.M., S.S. and D.S.; resources, L.S.; data curation, S.S.; writing—original draft preparation, N.M.M., S.S., D.S., M.H. and L.S.; writing—review and editing, N.M.M. and L.S.; visualization, N.M.M., S.S. and D.S.; supervision, L.S.; project administration, M.H.A.-R., A.A., C.M.C., Y.H.C., W.H., P.R.J., A.V.L., D.G.S., I.K.Z. and L.S.; funding acquisition, N.M.M., M.H.A.-R., A.A., C.M.C., Y.H.C., W.H., P.R.J., A.V.L., D.G.S., I.K.Z. and L.S. All authors have read and agreed to the published version of the manuscript.

Funding

This project has been funded in whole or in part with Federal funds from the National Cancer Institute, National Institutes of Health, under subcontract number 21X131F of Leidos Biomed’s prime, contract number 75N91019D00024, task order number 75N91019F00134; under award number U54CA209988; and by the Department of Energy under Award Number DE-SC0021630. The content of this publication does not necessarily reflect the views or policies of the funding agencies, nor does any mention of trade names, commercial products or organizations imply endorsement by the U.S. Government.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The TCGA datasets of BRCA project are available at https://portal.gdc.cancer.gov and https://xenabrowser.net/datapages/, accessed on 1 June 2021. The METABRIC dataset is available at https://www.cbioportal.org/datasets, accessed on 1 June 2021. Codes are available at our GitHub page https://github.com/ShahriyariLab/Mathematical-Model-of-Breast-Tumors-Progression-Based-on-Their-Immune-Infiltration, accessed on 1 June 2021.

Acknowledgments

We would like to thank the NCI-DOE 2020 Virtual Ideas Lab, Toward Building a Cancer Patient “Digital Twin” https://events.cancer.gov/cbiit/dtwin2020, accessed on 1 June 2021.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
TCGAThe Cancer Genome Atlas
METABRICMolecular Taxonomy of Breast Cancer International Consortium
HMGB1High mobility group box-1
LumALuminal A
LumBLuminal B
TNBCTriple negative breast cancer
IFN- γ Interferon gamma
HER2Human epidermal growth factor 2
EREstrogen receptor
DAMPDamage-associated molecular pattern
UCSCUniversity of Santa Cruz

Appendix A. Derivation of Sample Parameters

Appendix A.1. Steady State and Additional Assumptions

Since we were modeling the evolutions of tumors without any interventions, we assumed small tumors grow to large ones at a steady state. As there were both small and large tumors in each cluster, the averages of the smallest tumors and the largest tumors in each cluster, respectively, gave us the initial conditions and the steady state values for the cluster; see Table 2 and Table 3. However, we carried out our parameter estimation based on just the steady state values and independent of the initial conditions. This enabled us to validate the importance of our parameter values based on the data of every single patient in the same cluster and across clusters; see Appendices Appendix A.4 and Appendix A.5.
If we assume specific values of steady state for each variable, we obtain the following sample parameter set:
T N ,   T c ,   T h ,   T r ,   D N ,   D ,   M ,   M N ,   A ,   C ,   N ,   H ,   E ,   I L 6 ,   I L 10 ,   I L 12   and   I γ .
Under these assumptions, the system of Equations (1)–(17) provided us with 17 constraints on parameters. In order to solve this system, we needed to make additional assumptions about a selection of parameters. This ensured that the number of independent equations equaled the number of unknown parameters and that we could uniquely determine the values. We assumed that the carrying capacities A 0 , C 0 , and M 0 were proportional to the maximum values of their corresponding variables in each cluster. We also took the necrosis coefficient to be α N C = 0.5 . We also utilized the available research [71,121,122,123,124] to adopt the natural death/decay/degradation rates. Considering a cell or cytokine X, we estimated the death rate δ X = ln 2 / t 1 / 2 X using half-life measurements t 1 / 2 X .
More specifically, we found that naive T-cells’ half-life ranges between 1 and 8 years [71]. Taking it to be two years, we derived a death rate of 9.49 × 10 4 . Next, we found that the half-lives of cytotoxic T-cells and helper T-cells are 41 h and 3 days, respectively, resulting in δ T h = δ T r = 2.31 and δ T c = 0.406 [71]. The half-life of mature DCs is about 2 to 3 days, giving us an average of 2.5 days. Using this value for t 1 / 2 X , we obtained δ D = 0.277 . Intestinal macrophages have a half life ranging from 4 to 6 weeks. Considering the half life of five weeks, we calculated δ M = 0.0198 . According to [125], we found that the half life of adipocytes is 10 years, resulting in δ A = 2.8 × 10 3 . For HMGB1, there are varying results in [121]. We considered a half-life of 17 min and used this to find δ H = 18 . For the degradation rate of estrogen, we considered a 3 to 5 h half-life [122]. Using the average of 4 h, we calculated δ E = 4.16 . For IL-6, we assumed the half-life to be 15.5 h. Thus, we obtained δ I L 6 = 1.07 [71]. We took the half-lives of IL-10 and IL-12 to be 3.6 h [71] and 7.8 min [123], respectively. This resulted in δ I L 10 = 4.62 and δ I L 12 = 128 . Lastly, for IFN- γ , the estimated half-life was about 30 min, resulting in the rate δ γ = 33.3 . More concisely, we used the following the degradation rates in our computation:
δ T N = 9.49 × 10 4 δ T c = 0.406 δ T h = 0.231 δ T r = 0.231 δ D = 0.277 δ M = 0.0198 δ A = 2.8 × 10 3 δ H = 18 δ E = 4.16 δ I L 6 = 1.07 δ I L 10 = 4.62 δ I L 12 = 128 δ I γ = 33.3
Though this decreased the number of unknowns, there still existed a substantial number of parameters without values. Therefore, we next imposed heuristic assumptions based on activation, production, and inhibition rates. Let us look at those in detail.
We began by considering that the mean rate of doubling time of a breast cancer tumor is 193 ± 141 [126]. This gives an interval of 52 days to 334 days. Using this fact, we assumed that a faster doubling rate results from both innate cancer proliferation and proliferation caused by cytokine IL-6 and adipocytes. In the case of the faster double rate, we assumed that death is only innate. This results in
ln 2 52 1.3330 × 10 2 =   λ C + λ C I L 6 [ I L 6 mean ] + λ C A [ A mean ] δ C .
This contrasts with our slower doubling rate assumptions. We assumed only innate cancer proliferation. In this case, death rates included the effects of anti-cancer agents cytotoxic T-cells and IFN- γ , giving us
ln 2 334 2.0753 × 10 3 = λ C δ C T c [ T c mean ] + δ C I γ [ I γ mean ] + δ C .
We used the following mean values:
Table A1. Mean values of the variables used. The mean values were calculated from the combination of TCGA and METABRIC datasets for all of the clusters.
Table A1. Mean values of the variables used. The mean values were calculated from the combination of TCGA and METABRIC datasets for all of the clusters.
IL 6 mean A mean T c mean I γ mean
3.906 7.112 × 10 4 4.129 × 10 3 3.675
In these assumptions, we took I L 6 mean , A mean , T c mean , and I γ mean to be the average values of the corresponding variable across all patients. In our further assumptions, we used the maximum observable quantities for all the variables across all patients denoted by:
T N max , T c max , T h max , T r max , D N max , D max , M N max , M max , A max , C max , N max , H max , E max , I L 6 max , I L 10 max , I L 12 max , and I γ max .
The following maximum variables were calculated from the datasets:
Table A2. Maximum values of the variables. The maximum values were calculated from the combination of TCGA and METABRIC datasets for all of the clusters.
Table A2. Maximum values of the variables. The maximum values were calculated from the combination of TCGA and METABRIC datasets for all of the clusters.
T N max T c max T h max T r max D N max D max
4.270 × 10 4 1.599 × 10 4 2.495 × 10 4 1.174 × 10 4 6.507 × 10 3 8.274 × 10 3
M N max M max N max N max A max H max
4.577 × 10 4 5.247 × 10 4 4.258 × 10 5 1.019 × 10 5 2.612 × 10 5 1.160 × 10 1
I L 12 max I L 10 max E max I γ max I L 6 max
2.516 × 10 1 5.962 1.868 × 10 1 8.768 1.073 × 10 1
Assuming helper T-cells are primarily activated by dendritic cells, we can write
λ T h D [ D max ] = 200 λ T h H [ H max ] = 200 λ T h I L 12 [ I L 12 max ] = 200 λ T h E [ E max ] .
We took the inhibition of helper T-cells to be mostly accomplished by cytokine IL-10 and regulatory T-cells. In particular, we claim that these inhibitors are 20 times more effective than natural degradation.
δ T h T r [ T r max ] = δ T h I L 10 [ I L 10 max ] = 20 δ T h .
Next, we made the following assumptions about cytotoxic T-cells. The activation of cytotoxic T-cells by estrogen is much stronger than activation by IL-12 and dendritic cells. We declare estrogen to be 200 times as effective as IL-12, and we declare dendrites to be 100 times as effective as estrogen.
λ T c E [ E max ] = 100 λ T c D [ D max ] = 200 λ T c I L 12 [ I L 12 max ] .
We also state that IL-10 and regulatory T-cells are 20 times more effective than natural degradation of cytotoxic T-cells.
δ T c T r [ T r max ] = δ T c I L 10 [ I L 10 max ] = 20 δ T c .
For regulatory T-cells, we assumed that their activation by dendritic cells is four times stronger than their activation by estrogen. This is written as
λ T r D [ D max ] = 4 λ T r E [ E max ] .
We also claim that the activation of dendritic cells by HMGB1 and estrogen is twice as effective as by cancer cells. This gives us the following equality:
2 λ D C [ C max ] = λ D H [ H max ] = λ D E [ E max ] .
Next, we declare the inhibition of dendritic cells by cancer cells to be equivalent to the natural degradation rate. In other words,
δ D C [ C max ] = δ D .
Further, we assume that helper T-cells are the primary activator for macrophages. We express this assumption as
10 λ M I L 10 [ I L 10 max ] = 10 λ M I γ [ I γ max ] = 10 λ M I L 12 [ I L 12 max ] = λ M T h [ T h max ] = 10 λ M E [ E max ] .
For HMGB1, we declare the following equality to be true.
λ H D [ D max ] = λ H N [ N max ] = 10 λ H M [ M max ] = λ H T c [ T c max ] = λ H C [ C max ] .
Our next assumption sets the production of IL-12 to be equal among macrophages, dendritic cells, helper T-cells, and cytotoxic T-cells:
λ I L 12 M [ M max ] = λ I L 12 D [ D max ] = λ I L 12 T h [ T h max ] = λ I L 12 T c [ T c max ] .
For the production of IL-10, we let macrophages, dendritic cells, regulatory T-cells, helper T-cells, cytotoxic T-cells, and cancer cells work at equal rates:
λ I L 10 M [ M max ] = λ I L 10 D [ D max ] = λ I L 10 T r [ T r max ] = λ I L 10 T h [ T h max ] = λ I L 10 T c [ T c max ] = λ I L 10 C [ C max ] .
Next, we assumed that IFN- γ is mostly produced by cytotoxic T-cells. This gave us the following equality:
λ I γ T c [ T c max ] = 5 λ I γ T h [ T h max ] = 5 λ I γ D [ D max ] .
We also declare IL-6 to be equally produced by adipocytes, macrophages, and dendritic cells. Thus, we obtained the following:
λ I L 6 A [ A max ] = λ I L 6 M [ M max ] = λ I L 6 D [ D max ] .
Furthermore, we scaled the carrying capacities for cancer cells, adipocytes, and estrogen by 2, 2, and 1.5 , respectively.
C 0 = 2 × [ C max ] , A 0 = 2 × [ A max ] , and E 0 = 1.5 × [ E max ] .
We took the necrosis factor α N C to be 0.5 for all clusters.
α N C = 0.5 .
The proliferation rate of cancer cells by I L 6 was assumed to be twice that by adipocytes:
2 λ C A [ A max ] = λ C I L 6 [ I L 6 max ] .
We declare the general secretion of estrogen to be twenty times the production of estrogen by adipocytes:
20 λ E A [ A max ] = λ E .
We took the depreciation of naive dendritic cells to be equivalent to that of dendritic cells. Likewise, we took the depreciation of naive macrphage cells to be equivalent to that of macrophages.
δ D N = δ D
δ M N = δ M .
Finally, we assumed that the inhibition of cancer cells by IFN- γ is less effective than inhibition by cytotoxic cells and other causes:
δ C = δ C T c [ T c max ] = 6 δ C I γ [ I γ max ]
Altogether, these assumptions adequately allowed us to derive a sample parameter set.

Appendix A.2. Non-Dimensionalization

For a more stable numerical simulation, parameter estimation, and sensitivity analysis, we non-dimensionalized our parameters and variables with respect to each variable’s steady state value. Therefore, a non-dimensional variable [ X ¯ ] , non-dimensionalized as [ X ] / [ X ] , will take the value 1 at its steady state. Since the dimensional time scale does not introduce any inconveniences, we chose not to non-dimensionalize with respect to the time. As a result we got the following system:
d [ T h ¯ ] d t = λ ¯ T h H [ H ¯ ] + λ ¯ T h D [ D ¯ ] + λ ¯ T h I L 12 [ I L 12 ¯ ] + λ ¯ T h E [ E ¯ ] [ T N ¯ ] δ ¯ T h T r [ T r ¯ ] + δ ¯ T h I L 10 [ I L 10 ¯ ] + δ T h [ T h ¯ ] ,
d [ T c ¯ ] d t = λ ¯ T c E [ E ¯ ] + λ ¯ T c D [ D ¯ ] + λ ¯ T c I L 12 [ I L 12 ¯ ] [ T N ¯ ] δ ¯ T c T r [ T r ¯ ] + δ ¯ T c I L 10 [ I L 10 ¯ ] + δ T c [ T c ¯ ] ,
d [ T r ¯ ] d t = λ ¯ T r D [ D ¯ ] + λ ¯ T r E [ E ¯ ] [ T N ¯ ] δ T r [ T r ¯ ] ,
d [ T N ¯ ] d t = A ¯ T N λ ¯ T h H [ H ¯ ] + λ ¯ T h D [ D ¯ ] + λ ¯ T h I L 12 [ I L 12 ¯ ] + λ ¯ T h E [ E ¯ ] [ T N ¯ ] λ ¯ T c E [ E ¯ ] + λ ¯ T c D [ D ¯ ] + λ ¯ T c I L 12 [ I L 12 ¯ ] [ T N ¯ ] λ ¯ T r D [ D ¯ ] + λ ¯ T r E [ E ¯ ] + δ T N [ T N ¯ ] ,
d [ D N ¯ ] d t = A ¯ D N λ ¯ D C [ C ¯ ] + λ ¯ D H [ H ¯ ] + λ ¯ D E [ E ¯ ] [ D N ¯ ] δ D N [ D N ¯ ] ,
d [ D ¯ ] d t = λ ¯ D C [ C ¯ ] + λ ¯ D H [ H ¯ ] + λ ¯ D E [ E ¯ ] [ D N ¯ ] δ ¯ D C [ C ¯ ] + δ D [ D ¯ ] ,
d [ M N ¯ ] d t = A ¯ M ( λ ¯ M I L 10 [ I L 10 ¯ ] + λ ¯ M I γ [ I γ ¯ ] + λ ¯ M I L 12 [ I L 12 ¯ ] + λ ¯ M T h [ T h ¯ ] + λ ¯ M E [ E ¯ ] + δ M N ) [ M N ¯ ] ,
d [ M ¯ ] d t = λ ¯ M I L 10 [ I L 10 ¯ ] + λ ¯ M I γ [ I γ ¯ ] + λ ¯ M I L 12 [ I L 12 ¯ ] + λ ¯ M T h [ T h ¯ ] + λ ¯ M E [ E ¯ ] [ M N ¯ ] δ M [ M ¯ ] ,
d [ C ¯ ] d t = λ C + λ ¯ C I L 6 [ I L 6 ¯ ] + λ ¯ C A [ A ¯ ] 1 [ C ¯ ] C 0 ¯ [ C ¯ ] δ ¯ C T c [ T c ¯ ] + δ ¯ C I γ [ I γ ¯ ] + δ C [ C ¯ ] ,
d [ A ¯ ] d t = λ A [ A ¯ ] 1 [ A ¯ ] A 0 ¯ δ A [ A ¯ ] ,
d [ N ¯ ] d t = α ¯ N C δ ¯ C I γ [ I γ ¯ ] + δ ¯ C T c [ T c ¯ ] + δ C [ C ¯ ] δ N [ N ¯ ] ,
d [ H ¯ ] d t = λ ¯ H D [ D ¯ ] + λ ¯ H N [ N ¯ ] + λ ¯ H M [ M ¯ ] + λ ¯ H T c [ T c ¯ ] + λ ¯ H C [ C ¯ ] δ H [ H ¯ ] ,
d [ I L 12 ¯ ] d t = λ ¯ I L 12 M [ M ¯ ] + λ ¯ I L 12 D [ D ¯ ] + λ ¯ I L 12 T h [ T h ¯ ] + λ ¯ I L 12 T c [ T c ¯ ] δ I L 12 [ I L 12 ¯ ] ,
d [ I L 10 ¯ ] d t = λ ¯ I L 10 M [ M ¯ ] + λ ¯ I L 10 D [ D ¯ ] + λ ¯ I L 10 T r [ T r ¯ ] + λ ¯ I L 10 T h [ T h ¯ ] + λ ¯ I L 10 T c [ T c ¯ ] + λ ¯ I L 10 C [ C ¯ ] δ I L 10 [ I L 10 ¯ ] ,
d [ E ¯ ] d t = λ ¯ E A [ A ¯ ] + λ E [ E ¯ ] 1 [ E ¯ ] E 0 ¯ δ E [ E ¯ ] ,
d [ I γ ¯ ] d t = λ ¯ I γ T c [ T c ¯ ] + λ ¯ I γ T h [ T h ¯ ] + λ ¯ I γ D [ E ¯ ] [ D ¯ ] δ I γ [ I γ ¯ ] ,
d [ I L 6 ¯ ] d t = λ ¯ I L 6 A [ A ¯ ] + λ ¯ I L 6 M [ M ¯ ] + λ ¯ I L 6 D [ D ¯ ] δ I L 6 [ I L 6 ¯ ] .
Since we did not nondimensionalize with respect to time, the production rates λ C , λ A and λ E and the decay rates δ T h , δ T c , δ T r , δ T N , δ D N , δ D , δ M N , δ M , δ C , δ A , δ N , δ H , δ I L 12 , δ I L 10 , δ E , δ I γ , δ T h , and δ I L 6 were left unchanged. To non-dimensionalize the rest of the parameters, we used the following rules:
A ¯ T N = A T N [ T N ] , A ¯ D N = A D N [ D N ] , A ¯ M = A M [ M ] , α ¯ N C = α N C [ C ] [ N ] , C ¯ 0 = C 0 [ C ] , A ¯ 0 = A 0 [ A ] , E ¯ 0 = E 0 [ E ] .
Furthermore, for the reaction parameters of the form λ ¯ X Y used in terms of the form λ ¯ X Y [ Y ] [ Z ] , we used the general non-dimensinalization formula:
λ ¯ X Y = λ X Y [ Y ] [ Z ] [ X ] .
and finally, for inhibition rates of the form δ ¯ X Y used in terms of the form δ ¯ X Y [ Y ] [ Z ] , we used the general non-dimensinalization formula:
δ ¯ X Y = δ X Y [ Y ] .
Consequently, we got non-dimensionalized assumptions as follows:
ln 2 52 1.3330 × 10 2 = λ C + λ C I L 6 [ I L 6 mean ] [ I L 6 ] + λ C A [ A mean ] [ A ] δ C ,
ln 2 334 2.0753 × 10 3 = λ C δ C T c [ T c mean ] [ T c ] + δ C I γ [ I γ mean ] [ I γ ] + δ C ,
λ T h D [ D max ] [ D ] = 200 λ T h H [ H max ] [ H ] = 200 λ T h I L 12 [ I L 12 max ] [ I L 12 ] = 200 λ T h E [ E max ] [ E ] ,
δ T h T r [ T r max ] [ T r ] = δ T h I L 10 [ I L 10 max ] [ I L 10 ] = 20 δ T h , λ T r D [ D max ] [ D ] = 4 λ T r E [ E max ] E ,
2 λ D C [ C max ] [ C ] = λ D H [ H max ] [ H ] = λ D E [ E max ] [ E ] , δ D C [ C max ] [ C ] = δ D ,
10 λ M I L 10 [ I L 10 max ] [ I L 10 ] = 10 λ M I γ [ I γ max ] [ I γ ] = 10 λ M I L 12 [ I L 12 max ] [ I L 12 ] = λ M T h [ T h max ] [ T h ] = 10 λ M E [ E max ] [ E ] ,
50 δ C T c [ T c max ] [ T c ] = δ C I γ [ I γ max ] [ I γ ] ,
λ H D [ D max ] [ D ] = λ H N [ N max ] [ N ] = 10 λ H M [ M max ] [ M ] = λ H T c [ T c max ] [ T c ] = λ H C [ C max ] [ C ] ,
λ I L 12 M [ M max ] [ M ] = λ I L 12 D [ D max ] [ D ] = λ I L 12 T h [ T h max ] [ T h ] = λ I L 12 T c [ T c max ] [ T c ] ,
λ I L 10 M [ M max ] [ M ] = λ I L 10 D [ D max ] [ D ] = λ I L 10 T r [ T r max ] [ T r ] = λ I L 10 T h [ T h max ] [ T h ] = λ I L 10 T c [ T c max ] [ T c ] = λ I L 10 C [ C max ] [ C ] ,
λ I γ T c [ T c max ] [ T c ] = 5 λ I γ T h [ T h max ] [ T h ] = 5 λ I γ D [ D max ] [ D ] ,
λ I L 6 A [ A max ] [ A ] = λ I L 6 M [ M max ] [ M ] = λ I L 6 D [ D max ] [ D ] ,
2 λ C A [ A max ] [ A ] = λ C I L 6 [ I L 6 max ] [ I L 6 ] , 20 λ E A [ A max ] [ A ] = λ E ,
δ C = δ C T c [ T c max ] [ T c ] = 6 δ C I γ [ I γ max ] [ I γ ] .

Appendix A.3. Parameter Values

This section includes all the parameter values used in this project. Values are separated into two groups. First, we have degradation rates included from prior research (see Table A3). Second, we have scaling-dependent parameters (see Table 1) based on the assumptions made in Appendix A.1. In order to sufficiently analyze the dynamics of the tumor based on the estimated parameters, we emphasize that possible variations in assumptions and patient to patient values should be considered. These variations have been taken into account and addressed by performing sensitivity analysis and clustering patients based on their immune profiles. In their dimensional forms, the scaling-dependent parameters would also depend on the scaling constant α , in addition to derivation assumptions and patient data. Therefore, we list the objective non-dimensional values. The dimension of these parameters is per day.
Table A3. Half-lives and estimated death rate. Degredation and death rate taken or calculated from the given references.
Table A3. Half-lives and estimated death rate. Degredation and death rate taken or calculated from the given references.
ParameterValueReferenceParameterValueReference
δ T N 9.49 × 10 4 [71] δ H 18[71,121]
δ T c 0.231[71] δ E 4.16[122]
δ T h 0.406[71] δ I L 6 1.07[71]
δ T r 0.406[71] δ I L 10 4.62[71]
δ D 0.277[71] δ I L 12 128[123]
δ M 0.0198[71] δ I γ 33.3[71]
δ A 2.8 × 10 3 [125]
Table A4. Non-dimensional parameter values for each cluster.
Table A4. Non-dimensional parameter values for each cluster.
ParameterCluster 1Cluster 2Cluster 3Cluster 4Cluster 5
λ ¯ T h H 1.791 × 10 1 4.048 × 10 2 1.002 × 10 1 2.478 × 10 1 3.871 × 10 1
λ ¯ T h D 2.933 1.801 2.200 4.027 3.135
λ ¯ T h I L 12 1.569 × 10 1 2.286 × 10 2 6.850 × 10 2 2.623 × 10 1 3.919 × 10 1
λ ¯ T h E 1.756 × 10 1 2.633 × 10 2 1.063 × 10 1 2.795 × 10 1 4.696 × 10 1
λ ¯ T c E 6.022 3.297 4.332 8.420 7.669
λ ¯ T c D 5.028 × 10 3 1.128 × 10 2 4.482 × 10 3 6.064 × 10 3 2.560 × 10 3
λ ¯ T c I L 12 2.690 × 10 2 1.431 × 10 2 1.396 × 10 2 3.950 × 10 2 3.200 × 10 2
λ ¯ T r D 5.783 × 10 2 1.334 × 10 1 6.763 × 10 2 5.167 × 10 2 2.721 × 10 2
λ ¯ T r E 1.732 × 10 1 9.755 × 10 2 1.634 × 10 1 1.793 × 10 1 2.038 × 10 1
λ ¯ D C 3.346 × 10 2 6.768 × 10 2 4.237 × 10 2 2.540 × 10 2 2.505 × 10 2
λ ¯ D H 1.526 × 10 1 1.709 × 10 1 1.444 × 10 1 1.474 × 10 1 1.430 × 10 1
λ ¯ D E 1.497 × 10 1 1.112 × 10 1 1.532 × 10 1 1.663 × 10 1 1.734 × 10 1
λ ¯ M I L 10 2.499 × 10 3 1.201 × 10 3 1.978 × 10 3 2.535 × 10 3 3.128 × 10 3
λ ¯ M I γ 1.720 × 10 3 7.178 × 10 4 1.244 × 10 3 1.941 × 10 3 2.177 × 10 3
λ ¯ M I L 12 1.991 × 10 3 9.880 × 10 4 1.469 × 10 3 2.050 × 10 3 2.484 × 10 3
λ ¯ M T h 1.136 × 10 2 1.575 × 10 2 1.283 × 10 2 1.109 × 10 2 9.035 × 10 3
λ ¯ M E 2.229 × 10 3 1.138 × 10 3 2.279 × 10 3 2.184 × 10 3 2.977 × 10 3
λ ¯ C 5.756 × 10 2 4.250 × 10 2 5.669 × 10 2 3.427 × 10 2 3.959 × 10 2
λ ¯ C I L 6 4.345 × 10 4 8.494 × 10 4 3.381 × 10 4 5.530 × 10 3 4.595 × 10 3
λ ¯ C A 2.357 × 10 4 1.243 × 10 3 2.859 × 10 4 1.807 × 10 3 1.205 × 10 3
λ ¯ A 3.384 × 10 3 3.394 × 10 3 3.367 × 10 3 3.400 × 10 3 3.282 × 10 3
λ ¯ H D 1.458 2.113 1.529 1.239 7.319 × 10 1
λ ¯ H N 4.156 2.858 2.861 3.552 4.083
λ ¯ H M 9.720 × 10 1 8.862 × 10 1 1.421 7.101 × 10 1 7.883 × 10 1
λ ¯ H T c 3.605 4.620 4.014 7.244 6.063
λ ¯ H C 7.809 7.523 8.173 5.255 6.334
λ ¯ I L 12 M 5.253 × 10 1 4.698 × 10 1 6.480 × 10 1 3.738 × 10 1 4.730 × 10 1
λ ¯ I L 12 D 7.880 1.120 × 10 1 6.972 6.521 4.391
λ ¯ I L 12 T h 4.810 × 10 1 4.533 × 10 1 3.793 × 10 1 4.596 × 10 1 3.994 × 10 1
λ ¯ I L 12 T c 1.948 × 10 1 2.449 × 10 1 1.830 × 10 1 3.813 × 10 1 3.638 × 10 1
λ ¯ I L 10 M 1.197 1.155 1.606 9.978 × 10 1 1.252
λ ¯ I L 10 D 1.795 × 10 1 2.754 × 10 1 1.728 × 10 1 1.741 × 10 1 1.162 × 10 1
λ ¯ I L 10 T r 7.419 × 10 1 4.921 × 10 1 5.241 × 10 1 4.651 × 10 1 2.271 × 10 1
λ ¯ I L 10 T h 1.096 1.115 9.400 × 10 1 1.227 1.057
λ ¯ I L 10 T c 4.440 × 10 1 6.022 × 10 1 4.536 × 10 1 1.018 9.626 × 10 1
λ ¯ I L 10 C 9.615 × 10 1 9.806 × 10 1 9.235 × 10 1 7.384 × 10 1 1.006
λ ¯ E A 1.024 × 10 1 8.334 × 10 2 9.433 × 10 2 1.389 × 10 1 1.278 × 10 1
λ ¯ E 5.935 4.761 5.600 7.870 8.708
λ ¯ I γ T c 2.115 × 10 1 2.278 × 10 1 2.234 × 10 1 2.611 × 10 1 2.677 × 10 1
λ ¯ I γ T h 1.044 × 10 1 8.434 9.259 6.295 5.879
λ ¯ I γ D 1.711 2.084 1.702 8.931 × 10 1 6.465 × 10 1
λ ¯ I L 6 A 5.692 × 10 1 5.110 × 10 1 4.652 × 10 1 5.329 × 10 1 5.150 × 10 1
λ ¯ I L 6 M 4.355 × 10 1 4.513 × 10 1 5.460 × 10 1 4.573 × 10 1 5.079 × 10 1
λ ¯ I L 6 D 6.532 × 10 2 1.076 × 10 1 5.875 × 10 2 7.977 × 10 2 4.715 × 10 2
δ ¯ T h T r 7.562 × 10 1 6.086 × 10 1 5.961 × 10 1 6.522 × 10 1 2.427 × 10 1
δ ¯ T h I L 10 2.457 1.051 1.648 3.933 3.910
δ ¯ T c I L 10 4.319 1.847 2.896 6.913 6.871
δ ¯ T c T r 1.329 1.070 1.048 1.146 4.266 × 10 1
δ ¯ D C 5.876 × 10 2 7.271 × 10 2 6.297 × 10 2 6.209 × 10 2 6.444 × 10 2
δ ¯ D N 2.770 × 10 1 2.770 × 10 1 2.770 × 10 1 2.770 × 10 1 2.770 × 10 1
δ ¯ M N 1.980 × 10 2 1.980 × 10 2 1.980 × 10 2 1.980 × 10 2 1.980 × 10 2
δ ¯ C T c 4.399 × 10 3 5.276 × 10 3 4.936 × 10 3 8.052 × 10 3 6.762 × 10 3
δ ¯ C I γ 2.741 × 10 3 7.415 × 10 4 1.654 × 10 3 2.832 × 10 3 2.982 × 10 3
δ ¯ C 4.492 × 10 2 3.273 × 10 2 4.421 × 10 2 2.606 × 10 2 3.037 × 10 2
δ ¯ N 2.043 × 10 1 2.130 × 10 1 3.031 × 10 1 1.141 × 10 1 1.300 × 10 1
A ¯ T N 9.730 5.445 7.057 1.351 × 10 1 1.232 × 10 1
A ¯ D N 6.128 × 10 1 6.267 × 10 1 6.170 × 10 1 6.161 × 10 1 6.184 × 10 1
A ¯ M 3.960 × 10 2 3.960 × 10 2 3.960 × 10 2 3.960 × 10 2 3.960 · 10 2
α ¯ N C 3.924 5.497 5.966 3.089 3.240
C ¯ 0 9.428 7.619 8.798 8.923 8.597
A ¯ 0 5.795 5.713 5.937 5.666 6.813
E ¯ 0 3.161 6.958 3.649 2.045 1.862

Appendix A.4. Dynamics with Varying Initial Conditions

Figure A1. The dynamics with varying initial conditions. Subfigures (AE) show the dynamics of cells and cytokines with initial conditions from different patients in clusters 1, 2, 3, 4, and 5, respectively.
Figure A1. The dynamics with varying initial conditions. Subfigures (AE) show the dynamics of cells and cytokines with initial conditions from different patients in clusters 1, 2, 3, 4, and 5, respectively.
Jpm 11 01031 g0a1aJpm 11 01031 g0a1bJpm 11 01031 g0a1c

Appendix A.5. Dynamics of the Tumor Microenvironment with Cross-Cluster Initial Conditions

Figure A2. Dynamics with cross-cluster initial conditions: (A) parameters from cluster 1 and initial conditions from clusters 2, 3, 4, and 5; (B) parameters from cluster 2 and initial conditions from clusters 1, 3, 4, and 5; (C) parameters from cluster 3 and initial conditions from clusters 1, 2, 4, and 5; (D) parameters from cluster 4 and initial conditions from clusters 1, 2, 3, and 5; and (E) parameters from cluster 5 and initial conditions from clusters 1, 2, 3, and 4.
Figure A2. Dynamics with cross-cluster initial conditions: (A) parameters from cluster 1 and initial conditions from clusters 2, 3, 4, and 5; (B) parameters from cluster 2 and initial conditions from clusters 1, 3, 4, and 5; (C) parameters from cluster 3 and initial conditions from clusters 1, 2, 4, and 5; (D) parameters from cluster 4 and initial conditions from clusters 1, 2, 3, and 5; and (E) parameters from cluster 5 and initial conditions from clusters 1, 2, 3, and 4.
Jpm 11 01031 g0a2aJpm 11 01031 g0a2bJpm 11 01031 g0a2c

Appendix A.6. Bifurcation and Lyapunov Exponent for the Cancer ODE

To further investigate the effects of parameters on controlling the cancer development, we carried out bifurcation analysis on those involved in cancer ODE (9). This choice was made not only because these parameters were directly involved in cancer ODE, but also since they came out as the most sensitive parameters in Section 3.3. Our results show that regardless of minute differences in quantities, the qualitative behaviors are identical across all the clusters; see Figure A3.
In all clusters, for λ C , λ C I L 6 , λ C A , δ C T c , δ C I γ , and δ C , we can see the existence of a single equilibrium for a large range of parameter values (0 to 1). Increasing the values of λ C , λ C I L 6 , and λ C A caused the cancer population to reach higher steady state values. However, for λ C in clusters 1–5 and for λ C A in cluster 5, we observed interesting results at the beginning of their corresponding intervals. More specifically, very small values of λ C caused the cancer population to vanish for all clusters, whereas very small values of λ C A led to very small steady state (but not zero) values for cancer only in cluster 5; see Figure A3. On the contrary, parameters δ C T c , δ C I γ and δ C only pertain to a non-vanishing cancer population in a very small parameter space. For the most part, they caused the cancer population to disappear at the steady state; see Figure A3. The negative Lyapunov exponents show that all the steady state values are Lyapunov stable.
Figure A3. Bifurcation and Lyapunov exponent diagrams for cancer parameters: Subfigures show the bifurcation on top of the corresponding Lyapunov diagrams for clusters 1–5. Bifurcation was done for parameters λ C , λ C I L 6 , λ C A , δ C T c , δ C I γ , and δ C from the cancer ODE (9).
Figure A3. Bifurcation and Lyapunov exponent diagrams for cancer parameters: Subfigures show the bifurcation on top of the corresponding Lyapunov diagrams for clusters 1–5. Bifurcation was done for parameters λ C , λ C I L 6 , λ C A , δ C T c , δ C I γ , and δ C from the cancer ODE (9).
Jpm 11 01031 g0a3

Appendix A.7. Positivity

In order to prove that the system with positive parameters and positive initial conditions has positive solutions, we defined the following integrating factors for each of the state variables.
η T h ( t ) = exp 0 t δ T h T r [ T r ] + δ T h I L 10 [ I L 10 ] + δ T h d s , η T c ( t ) = exp 0 t δ T c T r [ T r ] + δ T c I L 10 [ I L 10 ] + δ T c d s , η T r ( t ) = exp δ T r t , η T N ( t ) = exp 0 t ( λ T h H [ H ] + λ T h D [ D ] + λ T h I L 12 [ I L 12 ] + λ T h E [ E ] + λ T c E [ E ] + λ T c D [ D ] + λ T c I L 12 [ I L 12 ] + λ T r D [ D ] + λ T r E [ E ] + δ T N ) d s ] , η D N ( t ) = exp 0 t λ D C [ C ] + λ D H [ H ] + λ D E [ E ] + δ D N d s , η D ( t ) = exp 0 t δ D C [ C ] + δ D d s , η M N ( t ) = exp 0 t λ M I L 10 [ I L 10 ] + λ M I γ [ I γ ] + λ M I L 12 [ I L 12 ] + λ M T h [ T h ] + λ M E [ E ] + δ M N d s , η M ( t ) = exp δ M t , η C ( t ) = exp 0 t δ C T c [ T c ] + δ C I γ [ I γ ] + δ C λ C + λ C I L 6 [ I L 6 ] + λ C A [ A ] 1 [ C ] C 0 d s , η N ( t ) = exp δ N t , η H ( t ) = exp δ H t , η I L 12 ( t ) = exp δ I L 12 t , η I L 10 ( t ) = exp δ I L 10 t , η E ( t ) = exp 0 t δ E λ E 1 [ E ] E 0 d s , η I γ ( t ) = exp δ I γ t , η I L 6 = exp δ I L 6 t .
Considering these integration factors, we rewrote our system of ODEs as:
d d t [ T h ] η T h = λ T h H [ H ] + λ T h D [ D ] + λ T h I L 12 [ I L 12 ] + λ T h E [ E ] [ T N ] η T h , d d t [ T c ] η T c = λ T c E [ E ] + λ T c D [ D ] + λ T c I L 12 [ I L 12 ] [ T N ] η T c , d d t [ T r ] η T r = λ T r D [ D ] + λ T r E [ E ] [ T N ] η T r , d d t [ T N ] η T N = A T N η T N , d d t [ D N ] η D N = A D N η D n , d d t [ D ] η D = λ D C [ C ] + λ D H [ H ] + λ D E [ E ] [ D N ] η D , d d t [ M N ] η M N = A M η M , d d t [ M ] η M = λ M I L 10 [ I L 10 ] + λ M I γ [ I γ ] + λ M I L 12 [ I L 12 ] + λ M T h [ T h ] + λ M E [ E ] [ M N ] η M , d d t [ C ] η C = 0 , d d t [ N ] η N = α N C δ C I γ [ I γ ] + δ C T c [ T c ] + δ C [ C ] η N , d d t [ H ] η H = λ H D [ D ] + λ H N [ N ] + λ H M [ M ] + λ H T c [ T c ] + λ H C [ C ] η H , d d t [ I L 12 ] η I L 12 = λ I L 12 M [ M ] + λ I L 12 D [ D ] + λ I L 12 T h [ T h ] + λ I L 12 T c [ T c ] η I L 12 , d d t [ I L 10 ] η I L 10 = λ I L 10 M [ M ] + λ I L 10 D [ D ] + λ I L 10 T h [ T r ] + λ I L 10 T h [ T h ] + λ I L 10 T c [ T c ] + λ I L 10 C [ C ] η I L 10 , d d t [ E ] η E = λ E A [ A ] η E d d t [ I γ ] η I γ = λ I γ T c [ T c ] + λ I γ T h [ T h ] + λ I γ D [ E ] [ D ] η I γ , d d t [ I L 6 ] η I L 6 = λ I L 6 A [ A ] + λ I L 6 M [ M ] + λ I L 6 D [ D ] η I L 6 .
The right-hand sides of all of the above equations are non-negative. This means that the terms [ X ] η X are non-decreasing for each state variable [ X ] . Therefore, for positive initial values, the dynamics stayed positive.

References

  1. About Breast Cancer. Available online: https://www.cancer.org/cancer/breast-cancer/about/how-common-is-breast-cancer.html (accessed on 20 May 2021).
  2. Harbeck, N.; Penault-Llorca, F.; Cortes, J.; Gnant, M.; Houssami, N.; Poortmans, P.; Ruddy, K.; Tsang, J.; Cardoso, F. Breast cancer. Nat. Rev. Dis. Prim. 2019, 5, 66. [Google Scholar] [CrossRef] [PubMed]
  3. Engstrøm, M.J.; Opdahl, S.; Hagen, A.I.; Romundstad, P.R.; Akslen, L.A.; Haugen, O.A.; Vatten, L.J.; Bofin, A.M. Molecular subtypes, histopathological grade and survival in a historic cohort of breast cancer patients. Breast Cancer Res. Treat. 2013, 140, 463–473. [Google Scholar] [CrossRef] [Green Version]
  4. Venkitaraman, R. Triple-negative/basal-like breast cancer: Clinical, pathologic and molecular features. Expert Rev. Anticancer Ther. 2010, 10, 199–207. [Google Scholar] [CrossRef]
  5. Bai, X.; Ni, J.; Beretov, J.; Graham, P.; Li, Y. Cancer stem cell in breast cancer therapeutic resistance. Cancer Treat. Rev. 2018, 69, 152–163. [Google Scholar] [CrossRef]
  6. Masoud, V.; Pagès, G. Targeted therapies in breast cancer: New challenges to fight against resistance. World J. Clin. Oncol. 2017, 8, 120. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Vagia, E.; Mahalingam, D.; Cristofanilli, M. The Landscape of Targeted Therapies in TNBC. Cancers 2020, 12, 916. [Google Scholar] [CrossRef] [Green Version]
  8. Nedeljković, M.; Damjanović, A. Mechanisms of Chemotherapy Resistance in Triple-Negative Breast Cancer—How We Can Rise to the Challenge. Cells 2019, 8, 957. [Google Scholar] [CrossRef] [Green Version]
  9. Kim, C.; Gao, R.; Sei, E.; Brandt, R.; Hartman, J.; Hatschek, T.; Crosetto, N.; Foukakis, T.; Navin, N.E. Chemoresistance Evolution in Triple-Negative Breast Cancer Delineated by Single-Cell Sequencing. Cell 2018, 173, 879–893.e13. [Google Scholar] [CrossRef] [Green Version]
  10. Place, A.E.; Jin Huh, S.; Polyak, K. The microenvironment in breast cancer progression: Biology and implications for treatment. Breast Cancer Res. 2011, 13, 1–11. [Google Scholar] [CrossRef] [Green Version]
  11. Korkaya, H.; Liu, S.; Wicha, M.S. Breast cancer stem cells, cytokine networks, and the tumor microenvironment. J. Clin. Investig. 2011, 121, 3804–3809. [Google Scholar] [CrossRef] [PubMed]
  12. Soysal, S.D.; Tzankov, A.; Muenst, S.E. Role of the Tumor Microenvironment in Breast Cancer. Pathobiology 2015, 82, 142–152. [Google Scholar] [CrossRef]
  13. Mittal, S.; Brown, N.J.; Holen, I. The breast tumor microenvironment: Role in cancer development, progression and response to therapy. Expert Rev. Mol. Diagn. 2018, 18, 227–243. [Google Scholar] [CrossRef] [PubMed]
  14. Wu, T.; Dai, Y. Tumor microenvironment and therapeutic response. Cancer Lett. 2017, 387, 61–68. [Google Scholar] [CrossRef] [PubMed]
  15. Palucka, K.; Banchereau, J. Cancer immunotherapy via dendritic cells. Nat. Rev. Cancer 2012, 12, 265–277. [Google Scholar] [CrossRef]
  16. Palucka, K.; Coussens, L.M.; O’Shaughnessy, J. Dendritic cells, inflammation, and breast cancer. Cancer J. (Sudbury, Mass.) 2013, 19, 511–516. [Google Scholar] [CrossRef] [Green Version]
  17. Tran Janco, J.M.; Lamichhane, P.; Karyampudi, L.; Knutson, K.L. Tumor-Infiltrating Dendritic Cells in Cancer Pathogenesis. J. Immunol. 2015, 194, 2985–2991. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  18. Fu, C.; Jiang, A. Dendritic Cells and CD8 T Cell Immunity in Tumor Microenvironment. Front. Immunol. 2018, 9, 3059. [Google Scholar] [CrossRef] [Green Version]
  19. Segovia-Mendoza, M.; Morales-Montor, J. Immune Tumor Microenvironment in Breast Cancer and the Participation of Estrogen and Its Receptors in Cancer Physiopathology. Front. Immunol. 2019, 10, 348. [Google Scholar] [CrossRef] [Green Version]
  20. Melssen, M.; Slingluff, C.L. Vaccines Targeting Helper T Cells for Cancer Immunotherapy. Curr. Opin. Immunol. 2017, 47, 85–92. [Google Scholar] [CrossRef] [PubMed]
  21. Ni, L.; Lu, J. Interferon gamma in cancer immunotherapy. Cancer Med. 2018, 7, 4509–4516. [Google Scholar] [CrossRef]
  22. Moon, B.I.; Kim, T.H.; Seoh, J.Y. Functional Modulation of Regulatory T Cells by IL-2. PLoS ONE 2015, 10, e0141864. [Google Scholar] [CrossRef]
  23. Lai, X.; Stiff, A.; Duggan, M.; Wesolowski, R.; Carson, W.E.; Friedman, A. Modeling combination therapy for breast cancer with BET and immune checkpoint inhibitors. Proc. Natl. Acad. Sci. USA 2018, 115, 5534–5539. [Google Scholar] [CrossRef] [Green Version]
  24. Obeid, E.; Nanda, R.; Fu, Y.X.; Olopade, O.I. The role of tumor-associated macrophages in breast cancer progression (review). Int. J. Oncol. 2013, 43, 5–12. [Google Scholar] [CrossRef] [Green Version]
  25. Bingle, L.; Brown, N.J.; Lewis, C.E. The role of tumour-associated macrophages in tumour progression: Implications for new anticancer therapies. J. Pathol. 2002, 196, 254–265. [Google Scholar] [CrossRef] [PubMed]
  26. Ali, H.; Provenzano, E.; Dawson, S.J.; Blows, F.; Liu, B.; Shah, M.; Earl, H.; Poole, C.; Hiller, L.; Dunn, J.; et al. Association between CD8+ T-cell infiltration and breast cancer survival in 12 439 patients. Ann. Oncol. 2014, 25, 1536–1543. [Google Scholar] [CrossRef]
  27. Ali, H.R.; Chlon, L.; Pharoah, P.D.P.; Markowetz, F.; Caldas, C. Patterns of Immune Infiltration in Breast Cancer and Their Clinical Implications: A Gene-Expression-Based Retrospective Study. PLoS Med. 2016, 13, e1002194. [Google Scholar] [CrossRef]
  28. Wang, Z.K.; Yang, B.; Liu, H.; Hu, Y.; Yang, J.L.; Wu, L.L.; Zhou, Z.H.; Jiao, S.C. Regulatory T cells increase in breast cancer and in stage IV breast cancer. Cancer Immunol. Immunother. 2012, 61, 911–916. [Google Scholar] [CrossRef]
  29. Shahriyari, L.; Komarova, N.L. Symmetric vs. asymmetric stem cell divisions: An adaptation against cancer? PLoS ONE 2013, 8, e76195. [Google Scholar] [CrossRef] [Green Version]
  30. Butner, J.D.; Wang, Z.; Elganainy, D.; Al Feghali, K.A.; Plodinec, M.; Calin, G.A.; Dogra, P.; Nizzero, S.; Ruiz-Ramírez, J.; Martin, G.V.; et al. A mathematical model for the quantification of a patient’s sensitivity to checkpoint inhibitors and long-term tumour burden. Nat. Biomed. Eng. 2021, 5, 297–308. [Google Scholar] [CrossRef] [PubMed]
  31. Shahriyari, L.; Komarova, N.L. The role of the bi-compartmental stem cell niche in delaying cancer. Phys. Biol. 2015, 12, 055001. [Google Scholar] [CrossRef] [PubMed]
  32. Shahriyari, L. Cell dynamics in tumour environment after treatments. J. R. Soc. Interface 2017, 14, 20160977. [Google Scholar] [CrossRef] [PubMed]
  33. Chamseddine, I.M.; Rejniak, K.A. Hybrid modeling frameworks of tumor development and treatment. Wiley Interdiscip. Rev. Syst. Biol. Med. 2019, 12, e1461. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. McKenna, M.T.; Weis, J.A.; Barnes, S.L.; Tyson, D.R.; Miga, M.I.; Quaranta, V.; Yankeelov, T.E. A Predictive Mathematical Modeling Approach for the Study of Doxorubicin Treatment in Triple Negative Breast Cancer. Sci. Rep. 2017, 7, 1–14. [Google Scholar] [CrossRef] [Green Version]
  35. Shahriyari, L.; Mahdipour-Shirayeh, A. Modeling dynamics of mutants in heterogeneous stem cell niche. Phys. Biol. 2017, 14, 016004. [Google Scholar] [CrossRef] [PubMed]
  36. Mehdizadeh, R.; Shariatpanahi, S.P.; Goliaei, B.; Peyvandi, S.; Rüegg, C. Dormant tumor cell vaccination: A mathematical model of immunological dormancy in triple-negative breast cancer. Cancers 2021, 13, 245. [Google Scholar] [CrossRef]
  37. Budithi, A.; Su, S.; Kirshtein, A.; Shahriyari, L. Data driven mathematical model of FOLFIRI treatment for colon cancer. Cancers 2021, 13, 2632. [Google Scholar] [CrossRef]
  38. Solís-Pérez, J.E.; Gómez-Aguilar, J.F.; Atangana, A. A fractional mathematical model of breast cancer competition model. Chaos Solitons Fractals 2019, 127, 38–54. [Google Scholar] [CrossRef]
  39. Wodarz, D.; Komarova, N. Dynamics of Cancer: Mathematical Foundations of Oncology; World Scientific Publishing Company: Singapore, 2014. [Google Scholar]
  40. Chakrabarti, A.; Verbridge, S.; Stroock, A.D.; Fischbach, C.; Varner, J.D. Multiscale models of breast cancer progression. Ann. Biomed. Eng. 2012, 40, 2488–2500. [Google Scholar] [CrossRef]
  41. Abernathy, K.; Abernathy, Z.; Baxter, A.; Stevens, M. Global dynamics of a breast cancer competition model. Differ. Equ. Dyn. Syst. 2020, 28, 791–805. [Google Scholar] [CrossRef]
  42. Knútsdóttir, H.; Pálsson, E.; Edelstein-Keshet, L. Mathematical model of macrophage-facilitated breast cancer cells invasion. J. Theor. Biol. 2014, 357, 184–199. [Google Scholar] [CrossRef] [PubMed]
  43. Jarrett, A.M.; Bloom, M.J.; Godfrey, W.; Syed, A.K.; Ekrut, D.A.; Ehrlich, L.I.; Yankeelov, T.E.; Sorace, A.G. Mathematical modelling of trastuzumab-induced immune response in an in vivo murine model of HER2+ breast cancer. Math. Med. Biol. A J. IMA 2019, 36, 381–410. [Google Scholar] [CrossRef] [PubMed]
  44. Roe-Dale, R.; Isaacson, D.; Kupferschmid, M. A Mathematical Model of Breast Cancer Treatment with CMF and Doxorubicin. Bull. Math. Biol. 2011, 73, 585–608. [Google Scholar] [CrossRef]
  45. Wei, H.C. Mathematical modeling of tumor growth: The MCF-7 breast cancer cell line. Math. Biosci. Eng. 2019, 16, 6512–6535. [Google Scholar] [CrossRef]
  46. Johnston, M.D.; Edwards, C.M.; Bodmer, W.F.; Maini, P.K.; Chapman, S.J. Mathematical modeling of cell population dynamics in the colonic crypt and in colorectal cancer. Proc. Natl. Acad. Sci. USA 2007, 104, 4008–4013. [Google Scholar] [CrossRef] [Green Version]
  47. Lewin, T.D.; Byrne, H.M.; Maini, P.K.; Caudell, J.J.; Moros, E.G.; Enderling, H. The importance of dead material within a tumour on the dynamics in response to radiotherapy. Phys. Med. Biol. 2020, 65, 015007. [Google Scholar] [CrossRef]
  48. Lewin, T.D.; Maini, P.K.; Moros, E.G.; Enderling, H.; Byrne, H.M. A three phase model to investigate the effects of dead material on the growth of avascular tumours. Math. Model. Nat. Phenom. 2020, 15, 22. [Google Scholar] [CrossRef] [Green Version]
  49. Friedman, A.; Hao, W. The Role of Exosomes in Pancreatic Cancer Microenvironment. Bull. Math. Biol. 2018, 80, 1111–1133. [Google Scholar] [CrossRef] [PubMed]
  50. Hao, W.; Crouser, E.D.; Friedman, A. Mathematical model of sarcoidosis. Proc. Natl. Acad. Sci. USA 2014, 111, 16065–16070. [Google Scholar] [CrossRef] [Green Version]
  51. Enderling, H.; Chaplain, M.A.J.; Anderson, A.R.A.; Vaidya, J.S. A mathematical model of breast cancer development, local treatment and recurrence. J. Theor. Biol. 2007, 246, 245–259. [Google Scholar] [CrossRef]
  52. Anderson, A.R.; Chaplain, M.A.; Newman, E.L.; Steele, R.J.; Thompson, A.M. Mathematical modelling of tumour invasion and metastasis. J. Theor. Med. 2000, 2, 129–154. [Google Scholar] [CrossRef] [Green Version]
  53. Dumitriu, I.E.; Baruah, P.; Valentinis, B.; Voll, R.E.; Herrmann, M.; Nawroth, P.P.; Arnold, B.; Bianchi, M.E.; Manfredi, A.A.; Rovere-Querini, P. Release of High Mobility Group Box 1 by Dendritic Cells Controls T Cell Activation via the Receptor for Advanced Glycation End Products. J. Immunol. 2005, 174, 7506–7515. [Google Scholar] [CrossRef] [Green Version]
  54. Bell, R.B.; Feng, Z.; Bifulco, C.B.; Leidner, R.; Weinberg, A.; Fox, B.A. 15-Immunotherapy. In Oral, Head and Neck Oncology and Reconstructive Surgery; Elsevier: Amsterdam, The Netherlands, 2018; pp. 314–340. [Google Scholar]
  55. Wang, K.; Vella, A.T. Regulatory T Cells and Cancer: A Two-Sided Story. Immunol. Investig. 2016, 45, 797–812. [Google Scholar] [CrossRef]
  56. Mohammad, I.; Starskaia, I.; Nagy, T.; Guo, J.; Yatkin, E.; Väänänen, K.; Watford, W.T.; Chen, Z. Estrogen receptor α contributes to T cell mediated autoimmune inflammation by promoting T cell activation and proliferation. Sci. Signal. 2018, 11. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  57. Ma, Y.; Shurin, G.V.; Peiyuan, Z.; Shurin, M.R. Dendritic cells in the cancer microenvironment. J. Cancer 2013, 4, 36–44. [Google Scholar] [CrossRef] [Green Version]
  58. Luo, C.Y.; Wang, L.; Sun, C.; Li, D.J. Estrogen enhances the functions of CD4+CD25+Foxp3+ regulatory T cells that suppress osteoclast differentiation and bone resorption in vitro. Cell. Mol. Immunol. 2011, 8, 50–58. [Google Scholar] [CrossRef] [PubMed]
  59. Tang, D.; Kang, R.; Zeh, H.J.; Lotze, M.T. High-mobility Group Box 1 [HMGB1] and Cancer. Biochim. Biophys. Acta 2010, 1799, 131–140. [Google Scholar] [CrossRef] [Green Version]
  60. Lee, S.Y.; Ju, M.K.; Jeon, H.M.; Jeong, E.K.; Lee, Y.J.; Kim, C.H.; Park, H.G.; Han, S.I.; Kang, H.S. Regulation of Tumor Progression by Programmed Necrosis. Oxidative Med. Cell. Longev. 2018, 2018, 3537471. [Google Scholar] [CrossRef] [Green Version]
  61. Apetoh, L.; Ghiringhelli, F.; Tesniere, A.; Criollo, A.; Ortiz, C.; Lidereau, R.; Mariette, C.; Chaput, N.; Mira, J.P.; Delaloge, S.; et al. The interaction between HMGB1 and TLR4 dictates the outcome of anticancer chemotherapy and radiotherapy. Immunol. Rev. 2007, 220, 47–59. [Google Scholar] [CrossRef]
  62. Aras, S.; Zaidi, M.R. TAMeless traitors: Macrophages in cancer progression and metastasis. Br. J. Cancer 2017, 117, 1583–1591. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  63. Sica, A.; Saccani, A.; Mantovani, A. Tumor-associated macrophages: A molecular perspective. Int. Immunopharmacol. 2002, 2, 1045–1054. [Google Scholar] [CrossRef]
  64. Tariq, M.; Zhang, J.; Liang, G.; Ding, L.; He, Q.; Yang, B. Macrophage Polarization: Anti-Cancer Strategies to Target Tumor-Associated Macrophage in Breast Cancer. J. Cell. Biochem. 2017, 118, 2484–2501. [Google Scholar] [CrossRef] [PubMed]
  65. Chanmee, T.; Ontong, P.; Konno, K.; Itano, N. Tumor-Associated Macrophages as Major Players in the Tumor Microenvironment. Cancers 2014, 6, 1670–1690. [Google Scholar] [CrossRef] [Green Version]
  66. Wu, Q.; Li, B.; Li, Z.; Li, J.; Sun, S.; Sun, S. Cancer-associated adipocytes: Key players in breast cancer progression. J. Hematol. Oncol. 2019, 12, 1–15. [Google Scholar] [CrossRef] [PubMed]
  67. Mao, Y.; Keller, E.T.; Garfield, D.H.; Shen, K.; Wang, J. Stroma Cells in Tumor Microenvironment and Breast Cancer. Cancer Metastasis Rev. 2013, 32, 303–315. [Google Scholar] [CrossRef] [Green Version]
  68. Chu, D.T.; Phuong, T.N.T.; Tien, N.L.B.; Tran, D.K.; Nguyen, T.T.; Thanh, V.V.; Quang, T.L.; Minh, L.B.; Pham, V.H.; Ngoc, V.T.N.; et al. The Effects of Adipocytes on the Regulation of Breast Cancer in the Tumor Microenvironment: An Update. Cells 2019, 8, 857. [Google Scholar] [CrossRef] [Green Version]
  69. Neel, J.C.; Humbert, L.; Lebrun, J.J. The dual role of TGFβ in human cancer: From tumor suppression to cancer metastasis. Int. Sch. Res. Not. 2012, 2012, 381428. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  70. Wang, Y.Y.; Attané, C.; Milhas, D.; Dirat, B.; Dauvillier, S.; Guerard, A.; Gilhodes, J.; Lazar, I.; Alet, N.; Laurent, V.; et al. Mammary adipocytes stimulate breast cancer invasion through metabolic remodeling of tumor cells. JCI Insight 2017, 2, e87489. [Google Scholar] [CrossRef] [Green Version]
  71. Kirshtein, A.; Akbarinejad, S.; Hao, W.; Le, T.; Su, S.; Aronow, R.A.; Shahriyari, L. Data Driven Mathematical Model of Colon Cancer Progression. J. Clin. Med. 2020, 9, 3947. [Google Scholar] [CrossRef]
  72. Rubartelli, A.; Lotze, M.T. Inside, outside, upside down: Damage-associated molecular-pattern molecules (DAMPs) and redox. Trends Immunol. 2007, 28, 429–436. [Google Scholar] [CrossRef]
  73. Li, G.; Liang, X.; Lotze, M.T. HMGB1: The Central Cytokine for All Lymphoid Cells. Front. Immunol. 2013, 4, 68. [Google Scholar] [CrossRef] [Green Version]
  74. Wang, S.; Zhang, Y. HMGB1 in inflammation and cancer. J. Hematol. Oncol. 2020, 13, 116. [Google Scholar] [CrossRef]
  75. Bonaldi, T.; Talamo, F.; Scaffidi, P.; Ferrera, D.; Porto, A.; Bachi, A.; Rubartelli, A.; Agresti, A.; Bianchi, M.E. Monocytic cells hyperacetylate chromatin protein HMGB1 to redirect it towards secretion. EMBO J. 2003, 22, 5551–5560. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  76. Tang, D.; Shi, Y.; Kang, R.; Li, T.; Xiao, W.; Wang, H.; Xiao, X. Hydrogen peroxide stimulates macrophages and monocytes to actively release HMGB1. J. Leukoc. Biol. 2007, 81, 741–747. [Google Scholar] [CrossRef]
  77. Semino, C.; Angelini, G.; Poggi, A.; Rubartelli, A. NK/iDC interaction results in IL-18 secretion by DCs at the synaptic cleft followed by NK cell activation and release of the DC maturation factor HMGB1. Blood 2005, 106, 609–616. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  78. Gougeon, M.L.; Bras, M. Natural killer cells, dendritic cells, and the alarmin high-mobility group box 1 protein: A dangerous trio in HIV-1 infection? Curr. Opin. HIV AIDS 2011, 6, 364–372. [Google Scholar] [CrossRef]
  79. Fan, X.; Zhang, H.; Cheng, Y.; Jiang, X.; Zhu, J.; Jin, T. Double roles of macrophages in human neuroimmune diseases and their animal models. Mediat. Inflamm. 2016, 2016, 8489251. [Google Scholar] [CrossRef] [Green Version]
  80. Hart, A.L.; Al-Hassi, H.O.; Rigby, R.J.; Bell, S.J.; Emmanuel, A.V.; Knight, S.C.; Kamm, M.A.; Stagg, A.J. Characteristics of intestinal dendritic cells in inflammatory bowel diseases. Gastroenterology 2005, 129, 50–65. [Google Scholar] [CrossRef]
  81. Iwasaki, A.; Kelsall, B.L. Freshly isolated Peyer’s patch, but not spleen, dendritic cells produce interleukin 10 and induce the differentiation of T helper type 2 cells. J. Exp. Med. 1999, 190, 229–240. [Google Scholar] [CrossRef] [Green Version]
  82. Couper, K.N.; Blount, D.G.; Riley, E.M. IL-10: The master regulator of immunity to infection. J. Immunol. 2008, 180, 5771–5777. [Google Scholar] [CrossRef] [PubMed]
  83. Trinchieri, G. Interleukin-10 production by effector T cells: Th1 cells show self control. J. Exp. Med. 2007, 204, 239–243. [Google Scholar] [CrossRef]
  84. Mufudza, C.; Sorofa, W.; Chiyaka, E.T. Assessing the Effects of Estrogen on the Dynamics of Breast Cancer. Comput. Math. Methods Med. 2012, 2012, e473572. [Google Scholar] [CrossRef] [Green Version]
  85. Simpson, E.R. Sources of estrogen and their importance. J. Steroid Biochem. Mol. Biol. 2003, 86, 225–230. [Google Scholar] [CrossRef]
  86. Nakaya, M.; Tachibana, H.; Yamada, K. Effect of estrogens on the interferon-γ producing cell population of mouse splenocytes. Biosci. Biotechnol. Biochem. 2006, 70, 47–53. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  87. Liu, T.; Clark, R.K.; McDonnell, P.C.; Young, P.R.; White, R.F.; Barone, F.C.; Feuerstein, G.Z. Tumor necrosis factor-alpha expression in ischemic neurons. Stroke 1994, 25, 1481–1488. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  88. Newman, A.M.; Steen, C.B.; Liu, C.L.; Gentles, A.J.; Chaudhuri, A.A.; Scherer, F.; Khodadoust, M.S.; Esfahani, M.S.; Luca, B.A.; Steiner, D.; et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat. Biotechnol. 2019, 37, 773–782. [Google Scholar] [CrossRef] [PubMed]
  89. Le, T.; Aronow, R.A.; Kirshtein, A.; Shahriyari, L. A review of digital cytometry methods: Estimating the relative abundance of cell types in a bulk of cells. Brief. Bioinform. 2021, 22, bbaa219. [Google Scholar] [CrossRef]
  90. Su, S.; Akbarinejad, S.; Shahriyari, L. Immune classification of clear cell renal cell carcinoma. Sci. Rep. 2021, 11, 4338. [Google Scholar] [CrossRef] [PubMed]
  91. Le, T.; Su, S.; Shahriyari, L. Immune classification of osteosarcoma. Math. Biosci. Eng. 2021, 18, 1879–1897. [Google Scholar] [CrossRef]
  92. The Cancer Genome Atlas Network. Comprehensive molecular portraits of human breast tumours. Nature 2012, 490, 61–70. [Google Scholar] [CrossRef] [Green Version]
  93. Curtis, C.; Shah, S.P.; Chin, S.F.; Turashvili, G.; Rueda, O.M.; Dunning, M.J.; Speed, D.; Lynch, A.G.; Samarajiwa, S.; Yuan, Y.; et al. The genomic and transcriptomic architecture of 2000 breast tumours reveals novel subgroups. Nature 2012, 486, 346–352. [Google Scholar] [CrossRef]
  94. Cerami, E.; Gao, J.; Dogrusoz, U.; Gross, B.E.; Sumer, S.O.; Aksoy, B.A.; Jacobsen, A.; Byrne, C.J.; Heuer, M.L.; Larsson, E.; et al. The cBio Cancer Genomics Portal: An Open Platform for Exploring Multidimensional Cancer Genomics Data: Figure 1. Cancer Discov. 2012, 2, 401–404. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  95. Goldman, M.J.; Craft, B.; Hastie, M.; Repečka, K.; McDade, F.; Kamath, A.; Banerjee, A.; Luo, Y.; Rogers, D.; Brooks, A.N.; et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat. Biotechnol. 2020, 38, 675–678. [Google Scholar] [CrossRef] [PubMed]
  96. Kim, Y.; Friedman, A. Interaction of Tumor with Its Micro-environment: A Mathematical Model. Bull. Math. Biol. 2010, 72, 1029–1068. [Google Scholar] [CrossRef] [PubMed]
  97. Heiss, F.; Winschel, V. Likelihood approximation by numerical integration on sparse grids. J. Econom. 2008, 144, 62–80. [Google Scholar] [CrossRef] [Green Version]
  98. Gerstner, T.; Griebel, M. Numerical integration using sparse grids. Numer. Algorithms 1998, 18, 209–232. [Google Scholar] [CrossRef]
  99. Punt, J. Chapter 4—Adaptive Immunity: T Cells and Cytokines. In Cancer Immunotherapy, 2nd ed.; Prendergast, G.C., Jaffee, E.M., Eds.; Academic Press: San Diego, CA, USA, 2013; pp. 41–53. [Google Scholar]
  100. Qiu, S.Q.; Waaijer, S.J.; Zwager, M.C.; de Vries, E.G.; van der Vegt, B.; Schröder, C.P. Tumor-associated macrophages in breast cancer: Innocent bystander or important player? Cancer Treat. Rev. 2018, 70, 178–189. [Google Scholar] [CrossRef] [Green Version]
  101. Keegan, T.H.; DeRouen, M.C.; Press, D.J.; Kurian, A.W.; Clarke, C.A. Occurrence of breast cancer subtypes in adolescent and young adult women. Breast Cancer Res. 2012, 14, R55. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  102. Werner, H.M.J.; Mills, G.B.; Ram, P.T. Cancer Systems Biology: A peek into the future of patient care? Nat. Rev. Clin. Oncol. 2014, 11, 167–176. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  103. Bekisz, S.; Geris, L. Cancer modeling: From mechanistic to data-driven approaches, and from fundamental insights to clinical applications. J. Comput. Sci. 2020, 46, 101198. [Google Scholar] [CrossRef]
  104. Heymann, M.F.; Heymann, D. Immune environment and osteosarcoma. In Osteosarcoma-Biology, Behavior and Mechanisms; InTech: London, UK, 2017; pp. 105–120. [Google Scholar]
  105. Griguolo, G.; Pascual, T.; Dieci, M.V.; Guarneri, V.; Prat, A. Interaction of host immunity with HER2-targeted treatment and tumor heterogeneity in HER2-positive breast cancer. J. ImmunoTherapy Cancer 2019, 7, 1–14. [Google Scholar] [CrossRef]
  106. Le, T.; Su, S.; Kirshtein, A.; Shahriyari, L. Data-Driven Mathematical Model of Osteosarcoma. Cancers 2021, 13, 2367. [Google Scholar] [CrossRef]
  107. Asano, Y.; Kashiwagi, S.; Goto, W.; Kurata, K.; Noda, S.; Takashima, T.; Onoda, N.; Tanaka, S.; Ohsawa, M.; Hirakawa, K. Tumour-infiltrating CD8 to FOXP3 lymphocyte ratio in predicting treatment responses to neoadjuvant chemotherapy of aggressive breast cancer. Br. J. Surg. 2016, 103, 845–854. [Google Scholar] [CrossRef] [Green Version]
  108. Yasuda, K.; Nirei, T.; Sunami, E.; Nagawa, H.; Kitayama, J. Density of CD4 (+) and CD8 (+) T lymphocytes in biopsy samples can be a predictor of pathological response to chemoradiotherapy (CRT) for rectal cancer. Radiat. Oncol. 2011, 6, 1–6. [Google Scholar] [CrossRef] [Green Version]
  109. Riemann, D.; Hase, S.; Fischer, K.; Seliger, B. Granulocyte-to-dendritic cell-ratio as marker for the immune monitoring in patients with renal cell carcinoma. Clin. Transl. Med. 2014, 3, 1–6. [Google Scholar] [CrossRef] [PubMed]
  110. García-Tuñón, I.; Ricote, M.; Ruiz, A.; Fraile, B.; Paniagua, R.; Royuela, M. Influence of IFN-gamma and its receptors in human breast cancer. BMC Cancer 2007, 7, 1–11. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  111. Gooch, J.L.; Herrera, R.E.; Yee, D. The role of p21 in interferon gamma-mediated growth inhibition of human breast cancer cells. Cell Growth Differ. Mol. Biol. J. Am. Assoc. Cancer Res. 2000, 11, 335–342. [Google Scholar]
  112. Ning, Y.; Riggins, R.; Mulla, J.; Chung, H.; Zwart, A.; Clarke, R. Interferon gamma restores breast cancer sensitivity to Fulvestrant by regulating STAT1, IRF1, NFkB, BCL2 family members and signaling to a caspase-dependent apoptosis. Mol Cancer Ther 2010, 9, 1274–1285. [Google Scholar] [CrossRef] [Green Version]
  113. Kopreski, M.S.; Lipton, A.; Harvey, H.A.; Kumar, R. Growth inhibition of breast cancer cell lines by combinations of anti-P185HER2 monoclonal antibody and cytokines. Anticancer Res. 1996, 16, 433–436. [Google Scholar]
  114. Bhat, P.; Leggatt, G.; Waterhouse, N.; Frazer, I.H. Interferon-γ derived from cytotoxic lymphocytes directly enhances their motility and cytotoxicity. Cell Death Dis. 2017, 8, e2836. [Google Scholar] [CrossRef] [Green Version]
  115. Zhao, Q.; Tong, L.; He, N.; Feng, G.; Leng, L.; Sun, W.; Xu, Y.; Wang, Y.; Xiang, R.; Li, Z. IFN-γ mediates graft-versus-breast cancer effects via enhancing cytotoxic T lymphocyte activity. Exp. Ther. Med. 2014, 8, 347–354. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  116. Gao, Y.; Chen, X.; He, Q.; Gimple, R.C.; Liao, Y.; Wang, L.; Wu, R.; Xie, Q.; Rich, J.N.; Shen, K.; et al. Adipocytes promote breast tumorigenesis through TAZ-dependent secretion of Resistin. Proc. Natl. Acad. Sci. USA 2020, 117, 33295–33304. [Google Scholar] [CrossRef]
  117. Voss, A.; Voss, J. A fast numerical algorithm for the estimation of diffusion model parameters. J. Math. Psychol. 2008, 52, 1–9. [Google Scholar] [CrossRef]
  118. Parra-Rojas, C.; Hernandez-Vargas, E.A. PDEparams: Parameter fitting toolbox for partial differential equations in python. Bioinformatics 2020, 36, 2618–2619. [Google Scholar] [CrossRef]
  119. Vyshemirsky, V.; Girolami, M. BioBayes: A software package for Bayesian inference in systems biology. Bioinformatics 2008, 24, 1933–1934. [Google Scholar] [CrossRef] [Green Version]
  120. Xun, X.; Cao, J.; Mallick, B.; Maity, A.; Carroll, R.J. Parameter estimation of partial differential equation models. J. Am. Stat. Assoc. 2013, 108, 1009–1020. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  121. Zandarashvili, L.; Sahu, D.; Lee, K.; Lee, Y.S.; Singh, P.; Rajarathnam, K.; Iwahara, J. Real-time Kinetics of High-mobility Group Box 1 (HMGB1) Oxidation in Extracellular Fluids Studied by in Situ Protein NMR Spectroscopy. J. Biol. Chem. 2013, 288, 11621–11627. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  122. Valley, C.C.; Solodin, N.M.; Powers, G.L.; Ellison, S.J.; Alarid, E.T. Temporal variation in estrogen receptor-? protein turnover in the presence of estrogen. J. Mol. Endocrinol. 2008, 40, 2334. [Google Scholar] [CrossRef] [Green Version]
  123. Lee, S.M.; Suen, Y.; Chang, L.; Bruner, V.; Qian, J.; Indes, J.; Knoppel, E.; van de Ven, C.; Cairo, M.S. Decreased interleukin-12 (IL-12) from activated cord versus adult peripheral blood mononuclear cells and upregulation of interferon- gamma, natural killer, and lymphokine-activated killer activity by IL- 12 in cord blood mononuclear cells. Blood 1996, 88, 945–954. [Google Scholar] [CrossRef]
  124. Scherer, P.E. The many secret lives of adipocytes: Implications for diabetes. Diabetologia 2019, 62, 223–232. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  125. Strawford, A.; Antelo, F.; Christiansen, M.; Hellerstein, M. Adipose tissue triglyceride turnover, de novo lipogenesis, and cell proliferation in humans measured with 2H2O. Am. J. Physiol.-Endocrinol. Metab. 2004, 286, E577–E588. [Google Scholar] [CrossRef] [Green Version]
  126. Ryu, E.B.; Chang, J.M.; Seo, M.; Kim, S.A.; Lim, J.H.; Moon, W.K. Tumour volume doubling time of molecular breast cancer subtypes assessed by serial breast ultrasound. Eur. Radiol. 2014, 24, 2227–2235. [Google Scholar] [CrossRef]
Figure 1. Interaction network. The main interaction network of cells and molecules in breast tumors, modeled in this paper. Variables of the model with their descriptions are given in Table 1.
Figure 1. Interaction network. The main interaction network of cells and molecules in breast tumors, modeled in this paper. Variables of the model with their descriptions are given in Table 1.
Jpm 11 01031 g001
Figure 2. Immune cell frequencies in each cluster. Clusters were obtained by applying K-means clustering to the percentages of 22 immune cell types in breast tumors. The most variant cells among clusters are shown in this figure.
Figure 2. Immune cell frequencies in each cluster. Clusters were obtained by applying K-means clustering to the percentages of 22 immune cell types in breast tumors. The most variant cells among clusters are shown in this figure.
Jpm 11 01031 g002
Figure 3. More Clinical features of each cluster. Subfigures (A,B), respectively, show box plots of patients’ age at diagnosis and survival months in each cluster. Subfigure (C) demonstrates Kaplan–Meier curves of overall survival probability across five clusters. Asterisks in the figures show the significance levels with Mann–Whitney-Wilcoxon (MWW) statistical test where, ns: no significance, ***: 0.0001 < p 0.001 , ****: p 0.0001 .
Figure 3. More Clinical features of each cluster. Subfigures (A,B), respectively, show box plots of patients’ age at diagnosis and survival months in each cluster. Subfigure (C) demonstrates Kaplan–Meier curves of overall survival probability across five clusters. Asterisks in the figures show the significance levels with Mann–Whitney-Wilcoxon (MWW) statistical test where, ns: no significance, ***: 0.0001 < p 0.001 , ****: p 0.0001 .
Jpm 11 01031 g003
Figure 4. Clinical features of each cluster. Subfigures (AD) show the percentage of patients with different subtypes of breast cancer, survival status, HER2 status, and ER status, respectively.
Figure 4. Clinical features of each cluster. Subfigures (AD) show the percentage of patients with different subtypes of breast cancer, survival status, HER2 status, and ER status, respectively.
Jpm 11 01031 g004
Figure 5. Dynamics of all variables. Dynamics of variables of the model over 3000 days. The different color lines describe the dynamics of different clusters.
Figure 5. Dynamics of all variables. Dynamics of variables of the model over 3000 days. The different color lines describe the dynamics of different clusters.
Jpm 11 01031 g005
Figure 6. Sensitivity analysis. Sensitivity level of the most sensitive parameters for cancer and total cell population, respectively.
Figure 6. Sensitivity analysis. Sensitivity level of the most sensitive parameters for cancer and total cell population, respectively.
Jpm 11 01031 g006
Figure 7. Sensitivity analysis. Other sensitive parameters for cancer and total cell population.
Figure 7. Sensitivity analysis. Other sensitive parameters for cancer and total cell population.
Jpm 11 01031 g007
Figure 8. Dynamics of cancer after the assumptions of the sensitive parameters were modified. Subfigures (AF) present dynamics after scaling the assumptions (24)–(29), respectively. The transparent region was the result of 10% perturbation of all the sensitive parameters from Section 3.3.
Figure 8. Dynamics of cancer after the assumptions of the sensitive parameters were modified. Subfigures (AF) present dynamics after scaling the assumptions (24)–(29), respectively. The transparent region was the result of 10% perturbation of all the sensitive parameters from Section 3.3.
Jpm 11 01031 g008
Table 1. Patient data correspondence with variables. Correspondence among the model variables and the gene expression data of the primary tumors and deconvolution results.
Table 1. Patient data correspondence with variables. Correspondence among the model variables and the gene expression data of the primary tumors and deconvolution results.
VariableNameData Used
T N Naive T-cellsCombination of CD4 naive and memory resting T-cells and resting NK cells
T h Helper T-cellsCombination of memory activated CD4 T-cells and follicular helper T-cells
T C Cytotoxic cellsCombination of CD8 T-cells and activated NK cells
T r Regulatory T-cellsRegulatory T-cells
D N Naive dendritic cellsNaive dendritic cells
DActivated dendritic cellsActivated dendritic cells
M N Naive MacrophagesCombination of Macrophages M0 and Monocytes
MMacrophagesCombination of M1 and M2 Macrophages
CCancer cellsEstimated
NNecrotic cellsEstimated
ACancer Associated AdipocytesAssumed to be twice the total number of immune cells
HHMGB1HMGB1 gene expression
I L 12 IL-12IL12A and IL12B gene expressions
I L 10 IL-10IL10 gene expression
EEstrogenESR1 and ESR2 gene expressions
I γ IFN- γ IFNG gene expressions
I L 6 IL-6IL6 gene expression
Table 2. Initial conditions of the model variables. Values of the initial conditions were obtained from the averages of the smallest tumors in METABRIC and TCGA data.
Table 2. Initial conditions of the model variables. Values of the initial conditions were obtained from the averages of the smallest tumors in METABRIC and TCGA data.
Cluster T N T h T C T r D N D M N M C
1 2.62 × 10 2 1.93 × 10 3 2.94 × 10 3 5.55 × 10 2 9.90 × 10 3 9.90 × 10 3 1.34 × 10 4 5.37 × 10 3 1.56 × 10 3
2 3.78 × 10 3 2.16 × 10 3 2.95 × 10 3 5.04 × 10 2 1.07 × 10 3 9.83 × 10 3 2.63 × 10 3 1.12 × 10 4 2.31 × 10 3
3 2.93 × 10 3 1.15 × 10 3 1.39 × 10 3 3.68 × 10 1 3.22 3.62 × 10 2 6.75 × 10 3 9.67 × 10 3 3.53 × 10 3
4 4.78 × 10 3 4.60 × 10 3 2.76 × 10 3 1.66 × 10 3 3.51 × 10 2 1.34 × 10 3 2.29 × 10 3 4.69 × 10 3 7.96 × 10 3
5 3.78 × 10 3 1.33 × 10 3 3.11 × 10 3 1.05 × 10 3 5.16 × 10 2 1.03 × 10 2 2.85 × 10 3 6.07 × 10 3 5.57 × 10 3
N A H IL 12 IL 10 E I γ IL 6
1 1.21 × 10 2 4.89 × 10 4 5.06 5.82 2.90 7.62 3.05 3.67
2 3.16 × 10 2 4.86 × 10 4 5.01 6.40 2.93 8.76 2.93 3.02
3 3.11 × 10 2 4.39 × 10 4 5.29 5.25 2.79 9.16 2.81 4.27
4 2.71 × 10 3 4.49 × 10 4 5.29 6.88 3.27 6.38 4.10 3.40
5 3.73 × 10 2 3.74 × 10 4 6.16 5.67 2.66 8.63 2.70 3.03
Table 3. Steady state values of the model variables. Large tumors in each cluster were grouped, and their average values were found for each variable, to be used in the parameter estimation.
Table 3. Steady state values of the model variables. Large tumors in each cluster were grouped, and their average values were found for each variable, to be used in the parameter estimation.
Cluster T N T h T C T r D N D M N M C
1 4.55 × 10 3 3.87 × 10 3 2.44 × 10 3 1.92 × 10 3 1.26 × 10 2 3.28 × 10 2 1.80 × 10 4 1.39 × 10 4 9.03 × 10 4
2 1.13 × 10 4 4.77 × 10 3 4.02 × 10 3 1.55 × 10 3 3.27 × 10 2 6.10 × 10 2 6.96 × 10 3 1.62 × 10 4 1.12 × 10 5
3 5.73 × 10 3 3.70 × 10 3 2.79 × 10 3 1.51 × 10 3 4.00 × 10 2 3.52 × 10 2 8.77 × 10 3 2.07 × 10 4 9.68 × 10 4
4 4.14 × 10 3 5.95 × 10 3 7.71 × 10 3 1.66 × 10 3 4.99 × 10 2 4.37 × 10 2 9.81 × 10 3 1.59 × 10 4 9.54 × 10 4
5 5.69 × 10 3 3.91 × 10 3 5.56 × 10 3 6.16 × 10 2 7.04 × 10 2 2.22 × 10 2 6.45 × 10 3 1.52 × 10 4 9.90 × 10 4
N A H IL 12 IL 10 E I γ IL 6
1 1.15 × 10 4 9.01 × 10 4 5.61 6.42 3.17 8.86 3.21 3.42
2 1.02 × 10 4 9.14 × 10 4 3.84 2.84 1.36 4.03 1.19 1.28
3 8.11 × 10 3 8.80 × 10 4 4.49 4.02 2.13 7.68 1.97 2.14
4 1.54 × 10 4 9.22 × 10 4 7.54 1.04 5.08 1.37 5.72 5.80
5 1.53 × 10 4 7.67 × 10 4 7.70 1.02 5.05 1.50 5.17 6.01
Table 4. Three important parameter values. Values of the parameters δ C T c , δ C I γ , and δ C without scaling and after scaling the assumption (24) by factors 0.2 and 5. These are the most effected parameters and correspond to the dynamics in Figure 8A.
Table 4. Three important parameter values. Values of the parameters δ C T c , δ C I γ , and δ C without scaling and after scaling the assumption (24) by factors 0.2 and 5. These are the most effected parameters and correspond to the dynamics in Figure 8A.
Clusters Without Scaling Scale = 0.2 Scale = 5
Cluster 1 δ C T c = 0.00440
δ C I γ = 0.00274
δ C = 0.04492
δ C T c = 0.00895
δ C I γ = 0.00111
δ C = 0.01827
δ C T c = 0.00124
δ C I γ = 0.00390
δ C = 0.06341
Cluster 2 δ C T c = 0.00528
δ C I γ = 0.00074
δ C = 0.03273
δ C T c = 0.00953
δ C I γ = 0.00027
δ C = 0.01182
δ C T c = 0.00163
δ C I γ = 0.00115
δ C = 0.05063
Cluster 3 δ C T c = 0.00494
δ C I γ = 0.00165
δ C = 0.04421
δ C T c = 0.01027
δ C I γ = 0.00069
δ C = 0.01839
δ C T c = 0.00137
δ C I γ = 0.00230
δ C = 0.06150
Cluster 4 δ C T c = 0.00805
δ C I γ = 0.00283
δ C = 0.02606
δ C T c = 0.01183
δ C I γ = 0.00083
δ C = 0.00766
δ C T c = 0.00310
δ C I γ = 0.00546
δ C = 0.05018
Cluster 5 δ C T c = 0.00676
δ C I γ = 0.00298
δ C = 0.03037
δ C T c = 0.01068
δ C I γ = 0.00094
δ C = 0.00960
δ C T c = 0.00238
δ C I γ = 0.00526
δ C = 0.05355
Table 5. Three important parameter values. Values of the parameters δ C T c , δ C I γ , and δ C without scaling and after scaling the assumption (25) by factors 0.2 and 5. These are the most effected parameters and correspond to the dynamics in Figure 8B.
Table 5. Three important parameter values. Values of the parameters δ C T c , δ C I γ , and δ C without scaling and after scaling the assumption (25) by factors 0.2 and 5. These are the most effected parameters and correspond to the dynamics in Figure 8B.
Cluster Without Scaling Scale = 0.2 Scale = 5
Cluster 1 δ C T c = 0.00440
δ C I γ = 0.00274
δ C = 0.04492
δ C T c = 0.00168
δ C I γ = 0.00104
δ C = 0.08556
δ C T c = 0.00652
δ C I γ = 0.00406
δ C = 0.013308
Cluster 2 δ C T c = 0.00528
δ C I γ = 0.00074
δ C = 0.03273
δ C T c = 0.00161
δ C I γ = 0.00023
δ C = 0.05002
δ C T c = 0.00967
δ C I γ = 0.00136
δ C = 0.01199
Cluster 3 δ C T c = 0.00494
δ C I γ = 0.00165
δ C = 0.04421
δ C T c = 0.00155
δ C I γ = 0.00052
δ C = 0.06927
δ C T c = 0.00879
δ C I γ = 0.00294
δ C = 0.01574
Cluster 4 δ C T c = 0.00805
δ C I γ = 0.00283
δ C = 0.02606
δ C T c = 0.00467
δ C I γ = 0.00164
δ C = 0.07563
δ C T c = 0.00941
δ C I γ = 0.00331
δ C = 0.00609
Cluster 5 δ C T c = 0.00676
δ C I γ = 0.00298
δ C = 0.03037
δ C T c = 0.00358
δ C I γ = 0.00158
δ C = 0.08036
δ C T c = 0.00823
δ C I γ = 0.00363
δ C = 0.00739
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Mohammad Mirzaei, N.; Su, S.; Sofia, D.; Hegarty, M.; Abdel-Rahman, M.H.; Asadpoure, A.; Cebulla, C.M.; Chang, Y.H.; Hao, W.; Jackson, P.R.; et al. A Mathematical Model of Breast Tumor Progression Based on Immune Infiltration. J. Pers. Med. 2021, 11, 1031. https://doi.org/10.3390/jpm11101031

AMA Style

Mohammad Mirzaei N, Su S, Sofia D, Hegarty M, Abdel-Rahman MH, Asadpoure A, Cebulla CM, Chang YH, Hao W, Jackson PR, et al. A Mathematical Model of Breast Tumor Progression Based on Immune Infiltration. Journal of Personalized Medicine. 2021; 11(10):1031. https://doi.org/10.3390/jpm11101031

Chicago/Turabian Style

Mohammad Mirzaei, Navid, Sumeyye Su, Dilruba Sofia, Maura Hegarty, Mohamed H. Abdel-Rahman, Alireza Asadpoure, Colleen M. Cebulla, Young Hwan Chang, Wenrui Hao, Pamela R. Jackson, and et al. 2021. "A Mathematical Model of Breast Tumor Progression Based on Immune Infiltration" Journal of Personalized Medicine 11, no. 10: 1031. https://doi.org/10.3390/jpm11101031

APA Style

Mohammad Mirzaei, N., Su, S., Sofia, D., Hegarty, M., Abdel-Rahman, M. H., Asadpoure, A., Cebulla, C. M., Chang, Y. H., Hao, W., Jackson, P. R., Lee, A. V., Stover, D. G., Tatarova, Z., Zervantonakis, I. K., & Shahriyari, L. (2021). A Mathematical Model of Breast Tumor Progression Based on Immune Infiltration. Journal of Personalized Medicine, 11(10), 1031. https://doi.org/10.3390/jpm11101031

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop