key: cord-0029199-13evk15i authors: Mohammad Mirzaei, Navid; Changizi, Navid; Asadpoure, Alireza; Su, Sumeyye; Sofia, Dilruba; Tatarova, Zuzana; Zervantonakis, Ioannis K.; Chang, Young Hwan; Shahriyari, Leili title: Investigating key cell types and molecules dynamics in PyMT mice model of breast cancer through a mathematical model date: 2022-03-16 journal: PLoS Comput Biol DOI: 10.1371/journal.pcbi.1009953 sha: 3093eede34755e45e5ce4f3b2a71e4f7816a541f doc_id: 29199 cord_uid: 13evk15i 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. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 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, basallike, 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 [4] [5] [6] . 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 [10] [11] [12] or its response to treatments [13] [14] [15] . 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 [18] [19] [20] [21] . 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 [22] [23] [24] [25] , gene targeting [26] [27] [28] , RNA interference [29] [30] [31] , 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 [34] [35] [36] , disease spread [37] [38] [39] , 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 [44] [45] [46] [47] , some investigate the mechanical properties of a cancerous tissue [48, 49] , and some others focus on modeling cancer response to different treatments [50] [51] [52] . The most desirable features of mathematical models of cancer, including stochastic [53-55] and deterministic models [56] [57] [58] , 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 longterm 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 [66 -70] . 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 [72] [73] [74] . 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. 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. T-cells. We categorize T-cells into 4 major sub-types: naive (T N ), helper (T h ), cytotoxic (T C ) and regulatory (T r ) 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. 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 Cytotoxic T-Cells (T c ). 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. Regulatory T-Cells (T r ). 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. Naive T-Cells (T N ). Combining Eqs (1)-(3) for the activation of naive T-cells and adding an independent naive T-cells production rate A T N , we get the following ODE for for naive T-cells Cancer cells [66, 89] and HMGB1 [67, [90] [91] [92] can activate dendritic cells. Moreover, cancer cells may promote natural death of dendritic cells in different ways [86] . Adding A D N as the production rate of naive dendritic cells, we get the following system of equations for naive dendritic cells (D N ) and activated dendritic cells (D) 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 Denoting naive macrophages by M N , activated macrophages by M, and the production rate of naive macrophages by A M , we can write the following system of equations for the dynamics of naive and activated macrophages. 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]/C 0 ), where C 0 is the maximum capacity. In addition, IL-6 promotes the proliferation of cancer cells [99] [100] [101] . 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. 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. 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. 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, [113] [114] [115] , natural killer (NK) cells which behave like cytotoxic T-cells [109, [116] [117] [118] , and cancer cells [66, 90, 91] . Therefore, the dynamics of HMGB1 is modeled by the following equation. IL-12 (IL 12 ). IL-12, stimulates the differentiation of naive T-cells into helper T-cells. Macrophages and dendritic cells secrete 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. IL-10 (IL 10 ). IL-10 is secreted by macrophages [93, 121, 122] , dendritic cells [86, [123] [124] [125] , 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. IL-6 (IL 6 ). 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]. 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 men- [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 t i 2 {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 191 192 is just the simplified ratio of cancer cells to necrotic cells 0.955 : 0.005. 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 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 C 0 and A 0 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 Minimize : 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 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 x to show non-dimensional variables. For a generic ODE system of the form where x ¼ hx 1 ; � � � ; x 15 i and θ = hθ 1 , � � �, θ 57 i are vectors of state variables and parameters of our model, respectively, the first order local sensitivity of a variable x j with respect to a parameter θ i is evaluated by To obtain the sensitivity vector, s = hs 1 , � � �, s M i, we use a direct differential method. That is, we differentiate Eq (18) with respect to θ i to get d dt and then we use a forward Euler discretization in time for Eq (20) to find s i ¼ @x=@y i . The sensitivity of each parameter in the neighborhood of a chosen parameter set O(θ) is then defined asŝ 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 steadystate 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. 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. 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 Investigating key cell types and molecules dynamics in PyMT mice model of breast cancer 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 Tcells, 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. parameters for cancer in all three mice. However, there are some differences when it comes to sensitivity for the total number of cells. 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 l IL 6 A ; d IL 6 , and l IL 6 M 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 l IL 10 C 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 l IL 10 C promotes the production of IL10, which can affect cancer through cytotoxic cells. The last set of parameters are C 0 for mouse 1, l T h H for mouse 2 and l MIL 10 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 T h 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 T N 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 M N [141] . Therefore, we have: 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 A D N for mouse 1, A M for mouse 2 and d T r 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 l CIL 6 again. As λ CA and l CIL 6 play a crucial role in cancer proliferation, they are also important in controlling the size of the tumor. 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. 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 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 , l CIL 6 , 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 , l CIL 6 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. 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 , l CIL 6 , δ 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 upregulated IL6, which build resistance to anti-VEGF therapy in breast cancer [155] . These suggest attractive therapy models with resistance terms for our future work. 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 ½X�, which is [X]/[X 1 ], equals 1. Accordingly, the following system is obtained: Since we are not non-dimensionalizing with respect to the time, the production rates, λ C and λ A , and the decay rates, 10 , and d IL 6 , are left unchanged. 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 C 0 = 2, A 0 = 2, and α NC = 1.5. CA: a cancer journal for clinicians Metastatic behavior of breast cancer subtypes Treatment of breast cancer Various types and management of breast cancer: an overview Cell dynamics in tumour environment after treatments Local recurrences and distant metastases after breast-conserving surgery and radiation therapy for early breast cancer Locoregional recurrence after breast cancer surgery: a systematic review by receptor phenotype. Breast cancer research and treatment Current challenges of metastatic breast cancer Mechanically stressed cancer microenvironment: Role in pancreatic cancer progression Cancer initiation and progression within the cancer microenvironment Role of cancer microenvironment in metastasis: focus on colon cancer Tumor microenvironment: recent advances in various cancer treatments Modulation of the tumor microenvironment for cancer treatment: a biomaterials approach The impact of tumor microenvironment on cancer treatment and its modulation by direct and indirect antivascular strategies Host microenvironment in breast cancer development: inflammatory cells, cytokines and chemokines in breast cancer progression: reciprocal tumor-microenvironment interactions Breast cancer stem cells, cytokine networks, and the tumor microenvironment Potential therapeutic implications of IL-6/IL-6R/gp130-targeting agents in breast cancer Proinflammatory cytokines in breast cancer: mechanisms of action and potential targets for therapeutics Cytokine network: new targeted therapy for pancreatic cancer. Current pharmaceutical design Cytokine therapy for cancer Transgenic mouse models in cancer research Mechanism-based cancer prevention approaches: targets, examples, and the use of transgenic mice Combination immunotherapy of primary prostate cancer in a transgenic mouse model using CTLA-4 blockade The interaction between HMGB1 and TLR4 dictates the outcome of anticancer chemotherapy and radiotherapy 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 The significant role of interleukin-6 and its signaling pathway in the immunopathogenesis and treatment of breast cancer Interleukin-6 signaling pathway in targeted therapy for cancer Transcriptomic dynamics of breast cancer progression in the MMTV-PyMT mouse model Constrained least-squares optimization for robust estimation of center of rotation LSimpute: accurate estimation of missing values in microarray data with least squares methods Sparse non-negative matrix factorizations via alternating non-negativity-constrained least squares for microarray data analysis Adaptive Immunity: T Cells and Cytokines Dendritic cells, inflammation and breast cancer CD133 mRNA transfected dendritic cells induces a coordinated cytotoxic and helper T cell responses against breast cancer stem cells. Molecular Therapy-Oncolytics Release of high mobility group box 1 by dendritic cells controls T cell activation via the receptor for advanced glycation end products The potential effect and mechanism of highmobility group box 1 protein on regulatory T cell-mediated immunosuppression Functional Modulation of Regulatory T Cells by IL-2 Regulatory T cells and breast cancer: implications for immunopathogenesis Regulatory T Cells and Cancer: A Two-Sided Story A survey on the role of interleukin-10 in breast cancer: A narrative Dendritic Cells and CD8 T Cell Immunity in Tumor Microenvironment Tumor-associated macrophages: a molecular perspective Increased drug resistance in breast cancer by tumorassociated macrophages through IL-10/STAT3/bcl-2 signaling pathway The role of tumor-associated macrophages in breast cancer progression (review) Macrophage Polarization: Anti-Cancer Strategies to Target Tumor-Associated Macrophage in Breast Cancer Tumor-Associated Macrophages as Major Players in the Tumor Microenvironment Cancer-associated adipocytes: key players in breast cancer progression Stroma Cells in Tumor Microenvironment and Breast Cancer IL-6 involvement in epithelial cancers. The Journal of clinical investigation The Effects of Adipocytes on the Regulation of Breast Cancer in the Tumor Microenvironment: An Update Local adipocytes enable estrogen-dependent breast cancer growth: Role of leptin and aromatase The dual role of TGFβ in human cancer: from tumor suppression to cancer metastasis. International Scholarly Research Notices CD4+ and CD8+ T cells have opposing roles in breast cancer progression and outcome Cancer-associated adipocytes exhibit an activated phenotype and contribute to breast cancer invasion. Cancer research Necrosis correlates with high vascular density and focal macrophage infiltration in invasive carcinoma of the breast HMGB1 as a therapeutic target in disease HMGB1: The Central Cytokine for All Lymphoid Cells High-mobility group box 1 (HMGB1) protein at the crossroads between innate and adaptive immunity HMGB1 in inflammation and cancer Development of hybrid small molecules that induce degradation of estrogen receptor-alpha and necrotic cell death in breast cancer cells. Cancer science Monocytic cells hyperacetylate chromatin protein HMGB1 to redirect it towards secretion Hydrogen peroxide stimulates macrophages and monocytes to actively release HMGB1 HMG-1 as a late mediator of endotoxin lethality in mice 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 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 Monocytes promote natural killer cell interferon gamma production in response to the endogenous danger signal HMGB1 Lactic acid bacterium potently induces the production of interleukin-12 and interferon-γ by mouse splenocytes Immune Tumor Microenvironment in Breast Cancer and the Participation of Estrogen and Its Receptors in Cancer Physiopathology Double roles of macrophages in human neuroimmune diseases and their animal models Role of cytokines Characteristics of intestinal dendritic cells in inflammatory bowel diseases 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 Monocyte-derived dendritic cells perform hemophagocytosis to fine-tune excessive immune responses IL-10: the master regulator of immunity to infection Interleukin-10 production by effector T cells: Th1 cells show self control Markers of dengue severity: a systematic review of cytokines and chemokines Malignant cells fuel tumor growth by educating infiltrating leukocytes to produce the mitogen Gas6 Adipocyte biology in breast cancer: From silent bystander to active facilitator Tumor necrosis factor-alpha expression in ischemic neurons Cancer associated fibroblasts express pro-inflammatory factors in human breast and ovarian tumors Determining cell type abundance and expression from bulk tissues with digital cytometry Transformed epithelial cells and fibroblasts/myofibroblasts interaction in breast tumor: a mathematical model and experiments Missing value estimation for DNA microarray gene expression data: local least squares imputation Partial least-squares discriminant analysis optimized by particle swarm optimization: application to 1H nuclear magnetic resonance analysis of lung cancer metabonomics Sensitivity analysis in chemical kinetics. Annual review of physical chemistry Numerical integration using sparse grids. Numerical algorithms Likelihood approximation by numerical integration on sparse grids Role of the tumor microenvironment in breast cancer Schrö der CP. Tumor-associated macrophages in breast cancer: Innocent bystander or important player? Cancer treatment reviews Macrophage IL-10 blocks CD8+ T cell-dependent responses to chemotherapy by suppressing IL-12 expression in intratumoral dendritic cells Macrophages impede CD8 T cells from reaching tumor cells and limit the efficacy of anti-PD-1 treatment 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 Annual report to the nation on the status of cancer, 1975-2008, featuring cancers associated with excess weight and lack of sufficient physical activity Targeting obesity-related adipose tissue dysfunction to prevent cancer development and progression Obesity hormone leptin: a new target in breast cancer? Leptin and cancer Multifaceted roles of interleukin-6 in adipocyte-breast cancer cell interaction Concise review: breast cancer stem cells: regulatory networks, stem cell niches, and disease relevance. Stem cells translational medicine Slug and Sox9 cooperatively determine the mammary stem cell state Interactions between cancer stem cells and their niche govern metastatic colonization Targeting exosomes from preadipocytes inhibits preadipocyte to cancer stem cell signaling in early-stage breast cancer. Breast cancer research and treatment Targeting AMP-activated protein kinase in adipocytes to modulate obesity-related adipokine production associated with insulin resistance and breast cancer cell proliferation Obesity promotes resistance to anti-VEGF therapy in breast cancer by up-regulating IL-6 and potentially FGF-2 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.