Skip to main content
Advertisement
  • Loading metrics

Investigating key cell types and molecules dynamics in PyMT mice model of breast cancer through a mathematical model

Abstract

The most common kind of cancer among women is breast cancer. Understanding the tumor microenvironment and the interactions between individual cells and cytokines assists us in arriving at more effective treatments. Here, we develop a data-driven mathematical model to investigate the dynamics of key cell types and cytokines involved in breast cancer development. We use time-course gene expression profiles of a mouse model to estimate the relative abundance of cells and cytokines. We then employ a least-squares optimization method to evaluate the model’s parameters based on the mice data. The resulting dynamics of the cells and cytokines obtained from the optimal set of parameters exhibit a decent agreement between the data and predictions. We perform a sensitivity analysis to identify the crucial parameters of the model and then perform a local bifurcation on them. The results reveal a strong connection between adipocytes, IL6, and the cancer population, suggesting them as potential targets for therapies.

Author summary

One of the outstanding challenges of the mathematical modeling of cancer progression is the existence of many unknown parameters. In this work, we develop a data-driven mathematical model of breast cancer progression by deriving a system of ordinary differential equations for the interaction networks of key cell types and molecules in breast tumors. To overcome the limitations of unknown parameters, we utilize a time course data of a PyMT mice model of breast cancer and estimate parameters using an optimization method. Although the predicted dynamics of cancer and necrotic cells using the obtained values of parameters are in good agreement with the data, the predicted values for a few other variables do not match the data. This might indicate that there are some other key interactions that have not been modeled, and/or there is a noise in the data. The sensitivity and bifurcation analyses show that the most important parameters in controlling the cancer cells population are the proliferation and death rates of cancer cells and adipocytes. These results are in agreement with some biological and clinical studies of breast cancer, which have reported a link between adipocytes and breast cancer progression.

Introduction

Breast cancer is known to be one of the most common types of cancer in women. In 2021, breast cancer accounted for 50% of all new diagnoses in the USA, with a projected death of 43,600 [1]. Breast cancer can be divided into different subtypes through molecular-level analysis of gene expression patterns. These subtypes are defined as luminal A (LumA), luminal B (LumB), luminal/human epidermal growth factor receptor 2 (HER2), HER2 enriched, basal-like, and triple-negative breast cancer (TNBC) nonbasal [2]. LumA is the most common type with the lowest mortality rate among other subtypes [3]. Most cancer treatments available today focus on killing cancer cells as well as removing them via surgery [46]. While these treatments may cure cancer, in some cases, cancer metastasizes to other areas of the body after treatments [7, 8]. Furthermore, it is estimated that 70-80% stage four metastatic breast cancer patients die within five years [9]. Understanding the biology of cancer as a whole intricate system of interactions is crucial for assessing the invasiveness of cancer and obtaining effective treatments.

The tumor microenvironment is a mixture of many cell types and molecules. The importance of the cells and molecules interaction networks within the microenvironment in tumor development has attracted many researchers from many disciplines; whether it is the cancer progression [1012] or its response to treatments [1315]. The interactions between immune cells, cancer cells, necrotic cells, and adipocytes result in interesting dynamics and lead to the secretion of important cytokines such as HMGB1, IL12, IL10, and IL6 [16, 17]. These cytokines are often subjects of targeted therapies [1821].

In-vivo investigation of the tumor microenvironment can be costly and straining for both patients and scientists. Therefore, mouse models have been one of the most popular methods of studying cancers’ initiation and progression. There are many different approaches such as transgenic [2225], gene targeting [2628], RNA interference [2931], and many more to create mouse models. The mammary specific polyomavirus middle T antigen overexpression mouse model (MMTV-PyMT) is the most popular mouse model used in cancer studies, especially in breast cancer which qualifies as a transgenic approach. These models show the signaling of receptor tyrosine kinases commonly activated in many human progressive tumors, including breast cancer [32, 33]. Therefore, they are very suitable as they closely mimic the development of breast cancer in humans.

Mathematical models have enabled scientists to gain insight into biological phenomenon, which are either obscure or costly to experiment. Modeling cardiovascular system [3436], disease spread [3739], muscle function [40, 41] and ocular disease [42, 43] are just a few examples of them. Cancer is among one of the most mathematically modeled diseases, some aim to answer questions about cancer as a whole system of biological and chemical interactions [4447], some investigate the mechanical properties of a cancerous tissue [48, 49], and some others focus on modeling cancer response to different treatments [5052]. The most desirable features of mathematical models of cancer, including stochastic [5355] and deterministic models [5658], are their ability to make good predictions, testing plausible biological hypotheses or generating clinically testable hypothesis. For example, a multiscale model of prostate cancer shows that low androgen levels may increase resistance to hormonal therapy and that treatment with 5α-reductase inhibitors may lead to more therapy-resistant cancer cells [59], and a data driven mathematical model predicts the response to FOLFIRI treatment for colon cancer patients [60]. Moreover, an agent-based model [61] of tumor progression indicate that while macrophages existence can increase the size of the tumor, an increase in their infiltration has a reverse effect. In another study, a hybrid agent-based model [62] of ductal carcinoma, which is a common type of breast cancer that starts in cells that line the milk ducts, finds that duct advance rates happen in two phases of an early exponential expansion, followed by a long-term steady linear expansion. Additionally, a free boundary mathematical model of the early detection of recurrences shows a relation between the size of the growing cancer and the total Serum uPAR mass in the cancer [63].

Simple models such as the logistic model and Gompertz model have helped in understanding the growth dynamics of the tumor, predicting the age of the tumor, etc. [64, 65]. These models are dependent on parameters such as proliferation and degradation rates which may depend on cells and molecules interactions in the tumor. Thus, such simple models have the capacity of expanding into more comprehensive models and include many factors through a system of ordinary differential equations (ODEs). This could give a better understanding of how the tumor develops to find optimal treatments [6670]. However, for mathematical models, obtaining parameter values is always a challenge due to the lack of biological data, especially time course dataset.

In this paper, We use a system of ODEs to describe the dynamics of some key players in the breast tumors microenvironment. We utilize a time course data set collected from three PyMT mice at four different stages of cancer’s progression [71] to estimate the parameters of our proposed mathematical model through an optimization approach. There are several methods for estimating model parameters, including Monte Carlo Hastings and steady-state assumption. Each of these methods has its benefits and drawbacks. Monte Carlo methods, in particular, necessitate a large number of simulations; therefore, they are extremely slow for large-scale problems. Although assuming steady-state evaluates the system parameters straightforwardly, some variables reach the steady-state prematurely and may not correspond with the data in some cases. Here, the least-squares optimization is used to handle parameter estimation in a feasible region for a breast cancer mouse model. Researchers have employed this method to evaluate the system parameters in various disciplines [7274]. We apply the least-square optimization method on the PyMT mouse time course data set to estimate the parameters of our model. We then assess the sensitivity of our model to its parameters using a sensitivity analysis based on a direct differential method. Finally, we locally investigate the effect of sensitive parameters using bifurcation plots. The results show interesting connections with important biological observations reported in the literature.

Materials and methods

Many cells and cytokines are involved in cancer development, but to avoid complexity, we only consider the most important ones give in Table 1. Fig 1 shows the interaction network between different cell types and molecules used in this paper.

thumbnail
Fig 1. Interaction network.

Diagram of interactions between different cell types and molecules in breast tumors for the mouse model, as it has been modeled in this paper. Variables of the model with their descriptions are given in Table 1.

https://doi.org/10.1371/journal.pcbi.1009953.g001

Cells and molecules interaction network—ODE model

T-cells.

We categorize T-cells into 4 major sub-types: naive (TN), helper (Th), cytotoxic (TC) and regulatory (Tr) T-cells; with naive T-cells being the only sub-type, which is not present in the tumor microenvironment and is mainly found in the lymphatic system [75]. Even though, introducing nonlinearity in ODEs for T-cells can prevent issues such as unlimited exponential growth, for simplicity, we just make activation rates of other T-cells proportional to the number of naive T-cells. In this way, we control the system with less complexity. Therefore, we describe the ODEs for helper, cytotoxic and regulatory T-cells prior to that of the naive T-cells.

CD4+ helper T-cells (Th).

Dendritic cells [76, 77], HMGB1 [78, 79], and IL-12 [66, 80] activate CD4+ helper T-cells, while regulatory T-cells [81, 82] and IL-10 [83, 84] inhibit them. Therefore, we use the following ODE to describe the dynamic of helper T-cells (1)

Cytotoxic T-Cells (Tc).

Dendritic cells [85, 86] and IL-12 activate naive CD8+ T-cells [66, 80]. On the other hand, regulatory T-cells [66, 82] and IL-10 [83, 84] suppress Cytotoxic T-Cells functionality. Due to the similarity between natural killer (NK) cells and Cytotoxic T-Cells in directly killing cancer cells, we assume this group includes both CD8+ T-cells and NK cells. Therefore, we model cytotoxic T-cells’ dynamics in the following way. (2)

Regulatory T-Cells (Tr).

Dendritic cells stimulate formation [87] and activation of regulatory T-cells [86, 88]. Hence, we have the following equation for the dynamics of T-reg cells. (3)

Naive T-Cells (TN).

Combining Eqs (1)–(3) for the activation of naive T-cells and adding an independent naive T-cells production rate , we get the following ODE for for naive T-cells (4)

Dendritic cells (D)

Cancer cells [66, 89] and HMGB1 [67, 9092] can activate dendritic cells. Moreover, cancer cells may promote natural death of dendritic cells in different ways [86]. Adding as the production rate of naive dendritic cells, we get the following system of equations for naive dendritic cells (DN) and activated dendritic cells (D) (5) (6)

Macrophages (M)

Macrophages have many phenotypes and can change them frequently. For simplicity, we avoid the break down of them into M1, M2, and other subsets, and we model all activated macrophages as a single variable denoted by M. Tumor associated macrophages (TAMs) are activated by IL-10 [9395]. Moreover, IL-12 activates M1 Macrophages [68, 93, 9698], while M2 macrophages are activated by helper T-cells secreted cytokines (IL-13 and IL-4) [93].

Denoting naive macrophages by MN, activated macrophages by M, and the production rate of naive macrophages by AM, we can write the following system of equations for the dynamics of naive and activated macrophages. (7) (8)

Cancer cells (C)

Traditionally, it is assumed that the growth rate of cancer cells is related to both the existing population and the available resources or space. Therefore, we model the proliferation of cancer cells by the logistic term [C](1 − [C]/C0), where C0 is the maximum capacity. In addition, IL-6 promotes the proliferation of cancer cells [99101]. Also, adipocytes, releasing metabolic substrates, promote proliferation of breast cancer cells [102, 103]. On the other hand, activated CD8+ T-cells kill cancer cells [66, 104, 105]. The dynamics of cancer cells is modeled by the following equation. (9)

Cancer associated adipocytes (A)

Adipocytes participate in a highly complex cycle orchestrated by cancer cells to promote tumor progression [106]. Therefore, after including them in Eq (9), for simplicity we consider an independent logistic model describing their dynamics. (10)

Necrotic cells (N)

Cells in the tumor microenvironment can die and turn into necrotic cells due to depletion of resources. This process is called necrosis and is known as a feature of tumors possessing an aggressive phenotype [107]. Necrotic cell death can replace other types of death for some cell types [45, 91]. As mentioned, activated CD8+ T-cells kill cancer cells [66, 104, 105] and we assume that death of cancer cells is the primary source of necrosis. Since a fraction of cancer cells can go through first becoming necrotic cells, the production rate of necrotic cells is modeled by the fraction (αNC) of dying cancer cells. (11)

Molecules

In this section, we describe ODEs that govern dynamics of molecules in our model.

HMGB1 (H).

High-mobility group box 1 (HMGB1) is known to be a prototypical damage-associated molecular pattern (DAMP) protein, which alarms the body about disturbances in homeostasis [108]. HMGB1 exert immune promoting activity by inducing angiogenesis, proliferation, and invasiveness of cancer cells [90]. HMGB1 is mainly produced by dendritic cells [90, 91, 109, 110], necrotic cells [66, 109, 111, 112], macrophages, [109, 113115], natural killer (NK) cells which behave like cytotoxic T-cells [109, 116118], and cancer cells [66, 90, 91].

Therefore, the dynamics of HMGB1 is modeled by the following equation. (12)

IL-12 (IL12).

IL-12, stimulates the differentiation of naive T-cells into helper T-cells. Macrophages and dendritic cells secrete IL-12 [66, 86, 97, 119]. Also, IL-12 is produced by Helper and cytotoxic T-cells [120]. We model the dynamics of IL-12 using the following equation. (13)

IL-10 (IL10).

IL-10 is secreted by macrophages [93, 121, 122], dendritic cells [86, 123125], T-reg cells [83, 126], CD4+ helper T-Cells [120, 127, 128], CD8+ cytotoxic T-cells [120, 126, 128], and cancer cells [87, 129]. Therefore, the dynamics of IL-10 is modeled in the following way. (14)

IL-6 (IL6).

The key cytokine that promotes the growth of cancer cells is IL-6 and is produced by cancer associated adipocytes [97, 99, 100, 130], macrophages [66, 93, 97, 98, 131, 132], and dendritic cells [66, 86, 120]. (15)

Mouse data analysis

For this study, we use the PyMT mice RNA-sequencing data available in the Gene Expression Omnibus (GEO) database as GSE76772 [71]. The PyMT gene expressions were acquired from 3 PyMT mice at four tumor progression stages: hyperplasia at week 6, adenoma/MIN at week 8, early carcinoma at week 10, and late carcinoma at week 12. The original study was designed to recognize gene expression similarities at different cancer stages. They used a directional RNA sequencing method to acquire the raw gene expression data. Later, they used statistical methods to remove transcriptionally inactive genes and get high confident normalized gene counts. We apply CIBERSORTx B-mode with the LM22 signature matrix [133] on the mentioned time-course gene expression data to estimate the relative abundance of each immune cell type in the tumor. Finally, we use expression values of genes encoding cytokines in the model and combined some immune cells to estimate the values of the model’s variables. Fig 2 shows the most variant immune cell frequencies for three mice at different time points. Since the deconvolution method only provides us the percentage of each immune cell type in primary tumors, we use the tumor size for each mouse to estimate the number of immune cells, cancer cells, and necrotic cells in each sample. Also, based on our observations, the average ratio of cancer cells to immune cells and necrotic cells is approximately 0.955:0.04:0.005 in mouse model breast tumors. Also, the epithelial cells density has been reported as 45 cells/mm3 in breast cancer [134]. Thus by choosing the scaling factor α = 45, the average density of cancer cells across all samples at each time is close to that value. So, we first calculate the total number of cells (TNC) for each mouse at a time point using where ti ∈ {6 weeks, 8 weeks, 10 weeks, 12 weeks} for i = 1, ⋯, 4. Using the TNC, we calculate the total number of cancer cells (TNCC), the total number of immune cells (TNIC), and the total number of necrotic cells (TNNC), using the ratio 0.955:0.04:0.005 and the following formulas See Table 2 for the values. The fraction is just the simplified ratio of cancer cells to necrotic cells 0.955 : 0.005.

thumbnail
Fig 2. Immune cell frequencies for each mouse at different time points.

Results acquired from deconvolution of gene expression data at each time point using CIBERSORTx B-mode.

https://doi.org/10.1371/journal.pcbi.1009953.g002

thumbnail
Table 2. Cells and molecules values for each mouse at different time points.

https://doi.org/10.1371/journal.pcbi.1009953.t002

Parameter estimation

For better stability, we perform the optimization process on the non-dimensionalized system, see the Non-dimensionalization appendix. In the following, for easier identification, we present matrix and vector quantities using boldface upper and lower case symbols, respectively.

Least-squares optimization.

The set of parameters, i.e., the coefficients, proliferation, and death rates, when no specified bounds are desired, can be evaluated by re-arranging the 15 ODEs expressed in Eqs (1)–(15), and solving a linear least-squares problem via the following formula (16) where A is the augmented matrix of mice’ data obtained by re-arranging, whose element ij is the ith observation of the jth variable, i.e., the value of each variable at specific time reported in Table 2. The rates are evaluated using a central finite difference method and are collected in vector b. The solution vector, θ, provides the approximation for the 57 model parameters. The values of the carrying capacities C0 and A0 are determined based on the data, thus are considered known quantities, which keeps the system linear. In this study the parameters of the system are all non-negative; however, the expression in Eq (16) does not enforce any bounds on the solution that may result in negative values for the parameters. To enforce non-negativity, the following linear least-squares optimization problem is solved (17) that finds the solutions in the feasible region. The parameters’ bounds are the only constraint imposed. In this study, we solve the optimization problem above by setting θmin = 10−5. In Eq (17), the solution is found for the minimum residual value in an iterative process where no inversion, such as the one in Eq (16), is needed. This approach of finding a solution using optimization (i.e., the iterative process) has been employed successfully in prior studies [135, 136]. See the Parameter values appendix where the optimized set of non-negative parameters are reported in Table 3.

Sensitivity analysis

Sensitivity analysis is generally used to assess the sensitivity of the model’s output to system parameters [137]. To identify the most crucial parameters affecting the dynamics of cancer and the total number of cells, we perform sensitivity analysis on these quantities.

We use to show non-dimensional variables. For a generic ODE system of the form (18) where and θ = 〈θ1, ⋯, θ57〉 are vectors of state variables and parameters of our model, respectively, the first order local sensitivity of a variable with respect to a parameter θi is evaluated by (19) To obtain the sensitivity vector, s = 〈s1, ⋯, sM〉, we use a direct differential method. That is, we differentiate Eq (18) with respect to θi to get (20) and then we use a forward Euler discretization in time for Eq (20) to find . The sensitivity of each parameter in the neighborhood of a chosen parameter set Ω(θ) is then defined as (21) This neighborhood is created by perturbing our original parameter set by 10% and the integration is carried out numerically with sparse grid points [138, 139].

In this study, because some of our state variables do not reach steady-state within an experimentally reasonable time interval, we use a direct differential method rather than a steady-state method. As mentioned before, the latest data point was extracted at 12 weeks. However, our observations show that we need to continue the simulation for much longer than experimentally feasible so that all variables reach the steady-state. For this reason, we use a direct differential approach up to 18 weeks to obtain the sensitivities.

Results

Dynamics

We use the optimized parameters from the least-squares optimization and substitute them in the system of ODEs to obtain the dynamics of the system variables. As mentioned before, mice data was collected at weeks 6, 8, 10, and 12 corresponding to 42, 56, 70, and 84 days. However, we shifted our time interval so that 6 weeks becomes the time origin (t = 0) and 12 weeks maps to 42 days. We also continue our ODE solutions further than 42 days (up to 126 days) to match our sensitivity analysis. Fig 3 shows a comparison between the solution of ODEs and mouse data for the time period of 126 days. For a clear comparison of the results, the solution of ODEs and mice data are shown in one plot.

thumbnail
Fig 3. Comparison of the dynamics and mouse data.

By solving the system of ODEs, the mouse data are compared to the values estimated using the optimal parameters obtained from the least-squares optimization. Time zero corresponds to 6 weeks which is the beginning of the mice data sampling.

https://doi.org/10.1371/journal.pcbi.1009953.g003

Despite the fact that the dynamics of cancer and necrotic cells fit the data well, a few predictions are less in accordance with the data points due to considerable disparities in the available data for some state variables. This can be remedied by more time course data decreasing the noise. However, it is considered a limitation of our study at this point.

Based on ODE simulations, naive T-cells for mice 1 and 3 show a quick decrease and then increase, while mouse 2 is strictly increasing. Helper T-cells for mice 1 and 3 quickly decrease to a very small steady-state, while for mouse 2, it gently increases. Cytotoxic cells behave the same way as the helper T-cells, and regulatory T-cells are generally decreasing with a negligible change in Mouse 2 compared to the other mice. Generally, for T-cells, our simulations are not in a good agreement with the data. Mouse 3 shows the best agreement in naive and cytotoxic T-cells, mouse 2 in regulatory T-cells and mouse 1 in helper T-cells.

The ODE results show that naive dendritic cells increase and then decrease in all 3 mice. Activated dendritic cells sharply decrease to a small steady-state in mice 1 and 3. In mouse 2, we can see a fluctuation in dendritic cells, but these changes are very small and negligible and eventually settles at a small value like the other two mice. Compared to T-cells, dendritic cells show a much better agreement with data in all three mice. The number of naive macrophages decrease in all three mice. Activated macrophages also decrease with similar rates in all three mice. Again, the simulation results are not in good agreement with the data. Mouse 1 shows the best agreement in naive macrophages and, mouse 3 in activated macrophages.

The predicted dynamics of cancer cells in all three mice follow the data very closely; dynamics in all three mice are similar, and they reach a steady-state value within 75 days. Given the closeness of reported data points in three mice, this similarity in their behavior is expected. Adipocytes’ dynamics also reach the steady-state values in all three mice. They do so by slowly decreasing in mouse 1 and slowly increasing in mouse 2. Mouse 3 undergoes a sharper increase before it converges to its steady-state value. Interestingly, in all three mice, the number of adipocytes converges to the same value. In mouse 3, we see a fair agreement between the dynamic of adipocytes and the data unlike the other two mice. Necrotic cells, like cancer cells, show a logistic growth and a perfect match with the data. We see almost the same behavior in all 3 mice reaching the same steady-state value. This is expected since the source of necrosis in the model is the death of cancer cells.

HMGB1 dynamics show a good agreement with the data, in all three cases. They start with a sharp increase from the initial value followed by a fast decrease. IL12 dynamics in mice 1 and 3 quickly decrease to a steady-state, but mouse 2 shows a mild, strictly increasing behavior. Mouse 1 and 3 also follow the data closely unlike mouse 2. IL10 decreases in all cases with a rather steeper slope in mouse 2. As a matter of fact, mouse 2 is the only mouse for which IL10 reaches its steady-state value within 125 days. Finally, IL6 decreases to a steady-state within the simulation time for all three mice, with mouse 3 showing the steepest descend followed by mouse 1. In general, we see a fair correspondence between dynamics and data in all 3 mice for both IL10 and IL 6. However, this correspondence is better in mice 1 and 3 than mouse 2.

Now we comment on the interactions, which we believe are responsible for the differences in the dynamics. For activated T-cell subtypes, we can see from equations Eqs (1) and (2) that dendritic cells, IL12 and HMGB1 are involved in their promotion, and regulatory T-cells and IL10 are involved in their inhibition process.

Helper T-cells are produced by activation of naive T-cells and inhibited by T-reg cell. For helper T-cells, we can see that mouse 2 has significantly fewer T-reg cells in the long run, and its IL12 levels are increasing, unlike the other two mice. These two behaviors contribute to a rise in the number of helper T-cells in this mouse. The same goes for cytotoxic cells. T-reg cell levels are only dependent on dendritic cells. Dendritic cells in Mice 1 and 3 start at high values contributing to some activation of T-reg cells and then a quick depletion. This explains the decreasing of T-reg cells in these mice. As for mouse 2, the number of dendritic cells starts low and stays low, resulting in the low numbers of T-reg cells and the sharp decrease in this case. Finally, Eq (5) shows that in the model, naive T-cells are produced via a constant rate, while activation of other T-cell subtypes contributes to their depletion. All three mice show a growing trend for naive T-cells. Given that the other subtypes are generally low or depleting (mice 1 and 3) or have a very gentle growth (mouse 2), the natural production of naive T-cells dominates the process.

Dendritic cells are activated and inhibited by cancer cells, see Eq (6). Therefore, its activation by HMGB1 can be the key to the observed behaviors. The sudden decrease in HMGB1 in all mice can be the reason that the decaying effects in Eq (6) have taken over so quickly. For the naive dendritic cells, the sudden surge results from the sudden drop in activated dendritic cells. However, natural decay takes over later.

By looking at Eq (7), we can see that macrophages get activated by IL10, IL12, and helper T-cells, and the only cause of death considered for them in this model is their natural decay. All the activators decrease in mice 1 and 3. In mouse 2, we have a slight increase in helper T-cells and IL12, and maybe that is why mouse 2 has the highest level of activated macrophages. However, at the end, these effects are not enough to win over the natural decay, and hence we see an overall decrease in all three cases. For naive macrophages, we can see that the death rate in Table 3, is an order of magnitude larger than the natural production rate. As a result, there is a decrease in naive macrophages in all three mice.

In the model, either cancer cells production happens independently or is promoted by adipocytes and IL6. They can also die naturally or by cytotoxic cells. As mentioned before, adipocytes in all three mice converge to the same value. Also, IL6 behavior is similar across all cases, while cytotoxic cells follow a different pattern in mouse 2. In fact, the increase in Tc observed in mouse 2 agrees with CD8 frequency shown in Fig 2. Adipocytes and IL6 play crucial roles in the cancer dynamics given their roles and the similarity between cancer dynamics and their dynamics. Adipocytes are modeled independently. They follow a logistic population model, and depending on their estimated growth and decay rate; they increase or decrease and then saturate. Also, since necrotic cells are produced as a result of cancer cells’ death, their dynamics are self-explanatory.

HMGB1 is produced proportional to the number of dendritic cells, necrotic cells, macrophages, cytotoxic cells, and cancer cells and decays naturally. Among the mentioned cell types, cancer and necrotic cells are the only ones whose numbers increase in time, and the rest of the cell types decrease. But, the parameter estimation shows that cancer cells and necrotic cells have the smallest production rates (two orders of magnitude smaller than the natural decay rate). Therefore, HMGB1 won’t be affected by them and decreases in time.

IL12 is produced proportional to the number of dendritic cells, macrophages, cytotoxic and helper T-cells and decays naturally. Higher levels of macrophages and increasing levels of helper and cytotoxic T-cells in mouse 2 is the reason for its mild increase, unlike mice 1 and 3.

IL10 is produced proportional to the number of dendritic cells, macrophages, cytotoxic, helper, and regulatory T-cells plus cancer cells and decays naturally. Given its large production rate by helper T-cells and its even larger decay rate, this cytokine is more significantly affected by these two parameters. Therefore, it decays quickly in all mice. However, we can see its steady-state value in mice 2 within the simulation time interval. This is mostly due to the increasing helper T-cells, which dampens the decreasing effect of its natural decay.

Finally, IL6 is produced by adipocytes, macrophages, and dendritic cells and is removed naturally. The dynamics of IL-6 is heavily depend on its production rate by dendritic cells, because its production by dendritic cells and its natural decay rate are orders of magnitude larger than its production rates by adipocytes and macrophages. Hence, we see a strict decreasing trend in all three mice.

Sensitivity analysis results

Figs 4 and 5 show the results of the sensitivity analysis. A positive value for a parameter means its increase will directly affect the population of cancer or the total number of cells depending on the plot title, and a negative value means the opposite. Fig 4 shows the same sensitive parameters for cancer in all three mice. However, there are some differences when it comes to sensitivity for the total number of cells.

thumbnail
Fig 4. Sensitivity results.

Sensitivity levels of the top 6 most sensitive parameters for cancer and total number of cells.

https://doi.org/10.1371/journal.pcbi.1009953.g004

thumbnail
Fig 5. Sensitivity results.

Other sensitive parameters for cancer and total number of cells.

https://doi.org/10.1371/journal.pcbi.1009953.g005

In all three mice, the natural decay rate of cancer cells is the most sensitive parameter. It is important to point out that calling it the natural decay rate of cancer is an abuse of terminology and is merely for convenience. In fact, in addition to natural death, δC includes the rate of cancer death caused by anything other than cytotoxic cells (which have been directly included in the model). This description engulfs a large set of biological reasons affecting the cancer population and is suitably recognized as the most sensitive parameter. Similarly, for λC being cancer proliferation rate promoted by anything other than adipocytes and IL6 (which have been directly included in the model). The next sensitive parameters for cancer are λCA, and λCIL6. As mentioned in the previous section, adipocytes and IL6 play a big role in cancer dynamics. Adding to that, we can see δA and λA are also among the most sensitive parameters and , and are at top of the rest of sensitive parameters, see Fig 5. These imply that controlling adipocytes and IL6 (in that order) might be a gateway to controlling cancer proliferation in these mice. The removal rate of cancer cells by cytotoxic cells is the tenth sensitive parameter. Even though this rate is directly involved in the cancer ODE, it is not as sensitive as adipocyte and IL6 related parameters mentioned before. This might be due to the very low expression of cytotoxic cells in the mouse model. Next, for mice 1 and 3, we have δM, and for mouse 2, we have as sensitive parameters. Macrophages affect the cancer population indirectly by producing cytokines like IL6, IL10, and IL2. IL6 directly affects cancer dynamics, while IL10 and IL12 do it by affecting cytotoxic cells. For Mouse 2, the parameter promotes the production of IL10, which can affect cancer through cytotoxic cells. The last set of parameters are C0 for mouse 1, for mouse 2 and for mouse 3. Cancer cells’ carrying capacity dictates how far they can go before depleting their resources and are explicitly involved in cancer ODE. Production of HMGB1 by Th can lead to cancer through several interactions, such as promoting the production of dendritic cells, which leads to more production of all cytokines, or through helper T-cells which are similarly cytokine producers. The fastest route that connects the production of IL10 by macrophages is the route that leads to cytotoxic cells. The more IL10, the more removal of cytotoxic cells and hence less death of cancer cells by cytotoxic cells.

Before discussing the sensitive parameters for total cells, it is important to note that TN is excluded from the total cell count, because they are not primarily present in the tumor microenvironment and are frequently detected in the circulation and lymphatic system [75]. Other T-cells get activated and infiltrate the tumor and can be found in copious amounts in tumors. Dendritic cells get activated inside of the tumor, and cancer cells, necrotic cells, and adipocytes are the other components of the breast tumor [140]. Finally, since most of the naive macrophages polarize into different phenotypes inside of the tumor, we include a 20% factor for MN [141]. Therefore, we have: (22) The total number of cells is an important measurement directly related to the size of the tumor. Sensitivity results show that parameters δA, λC, δN, λCA, and λA are included as top 5 sensitive parameters in all 3 mice in Fig 4. Four out of five of these parameters govern the population of adipocytes and cancer cells. These results are reasonable given that they have the largest populations among all other cells. However, the necrosis related parameter is not as straightforward, since we do not have a large number of necrotic cells in these tumors. If we track the influence of necrotic cells in the model, we see that they only contribute to the production of HMGB1. This cytokine is involved in the dynamics of helper T-cells and naive and activated dendritic cells. The sixth sensitive parameters are for mouse 1, AM for mouse 2 and for mouse 3. This difference is interesting since mouse 1 has the highest amount of naive dendritic cells, mouse 2 has the largest number of macrophages, and mouse 3 has the most regulatory T-cells. The rest of the sensitive parameters in Fig 5 are independent death rates or production rates of cells that are present in the tumor’s microenvironment and can be justified similar to above. However, we notice the presence of again. As λCA and play a crucial role in cancer proliferation, they are also important in controlling the size of the tumor.

Varying dynamics and bifurcation

This section further explores the effect of parameters on cancer dynamics. First, we perturb all sensitive parameters of cancer by 20% to see the collective effect of changing these parameters. Fig 6 shows that all three mice almost identically respond to these perturbations. This indicates good stability of the parameter estimations, especially since this has been acquired by the perturbation of parameters that have the biggest impact on the model dynamics.

thumbnail
Fig 6. Varying dynamics.

The transparent regions are acquired by 20% perturbation of all the sensitive parameters in Figs 4 and 5.

https://doi.org/10.1371/journal.pcbi.1009953.g006

Finally, we explore the local effect of the top 6 sensitive parameters (from Fig 4) for cancer on its dynamic. We do this by calculating the value of cancer at t = 42 days with respect to each sensitive parameter separately with the rest of the parameter values being fixed. As a reminder, we take 6 weeks which is the beginning of the mouse sampling, to be t = 0 days, and that makes t = 42 days corresponding to 12 weeks. As mentioned before, some state variables reach the steady-state very late; therefore, we limited the bifurcation points to the level of cancer at the last sampling time (12 weeks). The benefit of these results is that we can investigate the independent effect of large changes in single parameter values. We choose the interval [0, 0.2] for all the target parameters. Among sensitive parameters, λC has the largest estimated value of 0.0063, and the choice of the interval [0, 0.2] covers values 30 times larger than this value. Note that there is no limitation in choosing the interval, and this choice has only been made for better scaling and visual purposes (see Fig 7). In other words, the plots in Fig 7 are a zoomed-in and cropped version of more extensive bifurcation diagrams, as there are branches for negative values and chaos regions for large values that are not biologically realistic parameters’ value regimes.

thumbnail
Fig 7. Bifurcation plots.

Bifurcation of cancer values after 12 weeks with respect to the most sensitive parameters in Fig 4, namely (a)δC (b) λC (c) λCA (d) (e) δA and (f) λA.

https://doi.org/10.1371/journal.pcbi.1009953.g007

Again, all three mice show almost identical behaviors. By increasing death rate values such as δC and δA, we can significantly reduce the value of cancer cells at the last sampling time. As mentioned before, δC has a rather obscure meaning as it can be the death of cancer cells promoted by any reasons not directly included in the model. However, we can specifically see that removing adipocytes by significantly increasing their death rate leads to a notable reduction in the population of cancer cells at the last sampling time, see Fig 7a and 7e.

On the other hand, increasing the production rates such as λC, λCA, , and λA increases the cancer population at the last sampling time. Among these, λA has the smallest effect, but the others can cause the cancer population to reach double its last value in Fig 3. Also, λC has the same obscurity as δC, since it engulfs the production rate of cancer promoted by reasons other than what we have already included in the model. But λCA, and λA, specify that controlling the processes for which cancer production is promoted by IL6 and adipocytes or even reducing the production of these two can lead to a better result.

Discussion

In this study, we modeled the breast cancer progression in PyMT mice using a system of ODEs. Biologically, cancer is an intricate interaction network with many cells and molecules involved in its development. For the model, we identified key players and devised a simplified interaction network based on the available literature, see Fig 1. To further reduce the complexity of the model, we mostly used simple mass action kinetics and linear ODEs, except for cancer and adipocytes that follow a logistic growth model, see Eqs (1)–(15).

We acquired the gene expression data from the PyMT RNA-sequencing data available in the Gene Expression Omnibus (GEO) database as GSE76772 [71]. These data were collected from 3 PyMT mice after 6, 8, 10, and 12 weeks. We used CIBERSORTx B-mode to deconvolute the gene expression data and obtained each mouse’s time course data set. We used these data to estimate all the parameters in the ODE system.

Although the dynamics shown in Fig 3 are not an excellent fit for all the variables, ODE solutions closely follow a couple of important variables such as cancer and necrotic cells. The mismatch between the data and dynamics can be attributed to lacking sufficient biological information. Furthermore, we need to continue the simulations way past the feasible experimental time for all dynamics to reach the steady-state. However, we can see the steady-state values for cancer and a few more state variables by continuing the response evaluation up to 126 days (three times are longer than the last mouse data sampling).

The simulations indicate similar attributes for mice 1 and 3. Mouse 2 showed similar trends in most variables except for helper, cytotoxic, regulatory T-cells, activated dendritic cells, and IL12, see Fig 3. Maybe the most interesting observation was that cancer dynamics were almost identical in all three mice despite these differences. We argued that the similarity in adipocytes and IL6 dynamics in three mice dominates the discrepancies in other variables. This was simply a hunch based on the direct mathematical involvement of these two variables in the cancer ODE. We further confirmed this by looking at sensitivity levels of cancer to all the parameters of the model. The observed differences in variables of mouse 2 from mice 1 and 3 might be due to the differences in the percentages of macrophages’ subtypes (Fig 2), because a high level of M2 macrophages can suppress cytotoxic T cells and inhibit anti-tumor immunity [142, 143].

We carried out a sensitivity analysis based on a direct differential method, and the results showed us that the cancer population in all three mice is sensitive to δC, λC, λCA, , δA and λA in that order, see Fig 4. Commenting on the biological significance of δC and λC is rather difficult, since they can include the death and production of cancer promoted by cells and chemicals not included in the model. However, 3 out of 6 of these parameters are related to adipocytes, and looking at the rest of the sensitive parameters in Fig 5, we can also observe the importance of IL6 in cancer dynamics. We even investigated these parameters locally for much larger values through the bifurcation plots and saw regions for which cancer can be controlled through each of the sensitive parameters in Fig 5.

The link between obesity and breast cancer has been observed by many researchers. In 2007, about 7% of all new cases of cancer in women were related to obesity [144]. Obesity results in an elevated amount of adipose tissue (fat), and a direct relationship between excess fat and increased mortality rate in many types of cancer, including breast cancer, has been confirmed [145]. Prevention and medication approaches have been utilized to stop or reverse dysfunctional adipose tissue. Approaches such as weight loss strategies or medications such as metformin, statins, nonsteroidal anti-inflammatory drugs, and docosahexaenoic acid have been widely studied [146]. In addition, there are studies that suggest that leptin (a hormone produced by adipocytes) is involved in increasing breast cancer risk in postmenopausal women, and targeting it might be a key to controlling cancer in such patients [147, 148]. All of these confirm the importance of adipocytes in breast cancer development. Moreover, there are many studies recognizing IL6 as a key cytokine in progressive breast cancer, confirming that high levels of IL6 are related to poor breast cancer prognoses and showcasing its therapeutic significance in treating cancer patients [18, 19, 70]. Finally, it has been discussed that increased inflammation and IL-6 secretion in adipocytes, plus a hypoxic tumor microenvironment, creates an ideal opportunity for adipocyte-derived IL6 to promote angiogenesis [149]. So not only do adipocytes and IL6 independently contribute to poor breast cancer prognoses, but their combined effect has been acknowledged as a promoter of angiogenesis.

The approach chosen in this paper is one amongst many. There are many ways to model the interaction network. Many cells and cytokines have not been included in this study which have the potential to be integrated in our future studies. For example, our model does not consider resources such as oxygen or metabolites, macrophage heterogeneity, and the formation of blood vessels (angiogenesis). Similarly, we do not incorporate cancer stem cells, which are a tumor-initiating, self-renewing population typically resistant to therapeutics [150, 151]. The role of fibroblasts which tend to support the cancer cell niche [152], can also be considered in future iterations of the model. Furthermore, given the patient-to-patient heterogeneity, a mathematical model including more cell types and interactions mechanisms would require extensive time-course data and underlying parameters that describe these interactions. As mentioned in our manuscript, the lack of sufficient time-course data is a significant limitation of our study. Therefore, there is much room for improvement in expanding the interaction system, the validation phase, or even the dimension of the problem. Nevertheless, the current approach will guide our future studies to build targeted treatment models that focus on suppressing adipocytes and IL6. There are already studies targeting specific proteins or signaling patterns by adipocytes to control cancer [153, 154]. Also, high levels of adipocytes lead to up-regulated IL6, which build resistance to anti-VEGF therapy in breast cancer [155]. These suggest attractive therapy models with resistance terms for our future work.

Non-dimensionalization

We non-dimensionalize the system of ODEs by dividing each variable by its maximum value over all mice and time points for more stable numerical simulation, parameter estimation, and sensitivity analysis. As a result, the steady-state value of a non-dimensional variable , which is [X]/[X], equals 1. Accordingly, the following system is obtained: (23) (24) (25) (26) (27) (28) (29) (30) (31) (32) (33) (34) (35) (36) (37) Since we are not non-dimensionalizing with respect to the time, the production rates, λC and λA, and the decay rates, , , , , , δD, , δM, δC, δA, δN, δH, , , and , are left unchanged.

Parameter values

The values of the model parameters obtained from the least-squares optimization discussed in the section of Parameter estimation are reported in Table 3. The given constant values in the optimization process are C0 = 2, A0 = 2, and αNC = 1.5.

Acknowledgments

We thank the National Cancer Institute (NCI) Division of Cancer Biology and other organizers of the NCI 2020 Innovation Lab “Advancing Cancer Biology at the Frontiers of Machine Learning and Mechanistic Modeling” and the organizers of the NCI-DOE 2020 Virtual Ideas Lab, Toward Building a Cancer Patient “Digital Twin:” NCI, US Department of Energy, and the Frederick National Laboratory for Cancer Research. https://events.cancer.gov/cbiit/dtwin2020.

References

  1. 1. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2021. CA: a cancer journal for clinicians. 2021;71(1):7–33.
  2. 2. Kennecke H, Yerushalmi R, Woods R, Cheang MCU, Voduc D, Speers CH, et al. Metastatic behavior of breast cancer subtypes. Journal of clinical oncology. 2010;28(20):3271–3277. pmid:20498394
  3. 3. Engstrøm MJ, Opdahl S, Hagen AI, Romundstad PR, Akslen LA, Haugen OA, et al. Molecular subtypes, histopathological grade and survival in a historic cohort of breast cancer patients. Breast cancer research and treatment. 2013;140(3):463–473. pmid:23901018
  4. 4. Hortobagyi GN. Treatment of breast cancer. New England Journal of Medicine. 1998;339(14):974–984. pmid:9753714
  5. 5. Sharma GN, Dave R, Sanadya J, Sharma P, Sharma KK. Various types and management of breast cancer: an overview. J Adv Pharm Technol Res. 2010;1(2):109–126. pmid:22247839
  6. 6. Shahriyari L. Cell dynamics in tumour environment after treatments. Journal of the Royal Society Interface. 2017;14:20160977. pmid:28228541
  7. 7. Touboul E, Buffat L, Belkacémi Y, Lefranc JP, Uzan S, Lhuillier P, et al. Local recurrences and distant metastases after breast-conserving surgery and radiation therapy for early breast cancer. International Journal of Radiation Oncology* Biology* Physics. 1999;43(1):25–38. pmid:9989511
  8. 8. Lowery AJ, Kell MR, Glynn RW, Kerin MJ, Sweeney KJ. Locoregional recurrence after breast cancer surgery: a systematic review by receptor phenotype. Breast cancer research and treatment. 2012;133(3):831–841. pmid:22147079
  9. 9. Lim B, Hortobagyi GN. Current challenges of metastatic breast cancer. Cancer Metastasis Rev. 2016;35(4):495–514. pmid:27933405
  10. 10. Hadden M, Mittal A, Samra J, Zreiqat H, Sahni S, Ramaswamy Y. Mechanically stressed cancer microenvironment: Role in pancreatic cancer progression. Biochimica et Biophysica Acta (BBA)-Reviews on Cancer. 2020; p. 188418. pmid:32827581
  11. 11. Leong SP, Aktipis A, Maley C. Cancer initiation and progression within the cancer microenvironment. Clinical & experimental metastasis. 2018;35(5):361–367. pmid:29992410
  12. 12. Gout S, Huot J. Role of cancer microenvironment in metastasis: focus on colon cancer. Cancer Microenvironment. 2008;1(1):69–83. pmid:19308686
  13. 13. Wang J, Lei K, Han F. Tumor microenvironment: recent advances in various cancer treatments. Eur Rev Med Pharmacol Sci. 2018;22(12):3855–3864. pmid:29949179
  14. 14. Adjei IM, Blanka S. Modulation of the tumor microenvironment for cancer treatment: a biomaterials approach. Journal of functional biomaterials. 2015;6(1):81–103. pmid:25695337
  15. 15. Fokas E, McKenna WG, Muschel RJ. The impact of tumor microenvironment on cancer treatment and its modulation by direct and indirect antivascular strategies. Cancer and Metastasis Reviews. 2012;31(3):823–842. pmid:22825313
  16. 16. Ben-Baruch A. Host microenvironment in breast cancer development: inflammatory cells, cytokines and chemokines in breast cancer progression: reciprocal tumor–microenvironment interactions. Breast cancer research. 2002;5(1):1–6. pmid:12559043
  17. 17. Korkaya H, Liu S, Wicha MS, et al. Breast cancer stem cells, cytokine networks, and the tumor microenvironment. The Journal of clinical investigation. 2011;121(10):3804–3809. pmid:21965337
  18. 18. Heo TH, Wahler J, Suh N. Potential therapeutic implications of IL-6/IL-6R/gp130-targeting agents in breast cancer. Oncotarget. 2016;7(13):15460. pmid:26840088
  19. 19. Goldberg EJ, Schwertfeger LK. Proinflammatory cytokines in breast cancer: mechanisms of action and potential targets for therapeutics. Current drug targets. 2010;11(9):1133–1146.
  20. 20. Matsuo Y, Takeyama H, Guha S. Cytokine network: new targeted therapy for pancreatic cancer. Current pharmaceutical design. 2012;18(17):2416–2419. pmid:22372505
  21. 21. Tagawa M. Cytokine therapy for cancer. Current pharmaceutical design. 2000;6(6):681–699. pmid:10788604
  22. 22. Lampreht Tratar U, Horvat S, Cemazar M. Transgenic mouse models in cancer research. Frontiers in oncology. 2018;8:268. pmid:30079312
  23. 23. Hursting SD, Slaga TJ, Fischer SM, DiGiovanni J, Phang JM. Mechanism-based cancer prevention approaches: targets, examples, and the use of transgenic mice. Journal of the National Cancer Institute. 1999;91(3):215–225. pmid:10037099
  24. 24. Hurwitz AA, Foster BA, Kwon ED, Truong T, Choi EM, Greenberg NM, et al. Combination immunotherapy of primary prostate cancer in a transgenic mouse model using CTLA-4 blockade. Cancer research. 2000;60(9):2444–2448. pmid:10811122
  25. 25. Gingrich J, Greenberg N. A transgenic mouse prostate cancer model. Toxicologic pathology. 1996;24(4):502–504. pmid:8864193
  26. 26. Korac-Prlic J, Degoricija M, Vilović K, Haupt B, Ivanišević T, Franković L, et al. Targeting Stat3 signaling impairs the progression of bladder cancer in a mouse model. Cancer Letters. 2020;490:89–99. pmid:32659249
  27. 27. Floc’h N, Kinkade CW, Kobayashi T, Aytes A, Lefebvre C, Mitrofanova A, et al. Dual targeting of the Akt/mTOR signaling pathway inhibits castration-resistant prostate cancer in a genetically engineered mouse model. Cancer research. 2012;72(17):4483–4493. pmid:22815528
  28. 28. Zhao H, Richardson R, Talebloo N, Mukherjee P, Wang P, Moore A. uMUC1-targeting magnetic resonance imaging of therapeutic response in an orthotropic mouse model of colon cancer. Molecular imaging and biology. 2019;21(5):852–860. pmid:30793239
  29. 29. Zeng H, Wei W, Xu X. Chemokine (CXC motif) receptor 4 RNA interference inhibits bone metastasis in breast cancer. Oncology letters. 2014;8(1):77–81. pmid:24959222
  30. 30. Chang A, Le CP, Walker AK, Creed SJ, Pon CK, Albold S, et al. β2-Adrenoceptors on tumor cells play a critical role in stress-enhanced metastasis in a mouse model of breast cancer. Brain, behavior, and immunity. 2016;57:106–115. pmid:27321906
  31. 31. Ling X, Arlinghaus RB. Knockdown of STAT3 expression by RNA interference inhibits the induction of breast tumors in immunocompetent mice. Cancer research. 2005;65(7):2532–2536. pmid:15805244
  32. 32. Attalla S, Taifour T, Bui T, Muller W. Insights from transgenic mouse models of PyMT-induced breast cancer: recapitulating human breast cancer progression in vivo. Oncogene. 2021;40(3):475–491. pmid:33235291
  33. 33. Guy CT, Cardiff R, Muller WJ. Induction of mammary tumors by expression of polyomavirus middle T oncogene: a transgenic mouse model for metastatic disease. Molecular and cellular biology. 1992;12(3):954–961. pmid:1312220
  34. 34. Mohammad Mirzaei N, Weintraub WS, Fok PW. An integrated approach to simulating the vulnerable atherosclerotic plaque. American Journal of Physiology-Heart and Circulatory Physiology. 2020;319(4):H835–H846. pmid:32795179
  35. 35. Hao W, Friedman A. The LDL-HDL profile determines the risk of atherosclerosis: a mathematical model. PloS one. 2014;9(3):e90497. pmid:24621857
  36. 36. Chalmers AD, Cohen A, Bursill CA, Myerscough MR. Bifurcation and dynamics in a mathematical model of early atherosclerosis. Journal of mathematical biology. 2015;71(6):1451–1480. pmid:25732771
  37. 37. Baba IA, Yusuf A, Nisar KS, Abdel-Aty AH, Nofal TA. Mathematical model to assess the imposition of lockdown during COVID-19 pandemic. Results in Physics. 2021;20:103716. pmid:33520624
  38. 38. Chitnis N, Hyman JM, Cushing JM. Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model. Bulletin of mathematical biology. 2008;70(5):1272. pmid:18293044
  39. 39. Misra A, Singh V. A delay mathematical model for the spread and control of water borne diseases. Journal of theoretical biology. 2012;301:49–56. pmid:22586723
  40. 40. Woittiez R, Huijing P, Boom H, Rozendal R. A three-dimensional muscle model: a quantified relation between form and function of skeletal muscles. Journal of Morphology. 1984;182(1):95–113. pmid:6492171
  41. 41. Van den Bogert AJ, Geijtenbeek T, Even-Zohar O, Steenbrink F, Hardin EC. A real-time system for biomechanical analysis of human movement and muscle function. Medical & biological engineering & computing. 2013;51(10):1069–1077. pmid:23884905
  42. 42. Luke RA, Braun RJ, Driscoll TA, Begley CG, Awisi-Gyau D. Parameter estimation for evaporation-driven tear film thinning. Bulletin of Mathematical Biology. 2020;82:1–41. pmid:32506271
  43. 43. Zhu H, Chauhan A. A mathematical model for ocular tear and solute balance. Current eye research. 2005;30(10):841–854. pmid:16251121
  44. 44. Le T, Su S, Kirshtein A, Shahriyari L. Data-Driven Mathematical Model of Osteosarcoma. Cancers. 2021;13(10):2367. pmid:34068946
  45. 45. Kirshtein A, Akbarinejad S, Hao W, Le T, Su S, Aronow RA, et al. Data Driven Mathematical Model of Colon Cancer Progression. Journal of Clinical Medicine. 2020;9(12):3947. pmid:33291412
  46. 46. Rhodes A, Hillen T. A mathematical model for the immune-mediated theory of metastasis. Journal of theoretical biology. 2019;482:109999. pmid:31493486
  47. 47. Liao KL, Bai XF, Friedman A. The role of CD200–CD200R in tumor immune evasion. Journal of theoretical biology. 2013;328:65–76. pmid:23541619
  48. 48. Katira P, Bonnecaze RT, Zaman MH. Modeling the mechanics of cancer: effect of changes in cellular and extra-cellular mechanical properties. Frontiers in oncology. 2013;3:145. pmid:23781492
  49. 49. Mpekris F, Voutouri C, Papageorgis P, Stylianopoulos T. Stress alleviation strategy in cancer treatment: Insights from a mathematical model. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik. 2018;98(12):2295–2306.
  50. 50. Sanga S, Sinek JP, Frieboes HB, Ferrari M, Fruehauf JP, Cristini V. Mathematical modeling of cancer progression and response to chemotherapy. Expert review of anticancer therapy. 2006;6(10):1361–1376. pmid:17069522
  51. 51. Jain HV, Clinton SK, Bhinder A, Friedman A. Mathematical modeling of prostate cancer progression in response to androgen ablation therapy. Proceedings of the National Academy of Sciences. 2011;108(49):19701–19706. pmid:22106268
  52. 52. Le T, Su S, Shahriyari L. Investigating Optimal Chemotherapy Options for Osteosarcoma Patients through a Mathematical Model. Cells. 2021;10(8):2009. pmid:34440778
  53. 53. Shahriyari L, Komarova NL. Symmetric vs. asymmetric stem cell divisions: an adaptation against cancer? PloS one. 2013;8(10):e76195. pmid:24204602
  54. 54. Shahriyari L, Komarova NL. The role of the bi-compartmental stem cell niche in delaying cancer. Physical Biology. 2015;12(5):055001. pmid:26228740
  55. 55. Shahriyari L, Mahdipour-Shirayeh A. Modeling dynamics of mutants in heterogeneous stem cell niche. Physical Biology. 2017;14(1). pmid:28102174
  56. 56. Xu Y. A free boundary problem model of ductal carcinoma in situ. Discrete & Continuous Dynamical Systems—B. 2004;4(1):337–348.
  57. 57. Eladdadi A, Isaacson D. A mathematical model for the effects of HER2 overexpression on cell proliferation in breast cancer. Bulletin of mathematical biology. 2008;70(6):1707. pmid:18648887
  58. 58. Mohammad Mirzaei N, Su S, Sofia D, Hegarty M, Abdel-Rahman MH, Asadpoure A, et al. A Mathematical Model of Breast Tumor Progression Based on Immune Infiltration. Journal of Personalized Medicine. 2021;11(10):1031. pmid:34683171
  59. 59. Eikenberry SE, Nagy JD, Kuang Y. The evolutionary impact of androgen levels on prostate cancer in a multi-scale mathematical model. Biology direct. 2010;5(1):1–28.
  60. 60. Budithi A, Su S, Kirshtein A, Shahriyari L. Data Driven Mathematical Model of FOLFIRI Treatment for Colon Cancer. Cancers. 2021;13(11):2632. pmid:34071939
  61. 61. Norton KA, Jin K, Popel AS. Modeling triple-negative breast cancer heterogeneity: Effects of stromal macrophages, fibroblasts and tumor vasculature. Journal of theoretical biology. 2018;452:56–68. pmid:29750999
  62. 62. Butner JD, Fuentes D, Ozpolat B, Calin GA, Zhou X, Lowengrub J, et al. A multiscale agent-based model of ductal carcinoma in situ. IEEE Transactions on Biomedical Engineering. 2019;67(5):1450–1461. pmid:31603768
  63. 63. Hao W, Friedman A. Serum upar as biomarker in breast cancer recurrence: A mathematical model. PLoS One. 2016;11(4):e0153508. pmid:27078836
  64. 64. Vaghi C, Rodallec A, Fanciullino R, Ciccolini J, Mochel JP, Mastri M, et al. Population modeling of tumor growth curves and the reduced Gompertz model improve prediction of the age of experimental tumors. PLoS Comput Biol. 2020;16(2):e1007178. pmid:32097421
  65. 65. Strogatz SH. Nonlinear dynamics and chaos with student solutions manual: With applications to physics, biology, chemistry, and engineering. CRC press; 2018.
  66. 66. Lai X, Stiff A, Duggan M, Wesolowski R, Carson WE, Friedman A. Modeling combination therapy for breast cancer with BET and immune checkpoint inhibitors. Proceedings of the National Academy of Sciences of the United States of America;115(21):5534–5539. pmid:29735668
  67. 67. Apetoh L, Ghiringhelli F, Tesniere A, Criollo A, Ortiz C, Lidereau R, et al. The interaction between HMGB1 and TLR4 dictates the outcome of anticancer chemotherapy and radiotherapy. Immunological Reviews. 2007;220:47–59. pmid:17979839
  68. 68. Sun Y, Liu L, Zhou L, Yu S, Lan Y, Liang Q, et al. Tumor Microenvironment-Triggered Charge Reversal Polymetformin-Based Nanosystem Co-Delivered Doxorubicin and IL-12 Cytokine Gene for Chemo–Gene Combination Therapy on Metastatic Breast Cancer. ACS Applied Materials & Interfaces. 2020;12(41):45873–45890. pmid:32924511
  69. 69. Masjedi A, Hashemi V, Hojjat-Farsangi M, Ghalamfarsa G, Azizi G, Yousefi M, et al. The significant role of interleukin-6 and its signaling pathway in the immunopathogenesis and treatment of breast cancer. Biomedicine & Pharmacotherapy. 2018;108:1415–1424. pmid:30372844
  70. 70. Guo Y, Xu F, Lu T, Duan Z, Zhang Z. Interleukin-6 signaling pathway in targeted therapy for cancer. Cancer treatment reviews. 2012;38(7):904–910. pmid:22651903
  71. 71. Cai Y, Nogales-Cadenas R, Zhang Q, Lin JR, Zhang W, O’Brien K, et al. Transcriptomic dynamics of breast cancer progression in the MMTV-PyMT mouse model. Bmc Genomics. 2017;18(1):1–14. pmid:28212608
  72. 72. Chang LY, Pollard NS. Constrained least-squares optimization for robust estimation of center of rotation. Journal of biomechanics. 2007;40(6):1392–1400. pmid:16824530
  73. 73. Bø TH, Dysvik B, Jonassen I. LSimpute: accurate estimation of missing values in microarray data with least squares methods. Nucleic acids research. 2004;32(3):e34–e34. pmid:14978222
  74. 74. Kim H, Park H. Sparse non-negative matrix factorizations via alternating non-negativity-constrained least squares for microarray data analysis. Bioinformatics. 2007;23(12):1495–1502. pmid:17483501
  75. 75. Punt J. Chapter 4—Adaptive Immunity: T Cells and Cytokines. In: Prendergast GC, Jaffee EM, editors. Cancer Immunotherapy (Second Edition). second edition ed. San Diego: Academic Press; 2013. p. 41–53. Available from: https://www.sciencedirect.com/science/article/pii/B978012394296800004X.
  76. 76. Palucka K, Coussens LM, O’Shaughnessy J. Dendritic cells, inflammation and breast cancer. Cancer journal (Sudbury, Mass). 2013;19(6). pmid:24270350
  77. 77. Tay ASMS, Amano T, Edwards LA, John SY. CD133 mRNA transfected dendritic cells induces a coordinated cytotoxic and helper T cell responses against breast cancer stem cells. Molecular Therapy-Oncolytics. 2021. pmid:34485687
  78. 78. Dumitriu IE, Baruah P, Valentinis B, Voll RE, Herrmann M, Nawroth PP, et al. Release of high mobility group box 1 by dendritic cells controls T cell activation via the receptor for advanced glycation end products. The Journal of Immunology. 2005;174(12):7506–7515. pmid:15944249
  79. 79. Zhang Y, Yao Ym, Huang Lf, Dong N, Yu Y, Sheng Zy. The potential effect and mechanism of high-mobility group box 1 protein on regulatory T cell-mediated immunosuppression. Journal of Interferon & Cytokine Research. 2011;31(2):249–257. pmid:21087077
  80. 80. Bell RB, Feng Z, Bifulco CB, Leidner R, Weinberg A, Fox BA. In: 15—Immunotherapy. Elsevier; 2018. p. 314–340.
  81. 81. Moon BI, Kim TH, Seoh JY. Functional Modulation of Regulatory T Cells by IL-2. PLoS One. 2015;10(11):e0141864. pmid:26529512
  82. 82. Watanabe MAE, Oda JMM, Amarante MK, Voltarelli JC. Regulatory T cells and breast cancer: implications for immunopathogenesis. Cancer and Metastasis Reviews. 2010;29(4):569–579. pmid:20830504
  83. 83. Wang K, Vella AT. Regulatory T Cells and Cancer: A Two-Sided Story. Immunological Investigations. 2016;45(8):797–812. pmid:27603750
  84. 84. Sheikhpour E, Noorbakhsh P, Foroughi E, Farahnak S, Nasiri R, Neamatzadeh H. A survey on the role of interleukin-10 in breast cancer: A narrative. Reports of biochemistry & molecular biology. 2018;7(1):30. pmid:30324115
  85. 85. Fu C, Jiang A. Dendritic Cells and CD8 T Cell Immunity in Tumor Microenvironment. Frontiers in Immunology. 2018;9. pmid:30619378
  86. 86. Ma Y, Shurin GV, Peiyuan Z, Shurin MR. Dendritic cells in the cancer microenvironment. Journal of Cancer. 2013;4(1):36–44. pmid:23386903
  87. 87. Tran Janco JM, Lamichhane P, Karyampudi L, Knutson KL. Tumor-Infiltrating Dendritic Cells in Cancer Pathogenesis. The Journal of Immunology. 2015;194(7):2985–2991. pmid:25795789
  88. 88. Jonuleit H, Schmitt E, Steinbrink K, Enk AH. Dendritic cells as a tool to induce anergic and regulatory T cells. Trends in immunology. 2001;22(7):394–400. pmid:11429324
  89. 89. Aspord C, Pedroza-Gonzalez A, Gallegos M, Tindle S, Burton EC, Su D, et al. Breast cancer instructs dendritic cells to prime interleukin 13–secreting CD4+ T cells that facilitate tumor development. Journal of Experimental Medicine. 2007;204(5):1037–1047. pmid:17438063
  90. 90. Tang D, Kang R, Zeh HJ, Lotze MT. High-mobility Group Box 1 [HMGB1] and Cancer. Biochimica et biophysica acta;1799(1):131. pmid:20123075
  91. 91. Lee SY, Ju MK, Jeon HM, Jeong EK, Lee YJ, Kim CH, et al. Regulation of Tumor Progression by Programmed Necrosis. Oxidative Medicine and Cellular Longevity. pmid:29636841
  92. 92. Zhang Y, Liu Z, Hao X, Li A, Zhang J, Carey CD, et al. Tumor-derived high-mobility group box 1 and thymic stromal lymphopoietin are involved in modulating dendritic cells to activate T regulatory cells in a mouse model. Cancer Immunology, Immunotherapy. 2018;67(3):353–366. pmid:29116372
  93. 93. Aras S, Zaidi MR. TAMeless traitors: macrophages in cancer progression and metastasis. British Journal of Cancer. 2017;117(11):1583–1591. pmid:29065107
  94. 94. Sica A, Saccani A, Mantovani A. Tumor-associated macrophages: a molecular perspective. International Immunopharmacology. 2002;2(8):1045–54. pmid:12349942
  95. 95. Yang C, He L, He P, Liu Y, Wang W, He Y, et al. Increased drug resistance in breast cancer by tumor-associated macrophages through IL-10/STAT3/bcl-2 signaling pathway. Medical oncology. 2015;32(2):14.
  96. 96. Obeid E, Nanda R, Fu YX, Olopade OI. The role of tumor-associated macrophages in breast cancer progression (review). International Journal of Oncology. 2013;43(1):5–12. pmid:23673510
  97. 97. 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. Journal of Cellular Biochemistry. 2017;118(9):2484–2501. pmid:28106295
  98. 98. Chanmee T, Ontong P, Konno K, Itano N. Tumor-Associated Macrophages as Major Players in the Tumor Microenvironment. Cancers. 2014;6(3):1670–1690. pmid:25125485
  99. 99. Wu Q, Li B, Li Z, Li J, Sun S, Sun S. Cancer-associated adipocytes: key players in breast cancer progression. Journal of Hematology & Oncology; 12. pmid:31500658
  100. 100. Mao Y, Keller ET, Garfield DH, Shen K, Wang J. Stroma Cells in Tumor Microenvironment and Breast Cancer. Cancer metastasis reviews;32(0):303–315. pmid:23114846
  101. 101. Schafer ZT, Brugge JS, et al. IL-6 involvement in epithelial cancers. The Journal of clinical investigation. 2007;117(12):3660–3663. pmid:18060028
  102. 102. Chu DT, Phuong TNT, Tien NLB, Tran DK, Nguyen TT, Thanh VV, et al. The Effects of Adipocytes on the Regulation of Breast Cancer in the Tumor Microenvironment: An Update. Cells. 2019;8(8). pmid:31398937
  103. 103. Liu E, Samad F, Mueller BM. Local adipocytes enable estrogen-dependent breast cancer growth: Role of leptin and aromatase. Adipocyte. 2013;2(3):165–169. pmid:23991363
  104. 104. Neel JC, Humbert L, Lebrun JJ. The dual role of TGFβ in human cancer: from tumor suppression to cancer metastasis. International Scholarly Research Notices. 2012;2012.
  105. 105. Huang Y, Ma C, Zhang Q, Ye J, Wang F, Zhang Y, et al. CD4+ and CD8+ T cells have opposing roles in breast cancer progression and outcome. Oncotarget. 2015;6(19):17462. pmid:25968569
  106. 106. Dirat B, Bochet L, Dabek M, Daviaud D, Dauvillier S, Majed B, et al. Cancer-associated adipocytes exhibit an activated phenotype and contribute to breast cancer invasion. Cancer research. 2011;71(7):2455–2465. pmid:21459803
  107. 107. Leek R, Landers R, Harris A, Lewis C. Necrosis correlates with high vascular density and focal macrophage infiltration in invasive carcinoma of the breast. British journal of cancer. 1999;79(5):991–995. pmid:10070902
  108. 108. Xue J, Suarez JS, Minaai M, Li S, Gaudino G, Pass HI, et al. HMGB1 as a therapeutic target in disease. Journal of cellular physiology. 2021;236(5):3406–3419. pmid:33107103
  109. 109. Li G, Liang X, Lotze MT. HMGB1: The Central Cytokine for All Lymphoid Cells. Frontiers in Immunology. 2013;4. pmid:23519706
  110. 110. Bianchi ME, Manfredi AA. High-mobility group box 1 (HMGB1) protein at the crossroads between innate and adaptive immunity. Immunological reviews. 2007;220(1):35–46. pmid:17979838
  111. 111. Wang S, Zhang Y. HMGB1 in inflammation and cancer. Journal of Hematology & Oncology;13(1):116. pmid:32831115
  112. 112. Okuhira K, Demizu Y, Hattori T, Ohoka N, Shibata N, Nishimaki-Mogami T, et al. Development of hybrid small molecules that induce degradation of estrogen receptor-alpha and necrotic cell death in breast cancer cells. Cancer science. 2013;104(11):1492–1498. pmid:23992566
  113. 113. Bonaldi T, Talamo F, Scaffidi P, Ferrera D, Porto A, Bachi A, et al. Monocytic cells hyperacetylate chromatin protein HMGB1 to redirect it towards secretion. The EMBO Journal. 2003;22(20):5551–5560. pmid:14532127
  114. 114. Tang D, Shi Y, Kang R, Li T, Xiao W, Wang H, et al. Hydrogen peroxide stimulates macrophages and monocytes to actively release HMGB1. Journal of Leukocyte Biology. 2007;81(3):741–747. pmid:17135572
  115. 115. Wang H, Bloom O, Zhang M, Vishnubhakat JM, Ombrellino M, Che J, et al. HMG-1 as a late mediator of endotoxin lethality in mice. Science. 1999;285(5425):248–251. pmid:10398600
  116. 116. 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(2):609–616. pmid:15802534
  117. 117. Gougeon ML, Bras M. Natural killer cells, dendritic cells, and the alarmin high-mobility group box 1 protein: a dangerous trio in HIV-1 infection? Current Opinion in HIV and AIDS. 2011;6(5):364–372. pmid:21825870
  118. 118. DeMarco RA, Fink MP, Lotze MT. Monocytes promote natural killer cell interferon gamma production in response to the endogenous danger signal HMGB1. Molecular immunology. 2005;42(4):433–444. pmid:15607795
  119. 119. Kato I, Tanaka K, Yokokura T. Lactic acid bacterium potently induces the production of interleukin-12 and interferon-γ by mouse splenocytes. International journal of immunopharmacology. 1999;21(2):121–131. pmid:10230875
  120. 120. Segovia-Mendoza M, Morales-Montor J. Immune Tumor Microenvironment in Breast Cancer and the Participation of Estrogen and Its Receptors in Cancer Physiopathology. Frontiers in Immunology;10. pmid:30881360
  121. 121. Fan X, Zhang H, Cheng Y, Jiang X, Zhu J, Jin T. Double roles of macrophages in human neuroimmune diseases and their animal models. Mediators of inflammation. 2016;2016. pmid:27034594
  122. 122. Khan MM. Role of cytokines. In: Immunopharmacology. Springer; 2016. p. 57–92.
  123. 123. Hart AL, Al-Hassi HO, Rigby RJ, Bell SJ, Emmanuel AV, Knight SC, et al. Characteristics of intestinal dendritic cells in inflammatory bowel diseases. Gastroenterology. 2005;129(1):50–65. pmid:16012934
  124. 124. Iwasaki A, Kelsall BL. Freshly isolated Peyer’s patch, but not spleen, dendritic cells produce interleukin 10 and induce the differentiation of T helper type 2 cells. The Journal of experimental medicine. 1999;190(2):229–240. pmid:10432286
  125. 125. Ohyagi H, Onai N, Sato T, Yotsumoto S, Liu J, Akiba H, et al. Monocyte-derived dendritic cells perform hemophagocytosis to fine-tune excessive immune responses. Immunity. 2013;39(3):584–598. pmid:24035363
  126. 126. Couper KN, Blount DG, Riley EM. IL-10: the master regulator of immunity to infection. Journal of Immunology. 2008;180(9):5771–5777. pmid:18424693
  127. 127. Trinchieri G. Interleukin-10 production by effector T cells: Th1 cells show self control. Journal of Experimental Medicine. 2007;204(2):239–243. pmid:17296790
  128. 128. Lee YH, Leong WY, Wilder-Smith A. Markers of dengue severity: a systematic review of cytokines and chemokines. Journal of General Virology. 2016;97(12):3103–3119. pmid:27902364
  129. 129. Loges S, Schmidt T, Tjwa M, Van Geyte K, Lievens D, Lutgens E, et al. Malignant cells fuel tumor growth by educating infiltrating leukocytes to produce the mitogen Gas6. Blood, The Journal of the American Society of Hematology. 2010;115(11):2264–2273. pmid:19965679
  130. 130. Choi J, Cha YJ, Koo JS. Adipocyte biology in breast cancer: From silent bystander to active facilitator. Progress in lipid research. 2018;69:11–20. pmid:29175445
  131. 131. Liu T, Clark RK, McDonnell PC, Young PR, White RF, Barone FC, et al. Tumor necrosis factor-alpha expression in ischemic neurons. Stroke. 1994;25(7):1481–1488. pmid:8023366
  132. 132. Erez N, Glanz S, Raz Y, Avivi C, Barshack I. Cancer associated fibroblasts express pro-inflammatory factors in human breast and ovarian tumors. Biochemical and biophysical research communications. 2013;437(3):397–402. pmid:23831470
  133. 133. Newman AM, Steen CB, Liu CL, Gentles AJ, Chaudhuri AA, Scherer F, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nature biotechnology. 2019;37(7):773–782. pmid:31061481
  134. 134. Kim Y, Wallace J, Li F, Ostrowski M, Friedman A. Transformed epithelial cells and fibroblasts/myofibroblasts interaction in breast tumor: a mathematical model and experiments. Journal of mathematical biology. 2010;61(3):401–421. pmid:19902212
  135. 135. Kim H, Golub GH, Park H. Missing value estimation for DNA microarray gene expression data: local least squares imputation. Bioinformatics. 2005;21(2):187–198. pmid:15333461
  136. 136. Liu YF, Xu S, Gong H, Cui YF, Song DD, Zhou YP. Partial least-squares discriminant analysis optimized by particle swarm optimization: application to 1H nuclear magnetic resonance analysis of lung cancer metabonomics. Journal of Chemometrics. 2015;29(10):537–546.
  137. 137. Rabitz H, Kramer M, Dacol D. Sensitivity analysis in chemical kinetics. Annual review of physical chemistry. 1983;34(1):419–461.
  138. 138. Gerstner T, Griebel M. Numerical integration using sparse grids. Numerical algorithms. 1998;18(3):209–232.
  139. 139. Heiss F, Winschel V. Likelihood approximation by numerical integration on sparse grids. journal of Econometrics. 2008;144(1):62–80.
  140. 140. Soysal SD, Tzankov A, Muenst SE. Role of the tumor microenvironment in breast cancer. Pathobiology. 2015;82(3-4):142–152. pmid:26330355
  141. 141. Qiu SQ, Waaijer SJ, Zwager MC, de Vries EG, van der Vegt B, Schröder CP. Tumor-associated macrophages in breast cancer: Innocent bystander or important player? Cancer treatment reviews. 2018;70:178–189. pmid:30227299
  142. 142. Ruffell B, Chang-Strachan D, Chan V, Rosenbusch A, Ho CM, Pryer N, et al. Macrophage IL-10 blocks CD8+ T cell-dependent responses to chemotherapy by suppressing IL-12 expression in intratumoral dendritic cells. Cancer cell. 2014;26(5):623–637. pmid:25446896
  143. 143. Peranzoni E, Lemoine J, Vimeux L, Feuillet V, Barrin S, Kantari-Mimoun C, et al. Macrophages impede CD8 T cells from reaching tumor cells and limit the efficacy of anti–PD-1 treatment. Proceedings of the National Academy of Sciences. 2018;115(17):E4041–E4050. pmid:29632196
  144. 144. Polednak AP. Estimating the number of US incident cancers attributable to obesity and the impact on temporal trends in incidence rates for obesity-related cancers. Cancer detection and prevention. 2008;32(3):190–199. pmid:18790577
  145. 145. Eheman C, Henley SJ, Ballard-Barbash R, Jacobs EJ, Schymura MJ, Noone AM, et al. Annual report to the nation on the status of cancer, 1975-2008, featuring cancers associated with excess weight and lack of sufficient physical activity; 2012.
  146. 146. Gucalp A, Iyengar NM, Hudis CA, Dannenberg AJ. Targeting obesity-related adipose tissue dysfunction to prevent cancer development and progression. In: Seminars in oncology. vol. 43. Elsevier; 2016. p. 154–160.
  147. 147. Surmacz E. Obesity hormone leptin: a new target in breast cancer? Breast Cancer Research. 2007;9(1):1–2. pmid:17274833
  148. 148. Garofalo C, Surmacz E. Leptin and cancer. Journal of cellular physiology. 2006;207(1):12–22. pmid:16110483
  149. 149. Gyamfi J, Eom M, Koo JS, Choi J. Multifaceted roles of interleukin-6 in adipocyte–breast cancer cell interaction. Translational Oncology. 2018;11(2):275–285. pmid:29413760
  150. 150. Guo W. Concise review: breast cancer stem cells: regulatory networks, stem cell niches, and disease relevance. Stem cells translational medicine. 2014;3(8):942–948. pmid:24904174
  151. 151. Guo W, Keckesova Z, Donaher JL, Shibue T, Tischler V, Reinhardt F, et al. Slug and Sox9 cooperatively determine the mammary stem cell state. Cell. 2012;148(5):1015–1028. pmid:22385965
  152. 152. Malanchi I, Santamaria-Martínez A, Susanto E, Peng H, Lehr HA, Delaloye JF, et al. Interactions between cancer stem cells and their niche govern metastatic colonization. Nature. 2012;481(7379):85–89.
  153. 153. Gernapudi R, Yao Y, Zhang Y, Wolfson B, Roy S, Duru N, et al. Targeting exosomes from preadipocytes inhibits preadipocyte to cancer stem cell signaling in early-stage breast cancer. Breast cancer research and treatment. 2015;150(3):685–695. pmid:25783182
  154. 154. Grisouard J, Dembinski K, Mayer D, Keller U, Müller B, Christ-Crain M. Targeting AMP-activated protein kinase in adipocytes to modulate obesity-related adipokine production associated with insulin resistance and breast cancer cell proliferation. Diabetology & metabolic syndrome. 2011;3(1):1–7. pmid:21774820
  155. 155. Incio J, Ligibel JA, McManus DT, Suboj P, Jung K, Kawaguchi K, et al. Obesity promotes resistance to anti-VEGF therapy in breast cancer by up-regulating IL-6 and potentially FGF-2. Science translational medicine. 2018;10(432). pmid:29540614