# Oxygen, Angiogenesis, Cancer and Immune Interplay in Breast Tumor Micro-
Environment: A Computational Investigation

Navid Mohammad Mirzaei
Department of Epidemiology
Mailman School of Public Health
Columbia University
New York, NY 10032, USA
<EMAIL_ADDRESS>
&Panayotis G. Kevrekidis
Department of Mathematics and Statistics
University of Massachusetts
Amherst, MA 01003-4515, USA
<EMAIL_ADDRESS>
&Leili Shahriyari
Department of Mathematics and Statistics
University of Massachusetts
Amherst, MA 01003-4515, USA
<EMAIL_ADDRESS>Corresponding author.Corresponding author.

###### Abstract

Breast cancer is one of the most challenging global health problems among
women. This study investigates the intricate breast tumor microenvironment
(TME) dynamics utilizing data from Mammary-specific Polyomavirus Middle T
Antigen Overexpression mouse models (MMTV-PyMT). It incorporates Endothelial
Cells (ECs), oxygen, and Vascular Endothelial Growth Factors (VEGF) to examine
the interplay of angiogenesis, hypoxia, VEGF, and the immune cells in cancer
progression. We introduce an approach to impute the immune cell fractions
within the TME using single-cell RNA-sequencing (scRNA-seq) data from MMTV-
PyMT mice. We further quantify our analysis by estimating cell counts using
cell size data and laboratory findings from existing literature. Parameter
estimation is carried out via a Hybrid Genetic Algorithm (HGA). Our
simulations reveal various TME behaviors, emphasizing the critical role of
adipocytes, angiogenesis, hypoxia, and oxygen transport in driving immune
responses and cancer progression. The global sensitivity analyses highlight
potential therapeutic intervention points, such as VEGFs’ critical role in EC
growth and oxygen transportation and severe hypoxia’s effect on the cancer and
the total number of cells. The VEGF-mediated production rate of ECs shows an
essential time-dependent impact, highlighting the importance of early
intervention in slowing cancer progression. These findings align with the
observations from the clinical trials demonstrating the efficacy of VEGF
inhibitors and suggest a timely intervention for better outcomes.

_Keywords_ Breast cancer $\cdot$ Tumor Microenvironment $\cdot$ Angiogenesis

##  1 Introduction

Breast cancer is the most prevalent cancer among women globally [58], with 2.3
million new cases and 685,000 deaths reported in 2022 [5]. The estimated cases
for the United States in 2023 are 298,000, with 43,000 deaths [111]. Surgical
removal followed by radiotherapy is a standard treatment for non-metastatic
breast cancer, varying based on subtypes such as HR+, ERBB2+, and triple-
negative [108]. For example, endocrine therapy is prescribed for HR+ subtypes
and trastuzumab-based ERBB2-directed antibody therapy in addition to
chemotherapy for ERBB2+ subtypes and chemotherapy alone for triple-negative
subtypes. Understanding the complex cellular and molecular interplays of the
Tumor Micro-Environment (TME) is crucial in achieving more effective treatment
strategies [113].

The TME, an ecosystem including tumor cells, immune cells, fibroblasts, blood
vessels, cytokines, and the extracellular matrix, is an active promoter of
cancer progression [124]. New technologies are developed to study the
intricacies within the TME, such as spatial transcriptomics and multiplexed
proteomics [36]. Longitudinal in-vivo investigations of the TME are costly and
straining, so an alternative is using genetically engineered mouse models
[109, 131, 59].

The mammary-specific polyomavirus middle T antigen overexpression mouse model
(MMTV-PyMT) is one of the most popular mouse models since its discovery in
1992 [7, 11, 48, 119], developing spontaneous mammary tumors similar to human
breast cancers [24]. The histopathology and breast cancer biomarkers
expression in the MMTV-PyMT tumors is in line with late-stage human breast
cancers [1].

Despite mouse models’ convenience, making predictions about cells and cytokine
dynamics in the TME is still challenging. However, mathematical models can
help scientists make valuable estimations, predictions, or hypotheses based on
limited data. Some mathematical cancer models focus on a certain cell type
mutation [75, 74], some investigate a set of interactions in the TME [86, 71,
83], others aim to discover optimal treatment dosage [91, 12, 64], and some
describe it as an epidemiological problem and search for the contributing risk
factors [10, 66].

Among the crucial mechanisms in TME benefiting from mathematical modeling are
those leading to vessel formation (angiogenesis) and oxygen delivery. Chaplain
has introduced a generic tumor angiogenesis modelbased on the ECs, angiogenic
cytokine, and fibronectin interactions [18]; see also the earlier impactful
work of [3]. Harrington et al. designed a mathematical model to simulate
vessel formation in the cornea in the presence of promoter and inhibitor
cytokines [51], following the earlier work of [120]. Both these and other
related (see, e.g., the review of [99]) studies focus on the Vessels’ spatial
distribution and obtain valuable results agreeing with relevant biological
observations. Nevertheless, considerable space exists for further
developments, notably the angiogenic factors (so-called TAFs) and TME
interactions.

Angiogenesis, a critical process in tumor progression, plays a central role in
breast cancer development. Tumors require a constant blood supply to sustain
their growth, and angiogenesis ensures new blood vessel formation. The solid
TME often experiences hypoxia, a condition of insufficient oxygen levels.
Hypoxia triggers the VEGF upregulation, a key pro-angiogenic factor. VEGF
promotes EC proliferation and migration, leading to new blood vessel
formation. The interactions in the TME can play a significant role in
promoting or inhibiting angiogenesis. The intricate interplays between the
TME, angiogenesis, hypoxia, and VEGF underscore the significance of targeting
these pathways for therapeutic interventions. Numerous studies have
investigated these mechanisms, with seminal works including Folkman [41] and
Ferrara et al. [39].

In this paper, we propose an ODE model describing the interactions in the TME,
an extension of a previous study [85]. The extension concerns angiogenesis and
oxygen delivery mechanisms, i.e., we include three major players: ECs, oxygen,
and VEGF. While spatial interactions among these elements are not considered
for simplicity, we acknowledge their significance as a potential future
research direction. Here, we present a methodology to prepare scRNA-seq data
for an ODE system describing TME interactions. Using HGA, we estimate
parameters and discuss simulation outcomes, augmented by sensitivity analyses
to uncover influential mechanisms on the system, cancer cells, and total
immune cells. Our findings highlight the significance of pathways related to
angiogenesis, oxygenation, and hypoxia.

###  1.1 The model

Figure 1: Interaction network. The portion inside the dashed shape is an
addition to the network used by Mohammad Mirzaei et al. [85, 87]. Despite the
importance of oxygen in the network, no arrows are coming out of it. This is
because lack of oxygen (hypoxia) is the key factor, not oxygen itself.

The model in this paper extends the framework established by Mohammad Mirzaei
et al. [85]. In our adaptation, we have incorporated additional state
variables: VEGF, ECs, and Oxygen. These new elements are added to investigate
hypoxia and angiogenesis roles in tumor progression. For a detailed
representation of the interaction network, see Figure 1. The newly
incorporated interactions are distinctly marked for easy identification. Table
1 shows the model variables with their corresponding equations.

The naming conventions for our parameters are as follows:

$\displaystyle\lambda_{XY}$ The rate molecule/cell Y promotes the production
or proliferation of molecule/cell X $\displaystyle\delta_{XY}$ The rate
molecule/cell Y inhibits the production or proliferation of molecule/cell X
$\displaystyle\delta_{XYOx}$ The rate molecule/cell Y inhibits the production
of molecule/cell X proportional to the hypoxia function value
$\displaystyle\delta_{X}$ Natural decay/death rate of molecule/cell X
$\displaystyle A_{X}$ Innate production rate of the cell X
$\displaystyle\delta_{i}$ $\displaystyle\quad\text{For }i=1,\cdots,6.\text{
Oxygen consumption rate by}$ the six major cell types in the model

In this paper, parameters with overhead bars are in non-dimensionalized form.
Mathematically, we consider production processes, both for a single species or
stimulated by the presence of a second species. Additionally, we account for
degradation processes, either occurring naturally for a single species or
stimulated by the presence of or interaction with another species.

We used RNA-sequencing data from three MMTV-PyMT mice [14]. The proportions of
various immune cells were estimated using the digital cytometry technique
CIBERSORTx [92], with the signature matrix extracted from the Immunological
Genome Project [52]. We created a unique signature matrix for non-immune cell
populations using single-cell RNA sequencing data from the Tabula Muris
database [106]. Figure 2 show imputed fractions for immune and non-immune
cells. Additionally, we applied a laboratory protocol tailored for MMTV-PyMT
mice to quantify macrophages and cancer cells relative to tumor size. This
process approximated the total cell count, informed by imputed cell fractions,
and allowed the direct extraction of cytokine levels from the RNA-sequencing
data.

Variable | Name | Equation for their rate of changes
---|---|---
$E_{N}$ | ECs | $\begin{array}[]{l}\\\ \frac{d[\overline{E_{N}}]}{dt}=\overline{\lambda}_{VE_{N}}[\overline{V}]\Big{(}1-\frac{[\overline{E_{N}}]}{\overline{E_{N}}_{0}}\Big{)}[\overline{E_{N}]}-\overline{\delta}_{E_{N}IL_{12}}[\overline{IL_{12}}][\overline{E_{N}}]\\\ -\Big{(}{\delta}_{E_{N}}+\overline{\delta}_{E_{N}O_{x}}f_{2}([\overline{O_{x}}])\Big{)}[\overline{E_{N}}]\vspace{0.1in}\end{array}$
$V$ | VEGF | $\begin{array}[]{l}\\\ {{\frac{d[\overline{V}]}{dt}}}={{\overline{\lambda}_{VA}}[\overline{A}]+{\overline{\lambda}_{VM}}[\overline{M}]+{{\overline{\lambda}_{VCO_{x}}[\overline{C}]f_{1}(O_{x})}}+{{\overline{\lambda}_{VO_{x}}f_{1}([\overline{O_{x}}])}}}-{\delta_{V}[\overline{V}]\vspace{0.1in}}\end{array}$
$O_{x}$ | Oxygen | $\begin{array}[]{l}\\\ {{\frac{d[\overline{O_{x}}]}{d\tau}}}=\overline{\lambda}_{O_{x}E_{N}}([\overline{E_{N}}])-\overline{\delta}_{1}\left([\overline{T_{N}}]+[\overline{T_{h}}]+[\overline{T_{c}}]+[\overline{T_{r}}]\right)[\overline{O_{x}}]-\overline{\delta}_{2}\left([\overline{D_{N}}]+[\overline{D}]\right)[\overline{O_{x}}]\\\ -\overline{\delta}_{3}\left([\overline{M_{N}}]+[\overline{M}]\right)[\overline{O_{x}}]-\overline{\delta}_{4}[\overline{C}][\overline{O_{x}}]-\overline{\delta}_{5}[\overline{A}][\overline{O_{x}}]-\overline{\delta}_{6}[\overline{E_{N}}][\overline{O_{x}}]-\delta_{O_{x}}[\overline{O_{x}}]\vspace{0.1in}\end{array}$
$T_{N}$ | Naive T-cells | $\begin{array}[]{l}\\\ \frac{d[\overline{T_{N}}]}{dt}=\overline{A}_{T_{N}}-\left(\overline{\lambda}_{T_{h}H}[\overline{H}]+\overline{\lambda}_{T_{h}D}[\overline{D}]+\overline{\lambda}_{T_{h}IL_{12}}[\overline{IL_{12}}]\right)[\overline{T_{N}}]-\left(\overline{\lambda}_{T_{c}D}[\overline{D}]+\overline{\lambda}_{T_{c}IL_{12}}[\overline{IL_{12}}]\right)[\overline{T_{N}}]\\\ -\left(\overline{\lambda}_{T_{r}D}[\overline{D}]+\overline{\lambda}_{T_{r}Ox}f_{1}([\overline{Ox}])\right)[\overline{T_{N}}]-\left(\delta_{T_{N}}+\overline{\delta}_{T_{N}Ox}f_{2}([\overline{Ox}])\right)[\overline{T_{N}}]\vspace{0.1in}\end{array}$
$T_{h}$ | Helper T-cells | $\begin{array}[]{l}\\\ \frac{d[\overline{T_{h}}]}{dt}=\left(\overline{\lambda}_{T_{h}H}[\overline{H}]+\overline{\lambda}_{T_{h}D}[\overline{D}]+\overline{\lambda}_{T_{h}IL_{12}}[\overline{IL_{12}}]\right)[\overline{T_{N}}]-\left(\overline{\delta}_{T_{h}T_{r}}[\overline{T_{r}}]+\overline{\delta}_{T_{h}IL_{10}}[\overline{IL_{10}}]\right)[\overline{T_{h}}]\\\ -\left({{\delta}_{T_{h}}}+{{\overline{\delta}_{T_{h}O_{x}}f_{2}([\overline{O_{x}}])}}\right)[\overline{T_{h}}]\vspace{0.1in}\end{array}$
$T_{C}$ | Cytotoxic cells | $\begin{array}[]{l}\\\ \frac{d[\overline{T_{c}}]}{dt}=\left(\overline{\lambda}_{T_{c}D}[\overline{D}]+\overline{\lambda}_{T_{c}IL_{12}}[\overline{IL_{12}}]\right)[\overline{T_{N}}]-\left(\overline{\delta}_{T_{c}T_{r}}[\overline{T_{r}}]+\overline{\delta}_{T_{c}IL_{10}}[\overline{IL_{10}}]\right)[\overline{T_{c}}]\\\ -\left({{\delta}_{T_{c}}}+{{\overline{\delta}_{T_{c}O_{x}}f_{2}([\overline{O_{x}}])}}\right)[\overline{T_{c}}]\vspace{0.1in}\end{array}$
$T_{r}$ | Regulatory T-cells | $\begin{array}[]{l}\\\ \frac{d[\overline{T_{r}}]}{dt}=\left(\overline{\lambda}_{T_{r}D}[\overline{D}]+\overline{\lambda}_{T_{r}O_{x}}\ f_{1}([\overline{O_{x}}])\right)\ [\overline{T_{N}}]-\left({{\delta}_{T_{r}}}+{{\overline{\delta}_{T_{r}O_{x}}f_{2}([\overline{O_{x}}])}}\right)[\overline{T_{r}}]\vspace{0.1in}\end{array}$
$D_{N}$ | Naive DCs | $\begin{array}[]{l}\\\ \frac{d[\overline{D_{N}}]}{dt}=\overline{A}_{D_{N}}-\left(\overline{\lambda}_{DC}[\overline{C}]+\overline{\lambda}_{DH}[\overline{H}]+{{\overline{\lambda}_{DO_{x}}f_{1}([\overline{O_{x}}])}}\right)[\overline{D_{N}}]-\left({{\delta}_{D_{N}}}+{{\overline{\delta}_{D_{N}O_{x}}f_{2}([\overline{O_{x}}])}}\right)[\overline{D_{N}}]\vspace{0.1in}\end{array}$
$D$ | Activated DCs | $\begin{array}[]{l}\\\ \frac{d[\overline{D}]}{dt}=\left(\overline{\lambda}_{DC}[\overline{C}]+\overline{\lambda}_{DH}[\overline{H}]+{{\overline{\lambda}_{DO_{x}}f_{1}([\overline{O_{x}}])}}\right)[\overline{D_{N}}]-\left(\overline{\delta}_{DC}[\overline{C}]+{{\delta}_{D}}+{{\overline{\delta}_{D_{N}O_{x}}f_{2}([\overline{O_{x}}])}}\right)[\overline{D}]\vspace{0.1in}\end{array}$
$M_{N}$ | Naive Macrophages | $\begin{array}[]{l}\\\ \frac{d[\overline{M_{N}}]}{dt}=\overline{A}_{M_{N}}-\Big{(}\overline{\lambda}_{MIL_{10}}[\overline{IL_{10}}]+\overline{\lambda}_{MIL_{12}}[\overline{IL_{12}}]+\overline{\lambda}_{MT_{h}}[\overline{T_{h}}]+{\delta}_{M_{N}}\Big{)}[\overline{M_{N}}]\vspace{0.1in}\end{array}$
$M$ | Macrophages | $\begin{array}[]{l}\\\ \frac{d[\overline{M}]}{dt}=\left(\overline{\lambda}_{MIL_{10}}[\overline{IL_{10}}]+\overline{\lambda}_{MIL_{12}}[\overline{IL_{12}}]+\overline{\lambda}_{MT_{h}}[\overline{T_{h}}]\right)[\overline{M_{N}}]-{{\delta}_{M}}[\overline{M}]\vspace{0.1in}\end{array}$
$C$ | Cancer cells | $\begin{array}[]{l}\\\ \frac{d[\overline{C}]}{dt}=\left({\lambda}_{C}+\overline{\lambda}_{CIL_{6}}[\overline{IL_{6}}]+\overline{\lambda}_{CA}[\overline{A}]\right)\left(1-\frac{[\overline{C}]}{\overline{C_{0}}}\right)[\overline{C}]-\left(\overline{\delta}_{CT_{c}}[\overline{T_{c}}]+{{\overline{\delta}_{CT_{c}O_{x}}f_{1}([\overline{O_{x}}])}}[\overline{T_{c}}]\right)[\overline{C}]\\\ -\left({{\overline{\delta}_{CO_{x}}f_{2}([\overline{O_{x}}])}}+\delta_{C}\right)[\overline{C}]\vspace{0.1in}\end{array}$
$N$ | Necrotic cells | $\begin{array}[]{l}\\\ \frac{d[\overline{N}]}{dt}=\overline{\alpha}_{NC}\left(\overline{\delta}_{CT_{c}}[\overline{T_{c}}]+{{\overline{\delta}_{CT_{c}O_{x}}f_{1}([\overline{O_{x}}])}}[\overline{T_{c}}]+{{\overline{\delta}_{CO_{x}}f_{2}([\overline{O_{x}}])}}+\delta_{C}\right)[\overline{C}]-{\delta}_{N}[\overline{N}]\vspace{0.1in}\end{array}$
$A$ | Adipocytes | $\begin{array}[]{l}\\\ \frac{d[\overline{A}]}{dt}={\lambda}_{A}\left(1-\frac{[\overline{A}]}{\overline{A_{0}}}\right)[\overline{A}]-\delta_{A}[\overline{A}]\vspace{0.1in}\end{array}$
$H$ | HMGB1 | $\begin{array}[]{l}\\\ \frac{d[\overline{H}]}{dt}=\overline{\lambda}_{HD}[\overline{D}]+\overline{\lambda}_{HN}[\overline{N}]+\overline{\lambda}_{HM}[\overline{M}]+\overline{\lambda}_{HT_{c}}[\overline{T_{c}}]+\overline{\lambda}_{HC}[\overline{C}]-{\delta}_{H}[\overline{H}]\vspace{0.1in}\end{array}$
$IL_{12}$ | IL-12 | $\begin{array}[]{l}\\\ \frac{d[\overline{IL_{12}}]}{dt}=\overline{\lambda}_{IL_{12}M}[\overline{M}]+\overline{\lambda}_{IL_{12}D}[\overline{D}]-{\delta}_{IL_{12}}[\overline{IL_{12}}]\vspace{0.1in}\end{array}$
$IL_{10}$ | IL-10 | $\begin{array}[]{l}\\\ \frac{d[\overline{IL_{10}}]}{dt}=\overline{\lambda}_{IL_{10}M}[\overline{M}]+\overline{\lambda}_{IL_{10}D}[\overline{D}]+\overline{\lambda}_{IL_{10}T_{r}}[\overline{T_{r}}]+\overline{\lambda}_{IL_{10}T_{h}}[\overline{T_{h}}]+\overline{\lambda}_{IL_{10}T_{c}}[\overline{T_{c}}]+\overline{\lambda}_{IL_{10}C}[\overline{C}]\\\ +{{\overline{\lambda}_{IL_{10}O_{x}}f_{1}([\overline{O_{x}}])}}-{\delta}_{IL_{10}}[\overline{IL_{10}}]\vspace{0.1in}\end{array}$
$IL_{6}$ | IL-6 | $\begin{array}[]{l}\\\ \frac{d[\overline{IL_{6}}]}{dt}=\overline{\lambda}_{IL_{6}A}[\overline{A}]+\overline{\lambda}_{IL_{6}M}[\overline{M}]+{{\overline{\lambda}_{IL_{6}MO_{x}}[\overline{M}]f_{1}([\overline{O_{x}}])}}+\overline{\lambda}_{IL_{6}D}[\overline{D}]+{{\overline{\lambda}_{IL_{6}O_{x}}f_{1}([\overline{O_{x}}])}}\\\ -\delta_{IL_{6}}[\overline{IL_{6}}]\vspace{0.1in}\end{array}$

Table 1:  Model variables with their equations.

We employed a Hybrid Genetic Algorithm (HGA) for parameter estimation in our
model, followed by a global sensitivity analysis to determine parameters
significantly affecting overall system dynamics, cancer and total immune cell
counts. Perturbing these sensitive parameters allows us to assess and quantify
their impact on model outcomes, as detailed in the Methods section.

Figure 2: Cell fractions for (A)Immune cells and (B) Non-immune cells. The
proportions are calculated for three mice at week 6 (hyperplasia), week 8
(adenoma/MIN), week 10 (early carcinoma), and week 12 (late carcinoma).

##  2 Methods

### Model variables and their interactions

We use the mass action law to create the ODE system for the interactions among
variables as discussed in [84].

#### Hypoxia

Hypoxic effects take hours to impact the system [69], while the simulation
time for our problem and the associated data is on the order of weeks.
Therefore, we assume low oxygen levels entail hypoxic effects on cells for
each computational time step corresponding to one day. Usually, laboratory
cell cultures are in a hyperoxic condition with about 20.9% oxygenation [90].
However, the in-vivo normoxic oxygen level in breast tissue is about 8.5%, and
the hypoxic level is below 1% [126]. Although different T-cell types have
different hypoxic responses at different oxygen levels, most of these hypoxic
activities occur in the range 0%-1% [34]. In particular, at 0.1%-1%, there is
not enough oxygen in the environment to disintegrate HIF-1$\alpha$, so the
latter will accumulate and trigger some hypoxic functions and below 0.1% most
cell cycles yield, causing hypoxic death [34]. For simplicity, if a particular
cell or cytokine is promoted or inhibited via HIF pathways, which are not
included in our model, we use a rate independent of other state variables and
only dependent on oxygen levels. Some hypoxia outcomes through HIF pathways
included in our model are up-regulating the VEGF, IL-6, and IL-10 production
[94], promoting the maturation of the dendritic cells [68, 40], increasing the
regulatory T-cells differentiation [19, 26] and increasing the T-cells
cytotoxicity [37, 45]. These outcomes indirectly lead to reduced cytotoxic and
helper T-cell differentiation. Also, hypoxia promotes the macrophages’
differentiation into a pro-tumor subtype [27]. However, since we do not
include the macrophage subtypes in our model, we include this effect through
the IL-6 pathway. In other words, we assume hypoxia promotes IL-6 production
by macrophages, which has a strong pro-inflammatory effect.

In this paper, we assume blood vessel formation is proportional to the number
of endothelial cells (ECs), facilitating oxygen transportation. Thus, oxygen
serves as a source dependent on ECs. Regions with oxygen levels below a
critical threshold are considered hypoxic. We specify two hypoxic states: mild
hypoxia initiating the HIF pathways (between 0.1%-1% oxygen) and severe
hypoxia causing cell cycle arrest (below 0.1%). Cells can die through
mechanisms in the original model (via inhibitors or normoxic death) or undergo
hypoxic death. To capture the promotion or inhibition effects of hypoxia, we
define the following function:

$f_{i}(x)=\frac{(Ox_{crit}^{(i)})^{m}}{(Ox_{crit}^{(i)})^{m}+x^{m}},\qquad
i=1,2$ (1)

where $Ox_{crit}^{(i)}$ are the critical oxygen values and $m$ is the Hill-
like coefficient, which governs the decrease rate from the maximum to the
minimum function value. Notice that two functions are defined in (1),
corresponding to two critical oxygen values. The first critical value
corresponds to the transitioning threshold from normoxia to mild hypoxia, and
the second one marks the transition from mild to severe hypoxia.

#### T-cells

We classify T-cells into naive ($T_{N}$), helper ($T_{h}$), cytotoxic
($T_{C}$), and regulatory ($T_{r}$) types. Activated subtypes differentiate
from the naive subtype, introducing negative feedback terms proportional to
their level and corresponding activator levels. These terms serve as positive
feedback for the activated subtypes, indicating regulation by the available
number of naive subtypes. We employ a similar method for naive and activated
dendritic cells and macrophages.

We consider a constant intrinsic proliferation rate for naive T-cells. Naive
T-cells diminish due to differentiation or normoxic and hypoxic cell death.

Helper T-cells’ differentiation from naive T-cells is promoted by HMGB1 [32,
138], dendritic cells [115, 49], and IL-12 [60, 25]. They are inhibited by
regulatory T-cells [9, 135] and IL10 [104, 110]. We also consider normoxic and
hypoxic cell death for helper T-cells.

Similarly, cytotoxic T-cells are promoted by dendritic cells [43, 65] and
IL-12 [25, 82]. They are inhibited by regulatory T-cells [133, 9] and IL-10
[82, 105]. We consider normoxic and hypoxic cell death for these T-cells, too.

Like the other types, regulatory T-cell differentiation is promoted by
dendritic cells [98, 121]. Also, is mentioned in the “Hypoxia" subsection,
mild hypoxia can increase regulatory T-cell differentiation. Similarly to the
case of the other subtypes, we close the ODE by the normoxic and hypoxic cell
death.

Note, the hypoxic cell death is modeled via $f_{2}$, but the hypoxic promotion
of $T_{r}$ is modeled via $f_{1}$, corresponding to severe and mild hypoxia,
respectively. The corresponding ODEs in this case read:

$\displaystyle\frac{d[\overline{T_{N}}]}{dt}$
$\displaystyle=\overline{A}_{T_{N}}$
$\displaystyle-\left(\overline{\lambda}_{T_{h}H}[\overline{H}]+\overline{\lambda}_{T_{h}D}[\overline{D}]+\overline{\lambda}_{T_{h}IL_{12}}[\overline{IL_{12}}]\right)[\overline{T_{N}}]$
$\displaystyle-\left(\overline{\lambda}_{T_{c}D}[\overline{D}]+\overline{\lambda}_{T_{c}IL_{12}}[\overline{IL_{12}}]\right)[\overline{T_{N}}]$
$\displaystyle-\left(\overline{\lambda}_{T_{r}D}[\overline{D}]+\overline{\lambda}_{T_{r}Ox}f_{1}([\overline{Ox}])\right)[\overline{T_{N}}]$
$\displaystyle-\left(\delta_{T_{N}}+\overline{\delta}_{T_{N}Ox}f_{2}([\overline{Ox}])\right)[\overline{T_{N}}]$
(2) $\displaystyle\frac{d[\overline{T_{h}}]}{dt}$
$\displaystyle=X\left(\overline{\lambda}_{T_{h}H}[\overline{H}]+\overline{\lambda}_{T_{h}D}[\overline{D}]+\overline{\lambda}_{T_{h}IL_{12}}[\overline{IL_{12}}]\right)[\overline{T_{N}}]$
$\displaystyle-\left(\overline{\delta}_{T_{h}T_{r}}[\overline{T_{r}}]+\overline{\delta}_{T_{h}IL_{10}}[\overline{IL_{10}}]\right)[\overline{T_{h}}]$
$\displaystyle-\left({{\delta}_{T_{h}}}+{{\overline{\delta}_{T_{h}O_{x}}f_{2}([\overline{O_{x}}])}}\right)[\overline{T_{h}}],$
(3) $\displaystyle\frac{d[\overline{T_{c}}]}{dt}$
$\displaystyle=\left(\overline{\lambda}_{T_{c}D}[\overline{D}]+\overline{\lambda}_{T_{c}IL_{12}}[\overline{IL_{12}}]\right)[\overline{T_{N}}]$
$\displaystyle-\left(\overline{\delta}_{T_{c}T_{r}}[\overline{T_{r}}]+\overline{\delta}_{T_{c}IL_{10}}[\overline{IL_{10}}]\right)[\overline{T_{c}}]$
$\displaystyle-\left({{\delta}_{T_{c}}}+{{\overline{\delta}_{T_{c}O_{x}}f_{2}([\overline{O_{x}}])}}\right)[\overline{T_{c}}],$
(4) $\displaystyle\frac{d[\overline{T_{r}}]}{dt}$
$\displaystyle=\left(\overline{\lambda}_{T_{r}D}[\overline{D}]+\overline{\lambda}_{T_{r}O_{x}}\
f_{1}([\overline{O_{x}}])\right)\ [\overline{T_{N}}]$
$\displaystyle-\left({{\delta}_{T_{r}}}+{{\overline{\delta}_{T_{r}O_{x}}f_{2}([\overline{O_{x}}])}}\right)[\overline{T_{r}}],$
(5)

#### Dendritic cells

We model two dendritic cell types: naive ($D_{N}$) and activated ($D$)
dendritic cells. Like T-cells, the activated dendritic cells differentiate
from the naive ones. We also include normoxic and hypoxic cell death for both
subtypes. Like naive T-cells, we consider a constant intrinsic growth rate for
the naive dendritic cells. Their number decreases as they get activated or
die. The tumor burden produces danger signals that promote dendritic cell
maturation [76]. HMGB1 is another factor contributing to dendritic activation
[123, 15] and so is hypoxia, as explained in the "Hypoxia" subsection. Despite
increasing the dendritic cell maturation, cancer cells can release factors
inhibiting dendritic cell maturation and function [121].

$\displaystyle\frac{d[\overline{D_{N}}]}{dt}$
$\displaystyle=\overline{A}_{D_{N}}$
$\displaystyle-\left(\overline{\lambda}_{DC}[\overline{C}]+\overline{\lambda}_{DH}[\overline{H}]+{{\overline{\lambda}_{DO_{x}}f_{1}([\overline{O_{x}}])}}\right)[\overline{D_{N}}]$
$\displaystyle-\left({{\delta}_{D_{N}}}+{{\overline{\delta}_{D_{N}O_{x}}f_{2}([\overline{O_{x}}])}}\right)[\overline{D_{N}}],$
(6) $\displaystyle\frac{d[\overline{D}]}{dt}$
$\displaystyle=\left(\overline{\lambda}_{DC}[\overline{C}]+\overline{\lambda}_{DH}[\overline{H}]+{{\overline{\lambda}_{DO_{x}}f_{1}([\overline{O_{x}}])}}\right)[\overline{D_{N}}]$
$\displaystyle-\left(\overline{\delta}_{DC}[\overline{C}]+{{\delta}_{D}}+{{\overline{\delta}_{D_{N}O_{x}}f_{2}([\overline{O_{x}}])}}\right)[\overline{D}],$
(7)

#### Macrophages

We consider naive ($M_{N}$) and activated macrophages ($M$) in the model. For
simplicity, the different activated macrophages, such as M1 and M2 types, are
consolidated as $M$. Differentiation of these types along the lines (e.g., see
[13]) would be an interesting topic for further study. Furthermore, we do not
include hypoxic cell death for macrophages. Since these cell types tend to
accumulate in hypoxic regions, they quickly adapt to severe hypoxia [50]. Like
the other naive cell types, we model naive macrophages to proliferate at an
innate constant rate. Their activation is augmented by IL-10 [81, 89], IL-12
[82, 114] and helper T-cells [4]. Both macrophage subtypes die via normoxic
cell death. Activated macrophage polarity is defined through its interaction
with cancer cells, which is why they are often referred to as anti (M1) and
pro (M2) tumor macrophages [22]. However, as mentioned above, we do not
explicitly include these subtypes in our model, and their different
interaction with tumor cells is modeled through their cytokine pathways.

$\displaystyle\frac{d[\overline{M_{N}}]}{dt}$
$\displaystyle=\overline{A}_{M_{N}}-\Big{(}\overline{\lambda}_{MIL_{10}}[\overline{IL_{10}}]+\overline{\lambda}_{MIL_{12}}[\overline{IL_{12}}]$
$\displaystyle+\overline{\lambda}_{MT_{h}}[\overline{T_{h}}]+{\delta}_{M_{N}}\Big{)}[\overline{M_{N}}],$
(8) $\displaystyle\frac{d[\overline{M}]}{dt}$
$\displaystyle=\Big{(}\overline{\lambda}_{MIL_{10}}[\overline{IL_{10}}]+\overline{\lambda}_{MIL_{12}}[\overline{IL_{12}}]$
$\displaystyle+\overline{\lambda}_{MT_{h}}[\overline{T_{h}}]\Big{)}[\overline{M_{N}}]-{{\delta}_{M}}[\overline{M}],$
(9)

#### Cancer and Necrotic cells

Cancer cells (C) grow according to a logistic model with carrying capacity
$C_{0}$. Many factors can promote cancer proliferation, among which we have
explicitly included the IL-6 [30, 46] and adipocyte effects [132, 44]. The
other factors are considered via a constant rate $\lambda_{C}$. Cancer cells
can be cleared via cytotoxic cells [33, 101]. Also, as mentioned, T-cell
cytotoxicity increases through HIF pathways, which is modeled by the term
including $f_{1}([\overline{O_{x}}])$. Finally, we include hypoxic and
normoxic cell death mechanisms. Even though cancer cells proliferate much
faster than normal cells, factors such as physical space, fitness landscapes,
and nutrition impose a carrying capacity on their growth [134, 103]. We apply
the same rationale to the other two cell types (described in the following
sections), namely adipocytes and ECs.

A portion of cancer cell death turns into necrotic cells [70]. That portion is
modeled via the constant $\overline{\alpha}_{NC}$. We only consider regular
cell clearance for necrotic cells since they are dead cells and do not go
through hypoxia.

$\displaystyle\frac{d[\overline{C}]}{dt}$
$\displaystyle=\left({\lambda}_{C}+\overline{\lambda}_{CIL_{6}}[\overline{IL_{6}}]+\overline{\lambda}_{CA}[\overline{A}]\right)\left(1-\frac{[\overline{C}]}{\overline{C_{0}}}\right)[\overline{C}]$
$\displaystyle-\left(\overline{\delta}_{CT_{c}}[\overline{T_{c}}]+{{\overline{\delta}_{CT_{c}O_{x}}f_{1}([\overline{O_{x}}])}}[\overline{T_{c}}]\right)[\overline{C}]$
$\displaystyle-\left({{\overline{\delta}_{CO_{x}}f_{2}([\overline{O_{x}}])}}+\delta_{C}\right)[\overline{C}],$
(10) $\displaystyle\frac{d[\overline{N}]}{dt}$
$\displaystyle=\overline{\alpha}_{NC}\Big{(}\overline{\delta}_{CT_{c}}[\overline{T_{c}}]+{{\overline{\delta}_{CT_{c}O_{x}}f_{1}([\overline{O_{x}}])}}[\overline{T_{c}}]$
$\displaystyle+{{\overline{\delta}_{CO_{x}}f_{2}([\overline{O_{x}}])}}+\delta_{C}\Big{)}[\overline{C}]-{\delta}_{N}[\overline{N}],$
(11)

#### Adipocytes

Adipocytes contribute to cancer progression through intricate mechanisms.
Having that contribution included in (10) as
$\overline{\lambda}_{CA}[\overline{A}]$, for simplicity, we model the
adipocytes dynamics through a simple logistic model.

$\displaystyle\frac{d[\overline{A}]}{dt}$
$\displaystyle={\lambda}_{A}\left(1-\frac{[\overline{A}]}{\overline{A_{0}}}\right)[\overline{A}]-\delta_{A}[\overline{A}],$
(12)

#### Cytokines

HMGB1 is produced by dendritic cells, necrotic cells, macrophages, cytotoxic
T-cells, and cancer cells [15, 57]. IL-12 is secreted by dendritic cells and
macrophages [73]. IL-10 is a product of dendritic cells, macrophages, cancer
cells, and cytotoxic, helper, and regulatory T-cells [73]. Finally,
adipocytes, dendritic cells, and macrophages are the main IL-6 producers [20,
72]. Note we have a hypoxic promotion of IL-6 production by macrophages. As
mentioned in the "Hypoxia" section, this is to model the increase in pro-tumor
macrophage differentiation via mild hypoxia. Additionally, IL-10 and IL-6
increase through HIF pathways, as explained in the "Hypoxia" section. We
include a natural decay term for cytokine clearance.

$\displaystyle\frac{d[\overline{H}]}{dt}$
$\displaystyle=\overline{\lambda}_{HD}[\overline{D}]+\overline{\lambda}_{HN}[\overline{N}]+\overline{\lambda}_{HM}[\overline{M}]+\overline{\lambda}_{HT_{c}}[\overline{T_{c}}]$
$\displaystyle+\overline{\lambda}_{HC}[\overline{C}]-{\delta}_{H}[\overline{H}],$
(13) $\displaystyle\frac{d[\overline{IL_{12}}]}{dt}$
$\displaystyle=\overline{\lambda}_{IL_{12}M}[\overline{M}]+\overline{\lambda}_{IL_{12}D}[\overline{D}]-{\delta}_{IL_{12}}[\overline{IL_{12}}],$
(14) $\displaystyle\frac{d[\overline{IL_{10}}]}{dt}$
$\displaystyle=\overline{\lambda}_{IL_{10}M}[\overline{M}]+\overline{\lambda}_{IL_{10}D}[\overline{D}]+\overline{\lambda}_{IL_{10}T_{r}}[\overline{T_{r}}]$
$\displaystyle+\overline{\lambda}_{IL_{10}T_{h}}[\overline{T_{h}}]+\overline{\lambda}_{IL_{10}T_{c}}[\overline{T_{c}}]+\overline{\lambda}_{IL_{10}C}[\overline{C}]$
$\displaystyle+{{\overline{\lambda}_{IL_{10}O_{x}}f_{1}([\overline{O_{x}}])}}-{\delta}_{IL_{10}}[\overline{IL_{10}}],$
(15) $\displaystyle\frac{d[\overline{IL_{6}}]}{dt}$
$\displaystyle=\overline{\lambda}_{IL_{6}A}[\overline{A}]+\overline{\lambda}_{IL_{6}M}[\overline{M}]+{{\overline{\lambda}_{IL_{6}MO_{x}}[\overline{M}]f_{1}([\overline{O_{x}}])}}$
$\displaystyle+\overline{\lambda}_{IL_{6}D}[\overline{D}]+{{\overline{\lambda}_{IL_{6}O_{x}}f_{1}([\overline{O_{x}}])}}-\delta_{IL_{6}}[\overline{IL_{6}}].$
(16)

#### Newly added compartments: ECs, VEGF and Oxygen

ECs proliferate on their own at a very slow rate [88]. However, with VEGFs
present, this growth rate increases significantly [130, 16]; so, we utilize a
VEGF-catalyzed logistic growth for their population. In contrast, endogenous
angiogenesis inhibitors, such as IL-12, hinder EC’s proliferation [95]. There
are other TAFs and endogenous angiogenesis inhibitors, but we only include
VEGF and IL-12 here for simplicity. We consider normoxic and hypoxic cell
death for ECs. Macrophages, Adipocytes, and hypoxic Cancer cells secrete VEGF
[16, 42, 38]. Finally, the ECs, which form vessels, will be the oxygen source,
and all the cells consume oxygen. The consumption rates are grouped by cell
type to avoid complexity.

$\displaystyle{{\frac{d[\overline{E_{N}}]}{dt}}}$
$\displaystyle={\overline{\lambda}_{VE_{N}}[\overline{V}]\Big{(}1-\frac{[\overline{E_{N}}]}{\overline{E_{N}}_{0}}\Big{)}[\overline{E_{N}}]}-\overline{\delta}_{E_{N}IL_{12}}[\overline{IL_{12}}][\overline{E_{N}}]$
$\displaystyle-\Big{(}{\delta}_{E_{N}}+\overline{\delta}_{E_{N}O_{x}}f_{2}([\overline{O_{x}}])\Big{)}[\overline{E_{N}}],$
(17) $\displaystyle{{\frac{d[\overline{V}]}{dt}}}$
$\displaystyle={\overline{\lambda}_{VA}}[\overline{A}]+{\overline{\lambda}_{VM}}[\overline{M}]+{{\overline{\lambda}_{VCO_{x}}[\overline{C}]f_{1}(Ox)}}$
$\displaystyle+{{\overline{\lambda}_{VO_{x}}f_{1}([\overline{O_{x}}])}}-{\delta_{V}[\overline{V}]}$
(18) $\displaystyle{{\frac{d[\overline{Ox}]}{d\tau}}}$
$\displaystyle=\overline{\lambda}_{O_{x}E_{N}}([\overline{E_{N}}])-\overline{\delta}_{1}\left([\overline{T_{N}}]+[\overline{T_{h}}]+[\overline{T_{c}}]+[\overline{T_{r}}]\right)[\overline{O_{x}}]$
$\displaystyle-\overline{\delta}_{2}\left([\overline{D_{N}}]+[\overline{D}]\right)[\overline{O_{x}}]-\overline{\delta}_{3}\left([\overline{M_{N}}]+[\overline{M}]\right)[\overline{O_{x}}]$
$\displaystyle-\overline{\delta}_{4}[\overline{C}][\overline{O_{x}}]-\overline{\delta}_{5}[\overline{A}][\overline{O_{x}}]-\overline{\delta}_{6}[\overline{E_{N}}][\overline{O_{x}}]-\delta_{O_{x}}[\overline{O_{x}}]$
(19)

### Data preparation

We utilized RNA-sequencing data for three mice, accessible through the Gene
Expression Omnibus (GEO) database under the accession code GSE76772 [14],
encompassing four distinct tumor progression stages: hyperplasia at week 6,
adenoma/MIN at week 8, early carcinoma at week 10, and late carcinoma at week
12. This data was primarily used to compare gene expression patterns across
various cancer stages. Raw gene expression data was obtained through
directional RNA-sequencing, and statistical techniques were employed to
eliminate genes with low transcriptional activity. Then the refined data set
was normalized through Differential Expression analysis for Sequence count
data (DESeq) [2].

We use a digital cytometry method of CIBERSORTx B-mode [92] to infer the
different cell fractions from the RNA-sequencing data. For this task, we need
a signature matrix containing information about cell-specific gene expression
signatures [63, 6]. CIBERSORTx has an in-built signature matrix for homo
sapiens, but we need one for mice. So, we use the ImmGen signature matrix from
the Immunological Genome Project [52]. It contains gene expression information
for 207 immune cell types, most are certain immune cell subtypes consolidated
as one variable in our model. For non-immune cells, a good resource is the
Tabula Muris database [106]. They have a single-cell RNA-sequence (scRNA-seq)
atlas for all important mouse organs and tissues. This data consists of two
sets: data acquired from Fluorescence-activated cell sorting (FACS) and data
acquired from Droplets. Based on the above paper, the former shows a better
variety and is compatible with Immgen and CIBERSORTx. However, unlike the
ImmGen data, this data cannot be readily used as a signature matrix. First, we
must cross-reference the scRNA-seq count data with the FACS annotations to
convert the obscure cell barcodes to meaningful cell names. Then, we upload
our polished scRNA-seq data to CIBERSORTx to create a signature matrix. This
will create a signature matrix to infer the non-immune cell fractions from the
gene expression data.

We need a framework to estimate the number of cells from the calculated cell
fractions. Sun et al. share their procedure of isolating macrophages and
cancer cells for mouse breast cancer in a detailed protocol [116]. Their
samples are derived from the same mouse model used in the study from which we
obtained our data. According to their protocol, they can extract $4.5\pm
0.9\times 10^{7}$ cancer cells from 1 gram of tumor with more than 90%
precision. This number is obtained from two 0.5-gram tumors, each almost 10
millimeters in diameter. They mention that they avoid collecting necrotic
cells. They also collected $2.2\pm 0.6\times 10^{6}$ tumor-associated
macrophages (TAMs) from 1 gram of the same tumor with the same precision. In a
study serving as the basis for this protocol, Sun et al. also provide
information on tumor volumes at various stages [117]. According to their
paper, the tumor volume measures approximately 200 $mm^{3}$ during
hyperplasia, 400 $mm^{3}$ during adenoma, 700 $mm^{3}$ during early carcinoma,
and 1400 $mm^{3}$ during late carcinoma. These measurements are calculated in
caliper format $L\times W^{2}/2$, where L is the long diameter and W is the
short diameter. Assuming the tumor sections are close to spheres, we can take
L to be approximately equal to W. Hence, during hyperplasia, the tumor is 7.36
mm in diameter, 9.28 mm during adenoma, 11.18 mm during early carcinoma, and
14.09 mm during late carcinoma.

So if every 1 gram sample requires two tumors of almost 10 mm diameter to get
$4.5\pm 0.9\times 10^{7}$ cancer cells and $2.2\pm 0.6\times 10^{6}$ TAMs, we
will have the following estimates for our data:

  * •

$1.65\pm 0.33^{7}$ cancer cells and $8.09\pm 0.22\times 10^{5}$ TAMs during
hyperplasia.

  * •

$2.08\pm 0.41\times 10^{7}$ cancer cells and $1.02\pm 0.27\times 10^{6}$ TAMs
during adenoma/MIN.

  * •

$2.51\pm 0.50\times 10^{7}$ cancer cells and $1.22\pm 0.33\times 10^{6}$ TAMs
during early carcinoma.

  * •

$3.17\pm 0.63\times 10^{7}$ cancer cells and $1.54\pm 0.42\times 10^{6}$ TAMs
during late carcinoma.

In this study, the number of TAMs adds up to $M_{N}+M$ from equations (8) and
(9). With the cell fraction values available (Figure 2) and knowing the total
number of macrophages, we can promptly calculate the other immune cell
quantities.

As for the non-immune cells, we use:

$\sum_{i\in I}N_{i}V_{i}+\sum_{j\in J}M_{j}V_{j}+N_{C}V_{C}=\textrm{Total
Volume}.$ (20)

Where $N_{i}$ and $M_{j}$ are the numbers of $i$-th immune cell and $j$-th
non-immune cell, respectively. The sets $I$ and $J$ contain indices for immune
and non-immune cells (Adipocytes and ECs). The values $V_{i}$ and $V_{j}$ are
the single cell volume of the corresponding immune and non-immune cells. The
values $N_{C}$ and $V_{C}$ are the number and volume of cancer cells. Recall
that the reported tumor volumes are devoid of necrotic cells. Using table 2
with the tumor volumes mentioned, we can get the number of non-immune cells
using formula (20).

Tables 3 and 4 show the calculated number of cells along with the value of
cytokines inferred directly from the gene expression data. For the oxygen
levels, Ebbesen et al. report the following [34]:

  * •

10000-50000 ppm is normoxic (equivalent to 1-5% $O_{2}$ per kPa)

  * •

1000-10000 ppm is mild hypoxia when HIF pathways become significant.
(equivalent to 0.1-1% $O_{2}$ per kPa)

  * •

10-1000 ppm is severe hypoxia when cell death or cell cycle arrest is
triggered (equivalent to 0.001-0.1% $O_{2}$ per kPa)

Noting the oxygen molar mass is 16 g/mol, we can turn the above values into
mol/L values. Cai et al., whose data we use, speculate the decrease in the
expression of tricarboxylic acid cycle (TCA cycle) genes in later cancer
stages is due to hypoxia. Therefore, the levels given for oxygen in table 4
are picked such that the later stages go through hypoxic conditions.

All the data is non-dimensionalized before associating them with our model.
The non-dimensionalization is done by dividing each cell or cytokine value in
tables 3 and 4 by the largest value of that cell or cytokine across all
samples and all stages. Also, in simulations, we assume our week 6 is the
system’s initial state, corresponding to $t=0$.

Cell type | Length∗ ($\mu m)$ | Reference
---|---|---
Dendritic cells | 10-15 | Tasnim et al.[118]
Macrophages | 21 | Chitu et al. [23]
Granulocytes and Monocytes | 13-15 | Kornmann et al. [61]
Lymphocytes (i.e., T-cells, B-cells, and NKs) | 10-15 | O’Connell et al. [96]
Stem cells | 7-17 | Pillarisetti et al. [100]
ECs | 20 | Jiang et al. [56]
Luminal epithelial cells of mammary gland | 10-50 | Chen et al. [21]
Cancer cells | 12 | Laget et al. [62]
Adipose-tissue derived stem cells | 28 | Hoogduijn et al. [54]
Fibroblasts∗∗ | 10-20 | Turgay et al. [125]
Adipocytes | 20 | Hagberg et al. [47]
Table 2: Size of different cell types in mice. We use the midpoint in our
calculations for the lengths with a given range.
* We assume the major (L) and minor (W) diameter of the cells are roughly equal and use the caliper format $L\times W^{2}/2(=L^{3}/2)$ to get the volume.
** We have not imputed cell fractions for fibroblasts. However, since most
mammary stromal cells are adipocytes and fibroblasts, we use this fact to take
the intersection of their approximate size (20 $\mu m$) as the size of stromal
cells.

Sample: Stage | $T_{N}$ | $T_{h}$ | $T_{C}$ | $T_{r}$ | $D_{N}$ | $D$ | $M_{N}$ | $M$ | $C$ | $N^{*}$ | $A$ | $E_{N}$
---|---|---|---|---|---|---|---|---|---|---|---|---
M1: W6 | 1382527 | 663086 | 228615 | 92289 | 1227231 | 886534 | 389295 | 397821 | 13200000 | 6600000 | 416120 | 915268
M1: W8 | 1615524 | 641814 | 265967 | 176412 | 1574935 | 1057404 | 577007 | 710944 | 16700000 | 8350000 | 1582762 | 916916
M1: W10 | 2569211 | 689479 | 310204 | 292578 | 1636578 | 1328095 | 612205 | 942441 | 20100000 | 10050000 | 2763480 | 769340
M1: W12 | 2199247 | 937204 | 222401 | 147089 | 2324486 | 1686078 | 787834 | 1171025 | 25400000 | 12700000 | 5222452 | 3028739
M2: W6 | 1441964 | 529513 | 128158 | 79395 | 1115894 | 934727 | 450621 | 358170 | 16500000 | 8250000 | 461103 | 493492
M2: W8 | 1545885 | 583058 | 589238 | 407590 | 1152659 | 1411717 | 292503 | 455659 | 20800000 | 10400000 | 1080624 | 185467
M2: W10 | 1311444 | 478978 | 183104 | 69789 | 1243110 | 859137 | 430970 | 458725 | 30100000 | 15050000 | 2464786 | 931836
M2: W12 | 2297591 | 426244 | 127776 | 181325 | 1396820 | 1057303 | 559029 | 562339 | 31700000 | 15850000 | 5578533 | 1469427
M3: W6 | 1399198 | 430269 | 405214 | 79316 | 1021227 | 924575 | 395039 | 435495 | 19800000 | 9900000 | 333365 | 366687
M3: W8 | 1467427 | 507931 | 395476 | 81323 | 1268459 | 836283 | 439574 | 575209 | 24900000 | 12450000 | 1514664 | 785294
M3: W10 | 2161705 | 340691 | 272777 | 188039 | 1232740 | 1346494 | 554893 | 664238 | 25100000 | 12550000 | 2885384 | 1196948
M3: W12 | 3171264 | 963737 | 381975 | 162666 | 1909543 | 1177989 | 658718 | 884036 | 38000000 | 19000000 | 5437350 | 2937968

Table 3: Number of cells. The numbers in each cell are integer cutoffs of the
calculated values.
* We assume the number of necrotic cells is half of cancer cells. This is a crude approximation based on observations from the literature about the percentage of viable space in the TME occupied by necrotic cells [137].  Sample: Stage | $H$ | $IL_{12}$ | $IL_{10}$ | $IL_{6}$ | $V$ | $Ox$
---|---|---|---|---|---|---
M1: W6 | 1000 | 63 | 417 | 384 | 776 | 3.125
M1: W8 | 940 | 41 | 351 | 197 | 536 | 0.625
M1: W10 | 1103 | 3 | 404 | 184 | 494 | 0.0625
M1: W12 | 1050 | 2 | 455 | 193 | 798 | 0.000625
M2: W6 | 1182 | 0 | 450 | 196 | 775 | 3.125
M2: W8 | 932 | 16 | 723 | 35 | 721 | 0.625
M2: W10 | 945 | 0 | 429 | 162 | 313 | 0.0625
M2: W12 | 807 | 4 | 319 | 127 | 384 | 0.000625
M3: W6 | 1521 | 27 | 511 | 336 | 1081 | 3.125
M3: W8 | 1549 | 22 | 566 | 312 | 1054 | 0.625
M3: W10 | 957 | 3 | 349 | 229 | 565 | 0.0625
M3: W12 | 779 | 10 | 278 | 163 | 909 | 0.000625
Table 4: Cytokine levels. All the levels are, in fact, gene expression levels
except for oxygen. The unit for oxygen concentration is mol/L.

### Parameter estimation

We use an HGA scheme based on natural evolution for parameter estimation. We
utilize MATLAB’s HGA, which is included in the Global Optimization Toolbox.
First, we initialize a loss function that calculates the distance between the
predicted values and the measured data.

We then start with a randomly produced set of parameters called the initial
population. This can be constrained, biased, or a free random pick. For our
problem, we do not accept negative parameter values. However, after several
tests, we saw that the algorithm does not produce parameter values bigger than
2 for large upper bounds. So, we take the upper bound to be 2.

In each iteration, the algorithm leverages the individuals in the present
generation to construct the succeeding population. The creation of this new
population involves the following steps:

  * •

Assessing each member’s fitness in the current population yields the raw
fitness scores. This is done by utilizing the loss function.

  * •

Scaling these raw fitness scores into a more practical range of values,
referred to as expectation values.

  * •

Selecting parents and individuals chosen based on their expectation values.
MATLAB’s default selection method constructs a line where each parent is
associated with a line segment. The scaled value of the respective parent
determines the length of each segment. The algorithm advances along the line
in uniform steps, and at each step, it assigns a parent from the segment it
arrives at. This way, we can pick out the elite children.

  * •

Identifying certain individuals in the current population with lower fitness
as elite members carried over to the next population.

  * •

Generating offspring from the selected parents through random alterations to a
single parent (mutation) or combining the attributes of a pair of parents
(crossover).

  * •

Substituting the existing population with the newly created offspring to
constitute the subsequent generation.

  * •

The algorithm stops when one of the stopping criteria is met.

  * •

A local search algorithm starts using the best solution found by the genetic
algorithm as the initial guess. We use the widely known MATLAB function
fmincon, a gradient-based nonlinear constrained local minimizer.

The last step is only adopted in HGA. The choice of the local minimizer
depends on the problem. For the theory behind the genetic algorithm, we refer
the readers to Holland’s book [53]. Also, for a review of the genetic
algorithms and different hybrid versions, see the article by El-Mihoub et al.
[35].

Table 1 in the supplementary material shows the non-dimensional estimated
values of the parameters for mice 1,2 and 3. As mentioned, non-
dimensionalization is done by scaling each cell and cytokine level by the
highest value across all three mice and all time points. This will enable us
to compare the three mice and make the parameter estimation more stable. Note
that we do not non-dimensionalize the time scale.

### A discussion about the model’s identifiability

We recognize that the assumptions and connections we made with the literature
to obtain the data may be optimistic. For example, the cell sizes we used to
infer the number of non-immune cells can add unforeseen uncertainties to our
findings. However, our framework does result in plausible values of cells and
cytokines, and if all of the above can be reported under the same study, then
this can be a meaningful starting point for further studies and refinements.

Assuming all the model variables are practical outputs (i.e., are available
for all times), our model and its parameters are Full Input-State-Parameter
Observable (FISPO), a term coined by Villaverde et al. [129]. One can evaluate
the global structural identifiability of nonlinear systems by employing
concepts from differential geometry. This assessment primarily relies on
computing the rank of a specialized matrix known as the observability-
identifiability matrix, formed by utilizing Lie derivatives. If the parameters
are included in this computation as zero-dynamic variables, then the
observability analysis of the parameters will be equivalent to their
identifiability. For more details of these methods, see [129, 127, 128, 78,
77].

We use the MATLAB package STRIKE-GOLDD to ascertain the global structural
identifiability of our model’s parameters. As previously noted, STRIKE-GOLDD
declares all the model parameters structurally identifiable, assuming all our
variables can be regarded as practical outcomes (i.e., data is available for
them at all times).

A parameter can be globally structurally identifiable but not practically
identifiable. [136]. A parameter that minimizes the loss function (or
maximizes the likelihood function) but has an infinite confidence interval is
considered a practically non-identifiable parameter [102]. To investigate
these parameters within the estimated range, we create synthetic data by
perturbing the solutions in Figure 3 by 5%. Then we calculate the profile
likelihood (see [102]) of each parameter, i.e., $PL(\theta_{i})$. This is done
by keeping the corresponding $\theta_{i}$ fixed and minimizing the loss
function with respect to the rest of the parameters ($\theta_{j}$ where $j\neq
i$). Next, we vary $\theta_{i}$ for each calculated profile likelihood and
record the changes in the profile likelihood function. We then calculate a
threshold for each set of estimated parameters by
$L(\boldsymbol{\theta})+\Delta_{\alpha}$, where $L$ is the loss function used
for parameter estimation and $\Delta_{\alpha}$ is the $\alpha$-quantile of the
$\chi^{2}$ distribution with $df=1$ for pointwise confidence intervals.
Profile likelihoods not exceeding this threshold will be labeled as
practically non-identifiable within the specified parameter range.
Additionally, parameters with very flat profile likelihood functions may have
local structural identifiability issues around the estimated value. Refer to
Figures 1, 2, and 3 in the supplementary materials for results. Within the
investigated parameter neighborhoods ($[0,2*\theta_{i}]$), mouse 1 has 41
practically non-identifiable parameters, with 31 displaying visibly flat
profiles, potentially indicating local structural non-identifiability. These
numbers are lower for mice 2 and 3, which aligns with their Principle
Component Analysis (PCA) count of sensitive variables. Mouse 2 and 3 have 35
parameters that are practically non-identifiable, with 12 from mouse 2 and 13
from mouse 3 having visibly flat profile likelihoods within the studied
parameter range. This shows these parameters are likely to remain non-
identifiable even with increased data points or measurement quality. Notably,
none of these parameters with flat likelihood profiles are identified as
sensitive parameters for each mouse, suggesting their values won’t
significantly impact system robustness. However, we acknowledge limitations in
this approach, as it relies on synthetic data and assumes all variables are
always observable. However, the method can be adapted to accommodate feasible
experimental situations.

### Sensitivity analysis

Our study employs a Principal Component Global Sensitivity Analysis (PCGSA) to
scrutinize the parameter’s impact on our model. To accomplish this, we utilize
the Matlab package PeTTSy. The methodology for calculating sensitivity
measures is derived from the approach outlined by the PeTTSy developers in
their seminal paper [31]. We provide a summary of their technique here. We
strongly encourage readers to consult the original study for a more
comprehensive understanding and in-depth insights, as our description serves
only as an overview of their detailed method.

Suppose $X_{m}(t,{\bf k})$ is a solution of an ODE system with $t$ being the
time variable and ${\bf k}=(k_{1},k_{2},\cdots,k_{N})$ being the parameters
vector. So the i-th solution derivative with respect to the j-th parameter at
a time $t$ is given by $\frac{\partial X_{m}(t,{\bf k})}{\partial k_{j}}(t)$.
A Singular Value Decomposition (SVD) analysis is carried out to calculate the
global sensitivity. This is to investigate the linearization of the mapping
from perturbed parameters $\delta{\bf k}\in R^{N}$ to the changes they cause
in the solution, namely $\delta{\bf X}$. Notice the mapping is to an infinite
dimensional space of smooth functions. PeTTSy approximates these functions by
evaluating them on a fine time grid (with $n$ time steps), which results in
high-dimensional vectors. The final product will be a large matrix
$M=\frac{\partial{\bf X}}{\partial{\bf k}}$ with $nK$ rows and $N$ columns,
where $K$ is the number of our state variables. Like the state variables, we
can include any other measurements, such as the sum of all the cells or the
sum of all immune cells, as a single variable for the sensitivity analysis.

The issue is that this matrix is relatively large, and investigating the
parameters effect on the time-variant of a large set of solutions can be
tedious. This is where the SVD and PCA are helpful. They help us focus only on
a handful of Principal Components (PCs) to understand the overall parameters
effect on the whole model. We start with a singular value decomposition of
$M$.

$M=U\Sigma V^{T}$ (21)

Where $U$ is an $nK\times nK$ orthogonal matrix, $\Sigma$ is an $nK\times N$
rectangular diagonal matrix with singular values $\sigma_{i}$ on its diagonal,
and $V$ is a matrix whose columns form an orthonormal basis for the parameter
space. The principal global sensitivity values are $S_{ij}=\sigma_{i}V_{ji}$.
They show

$\delta{\bf X}=\sum_{i,j}S_{ij}\delta k_{j}\sigma_{i}U_{i}+O(||\delta
k^{2}||),$ (22)

where $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{N}\geq 0$ are the
singular values and $U_{i}$ (the columns of $U$) are PCs.

For the parameter sensitivity plots, we use the values
$S_{ij}=\sigma_{i}V_{ji}$ for each PC $U_{i}$. To determine how many PCs are
needed to get a good sense of the system’s overall sensitivity to parameters,
we utilize a singular spectrum plot. We filter out the PCs corresponding to
insignificant singular values.

We also investigate the sensitivity of cancer cells and the total number of
immune cells to the model parameters. Given that we only get one significant
singular value for these two cases, the PC sensitivity and the classical
notion of sensitivity will be synonymous according to the relationship (22).

##  3 Results

### Dynamics

We solve the non-dimensional ODEs (2)-(19) corresponding to Figure 1 with
parameters given in supplementary materials (Table 1). The results are shown
in Figure 3; the corresponding findings illustrate a good agreement with the
data. The mouse data is cut off after 12 weeks due to the subject’s
euthanizing, but we follow the dynamics until day 150. Moreover, the tumor
samples are collected starting at week 6, which in the models marks $t=0$. Due
to the system’s complexity, we comment primarily on the dynamics related to
the model extension, namely the ones involving ECs, VEGF, and oxygen. Further
investigation will be carried out through sensitivity analysis.

Figure 3: Cells and cytokine dynamics. The x-axes show time in days. The
y-axes are the dimensionless levels.

Figure 3 shows a qualitative resemblance between the macrophage and VEGF
dynamics. Mouse 1 and 2 have the highest and lowest number of macrophages,
respectively, consistent with adipocyte numbers, especially after $t=40$.
Notice that this is true in numerical observations, even though not manifested
in the experimental data for earlier times. The causal effect is expected as
macrophages and adipocytes are primary VEGF producers in our model, besides
hypoxic cancer cells and mild hypoxia-induced pathways. Initially, VEGF levels
decrease due to higher oxygen levels delaying hypoxic cancer cells and other
hypoxia mechanisms contributing to its production. However, despite this
decrease, endothelial cells (ECs) increase in mice 2 and 3, attributed to
lower initial IL12 levels in these subjects compared to mouse 1. Consequently,
even with a higher inhibition rate, ECs do not significantly decrease due to
their low IL12 levels.

Having the lowest endothelial cells (ECs), mouse 2, experiences the most rapid
decline in oxygen levels, maintaining the lowest oxygen values throughout the
simulation. However, during the time-period from 0 to 20, where oxygen
decrease is significant, the blue curve representing $E_{N}$ for mouse 2
surpasses that of mouse 3. This phenomenon is attributed to IL12 levels as
discussed earlier. Despite starting with more oxygen producers ($E_{N}$),
mouse 2 experiences a more pronounced oxygen decrease compared to mouse 3 due
to a higher number of oxygen-consuming cell types during this short interval.
In addition to ECs, mouse 2 has more T-cells, activated dendritic cells, naive
macrophages, and adipocytes during this period, see Figure 3.

Mouse 3 has the most cancer cells, followed by Mouse 1 for long enough times.
This could be influenced by various molecular, cellular, and hypoxic pathways,
with cytotoxic T-cells enhanced by HIF pathways playing a significant role.
Mouse 2, having the highest cytotoxic T-cells and the lowest adipocytes, shows
the least cancer cell growth. The next section will explore more analysis of
the key pathways through sensitivity analysis.

### Sensitivity Analysis Results

We analyze the sensitivity of the entire system, cancer cells, and total
immune cell count ($T_{N}+T_{h}+T_{C}+T_{r}+D_{N}+D+M_{N}+M$) to model
parameters. Due to numerous state variables, we conduct a PCGSA for the first
analysis. Classical global sensitivity and first PC sensitivity are equivalent
for the latter two. First, we conduct a comprehensive PCGSA of the entire
system relative to its parameters. With 89 parameters and 18 variables,
analyzing all 18 PCs may be redundant. To determine the necessary PCs, we
perform a singular spectrum analysis. Figure 4 displays each mouse’s top 10
normalized singular values. In mouse 1, the eighth singular value is
approximately 20% of the first. For mice 2 and 3, this occurs around the fifth
and sixth singular values, respectively. Therefore, we focus on sensitivity to
variations in the top 8 PCs for each mouse.

Figure 4: Singular spectrum. Top 10 largest singular values for Mouse 1,2 and
3.

Figure 5 shows state variables more sensitive to parameter variations when
different PCs are considered, illustrating $\sigma_{i}U_{i}$ coefficients in
(22). Only the variables whose sensitivity exceeds 20% of the global maximum
sensitivity are included. There are sensitive state variables for all 8 PCs
considered in the analysis for mouse 1. However, mice 2 and 3 only have
variables exceeding 20% of the global maximum up to the fifth and sixth PCs,
respectively. This was expected based on the observations in Figure 4. We
notice the dominant presence of immune cells for each PC in all three
subjects. Cells such as regulatory and helper T-cells, macrophages, and
dendritic cells are among the most repeated variables. These cells are all
involved in producing IL10, the most sensitive cytokine in the plot. This plot
can also help in designing experiments for validation. The variables shown in
this plot and, more specifically, the ones shared in all three subjects are
appropriate choices for experimental observations.

Figure 5: State variables sensitivity heat maps with respect to time. Only
variables above 20% of the maximum global sensitivity are shown for each
mouse. For each PC, the represented variables are sorted according to their
scaled ratio of change. Figure 6: PCGSA of the top 10 most sensitive
parameters for (A) Mouse 1, (B) Mouse 2, and (C) Mouse 3.

Figure 6 shows the PCGSA result for mice 1, 2, and 3. This plot has to be
accompanied by Figure 5 for clearer biological inferences. Parameters present
in all three plots are $\lambda_{VE_{N}}$, $A_{M_{N}}$ and $A_{D_{N}}$.
Changing $\lambda_{VE_{N}}$ inflicts the highest variation in the first PC in
all three mice. Figure 5 shows immune cells such as regulatory and helper
T-cells, macrophages, and dendritic cells are the most susceptible to this
parameter variation. This parameter describes the promotion rate of EC
proliferation by VEGF. One can connect this to the rate at which vessels are
formed in the microenvironment. Also, $\lambda_{VE_{N}}$ involves two
variables from the new model extension; see Figure 1.

Parameters $A_{M_{N}}$ and $A_{D_{N}}$ are the innate naive macrophages and
naive dendritic cell proliferation rates. For mouse 1, $A_{M_{N}}$ affects PC
6 and 3, with naive macrophages and HMGB1 receiving the highest impact. In
mice 2 and 3, PCs 1, 2, 3, and 5 are affected the most, with naive and
activated macrophages and helper and regulatory T-cells being the most
sensitive variables to this parameter. Finally, besides the dendritic cell
subtypes, helper and naive T-cells and activated and naive macrophages are
included in the PCs, which are significantly influenced by $A_{D_{N}}$
variation. However, these effects are expected, as they control the production
of two crucial naive immune cells and, consequently, their activated form.
These cell types are involved in some of the most central interactions within
the system.

Oxygen-related parameters $Ox_{crit}^{(2)}$, $\lambda_{OxE_{N}}$, and
$\delta_{5}$ are among the top ten most sensitive parameters in mice 1 and 3.
All three show significant sensitivity for the first PC, greatly affecting
helper and regulatory T-cells, macrophages, and dendritic cells. The parameter
$Ox_{crit}^{(2)}$ decides the transitioning threshold from mild to severe
hypoxia, $\lambda_{OxE_{N}}$ describes how much oxygen is provided by
ECs/vessels, and $\delta_{5}$ is the oxygen uptake rate by adipocytes. This is
interesting since in our other studies, parameters related to adipocytes had a
considerable effect on the model [85, 86]. Adipocytes are abundant in breast
tissue, and adiposity is considered a well-known risk factor [112]. Besides
affecting cancer directly, this suggests adipocytes also affect the
microenvironment immune profile through oxygen-related pathways. Moreover,
$\lambda_{A}$ and $\delta_{A}$, which are adipocyte production and death
rates, are placed among sensitive parameters; see Figure 6 and tables 2,3 and
4 in the supplementary material.

We should point out that IL-12 is affecting mouse 2’s first PC by inhibiting
ECs ($\delta_{E_{N}IL_{12}}$). It has a similar effect on mouse 3 by promoting
the cytotoxic cells ($\lambda_{T_{C}IL_{12}}$). So, despite not being among
the sensitive state variables in Figure 5, IL-12 bears a considerable impact
on the immune profile through the PC1 sensitive cells such as $T_{h}$,
$M_{N}$, $M$ and $T_{r}$. Finally, in mouse 2, the VEGF production by
macrophages ($\lambda_{VM}$) seems to have a similar effect through PC1 on the
system.

Overall, the PCGSA shows the importance and intricacy of the immune system
interactions in the TME. All parameter variations, especially parameters
involving vasculature and oxygen delivery from our model extension,
substantially influence crucial immune cells. We now examine whether the
previously discussed parameters significantly influence diagnostic metrics,
specifically the cancer and the total immune cell counts. Figure 7 illustrates
the sensitivity analysis focusing on these crucial measures.

As expected, cancer levels are sensitive to parameters directly involved in
the equation for their evolution (see Eq. (10)) such as $\delta_{C}$,
$\lambda_{C}$, $\lambda_{CA}$ and $\lambda_{CIL_{6}}$. This result agrees with
the findings from the simpler model [85]. Also, when calculating the total
immune cells’ sensitivity to parameters, immune cells’ innate production and
death rates appear among the most sensitive parameters.

In Figure 7, $\lambda_{VE_{N}}$, $Ox_{crit}^{(2)}$ and $\lambda_{OxE_{N}}$ are
the most repeated, warranting particular attention in our study. They were
present in the PCA and are, importantly, a part of the model extension; see
Figure 1. The results show that an increase in $\lambda_{VE_{N}}$ leads to an
increase in all subjects’ immune cells and cancer cells. This might be due to
less hypoxic death for cancer cells or more pro-tumor cytokine production by
immune cells such as macrophages. Increasing $Ox_{crit}^{(2)}$ causes a
quicker transition to severe hypoxia, decreasing the total number of immune
and cancer cells. This effect can be due to several factors, such as an
increase in cytotoxicity or a decrease in IL-6 production. Note that the
result for cancer sensitivity to $Ox_{crit}^{(2)}$ is not included in Figure 7
for mouse 2 since it is the 13th most sensitive parameter. Nevertheless, the
value is negative like the others.

The oxygen delivery rate by the ECs/vessels ($\lambda_{OxE_{N}}$) is
positively related to the total number of immune and cancer cells. Again, this
is not included in the figure for mouse 2 since it is the 11th most sensitive
parameter. The model extension effect on cancer cells through these parameters
and the immune system showcases the importance of angiogenesis and oxygenation
in the TME. This, in turn, introduces a limitation to our work regarding the
cancer sensitivity to $Ox_{crit}^{(2)}$. Even though the hypoxic contribution
to cancer proliferation through HIF pathways is included, this study considers
purely oxidative phosphorylation for cells and not glycolysis as a metabolism
approach. However, the metabolic switch to glycolysis under hypoxic stress is
a common adaptive response from cancer cells and a challenge in cancer
research [29]. Therefore, increasing $Ox_{crit}^{(2)}$ might not necessarily
lead to fewer cancer cells if it induces a quicker metabolic change. These
intricate mechanisms require deeper analysis and a more complicated model,
which we defer to future research.

The relationship between adipocyte-related parameters and cancer and total
immune cells can be explained by investigating the sensitivity value signs for
$\lambda_{A}$, $\delta_{A}$, $\delta_{5}$, and $\lambda_{CA}$. In summary, if
we increase the adipocyte number through increasing $\lambda_{A}$ or
decreasing $\delta_{A}$, the total number of immune cells decreases, and of
cancer cells increases. Increasing adipocytes greatly increases oxygen
consumption, which is lethal for all cells affected by the hypoxic cell death,
including immune cells. However, this effect is foiled for cancer cells by the
direct adipocytes effect on cancer proliferation through $\lambda_{CA}$.

An increase in VEGF leads to more immune cells and more cancer cells. This can
be done through decreasing $\delta_{V}$ (VEGF decay rate) in mice 1 and 2 and
increasing $\lambda_{VM}$ (VEGF production by macrophages) in mouse 2. The
same effects occur for mouse 3 but are not included among the top 10 in figure
7.

Figure 7: Sensitivity analysis with respect to cancer and total immune cells.
Top 10 most sensitive parameters with respect to (left) cancer cells, (right)
total immune cells.

We also investigate the effect of IL-12-related parameters. Inhibiting the
IL-12 production (increasing $\delta_{IL12}$) leads to more immune cells and
cancer cells in mice 2 and 3. Less IL-12 means less cytotoxic T-cell
activation, a major cancer cell inhibitor. Also, less IL-12 means more EC
proliferation (see $\delta_{E_{N}IL_{12}}$), which results in more oxygen,
promoting the proliferation of all immune and cancer cells. See the effect of
its increased production by dendritic cells ($\lambda_{IL_{12}D}$) on cancer
and immune cells in mice 2 and 3 in Figure 7. Lastly, $\lambda_{T_{C}IL_{12}}$
is negatively associated with cancer cells and total immune cells in mouse 3.
This is the cytotoxic T-cell activation rate by IL-12, one of the primary
cancer cell inhibitors.

We finally explore the effect of perturbing $\lambda_{VE_{N}}$,
$Ox_{crit}^{(2)}$ and $\lambda_{OxE_{N}}$ on the cancer dynamics and total
immune cells given their presence in all three sensitivity analyses for three
mice. These parameters are associated with our model extension and include the
effect of all three new compartments: ECs, VEGF, and oxygen. Figure 8 shows
the effect of the aforementioned parameters’ perturbation on cancer and total
immune cell dynamics. For $\lambda_{VE_{N}}$, the largest impact is between
$t=30$ and $t=80$, which wears off as cancer approaches a steady state. Cancer
steady-state translates into a state in which cancer cannot physically grow
any further due to space or resources depletion. The impact becomes visible
for the total immune cells between times 20 and 30. For mouse 1, it is a
short-lived effect wearing off around time 70. For mouse 3, it takes longer
but eventually tapers after time 100, and Mouse 2 shows a consistent regime
throughout the simulation time. This observation implies controlling cancer
through controlling angiogenesis and VEGF inhibitors should happen in earlier
stages to have a positive impact. Once cancer reaches its steady state, these
variations do not leave a significant impact. It has been discussed in the
literature for humans and mice that VEGF therapy for late-stage breast cancer
has not been correlated with better prognoses [8].

Perturbing $\lambda_{OxE_{N}}$ and $Ox_{crit}^{(2)}$ results in a sustainable
but less significant effect. Especially in mice 1 and 3, the impact region for
total immune cells increases in time. Even though controlling hypoxia or
oxygen delivery is not medically feasible, numerous studies focus on targeting
the HIF pathways for treatment [79, 97, 107].

Figure 8: Effect of parameter perturbation on dynamics. The effect of
perturbing $\lambda_{VE_{N}}$ (top), $\lambda_{OxE_{N}}$ (middle) and
$Ox_{crit}^{(2)}$ (bottom) on dynamics of cancer cells and total immune cells.

##  4 Discussion

In this study, we extended an existing ODE model, as detailed in [85], to
investigate the mechanisms that promote or are angiogenesis outcomes within
the PyMT-MMTV mouse breast TME. We added three important elements: ECs, VEGF,
and oxygen. ECs, forming the blood vessels’ infrastructure, are crucial for
transporting metabolites like oxygen. VEGF, a key cytokine, initiates and
accelerates angiogenesis. It is secreted by various cells, including
macrophages, adipocytes, and hypoxic cancer cells. Additionally, ECs’ activity
is modulated by IL-12. By incorporating these components, we effectively
linked them with the previously developed model [85], thereby investigating
the effect of angiogenic processes and hypoxia in breast TME.

Using PyMT-MMTV mice model data and HGA parameter estimation, we approximated
the ODE model parameters and simulated results for 150 days. Significant
differences emerged in dynamics: mouse 1 showed the highest cancer levels,
while mouse 2 showed the lowest levels attributed to cytotoxicity and
adiposity. Other notable interactions included high levels of helper T-cells
and dendritic cells in mouse 1, linked to HMGB1, and elevated activated
macrophages due to IL-10 and helper T-cells. Also, the number of ECs
correlated with VEGF levels. However, some pathways couldn’t be directly
inferred from the interaction network, requiring global sensitivity analysis
experiments for further investigation.

Through sensitivity analyses, we observed that similar to a previous study,
adipocyte-related parameters (production and death rates) were crucial to the
system [85]. Additionally, the oxygen consumption by adipocytes was found to
be sensitive. Several studies have highlighted adipocytes’ effect on hypoxia
in the TME. Nieman et al. argue that in individuals with excessive adipocytes,
the adipose tissue is associated with chronic hypoxia, contributing to tumor
progression through HIF pathways [93]. Trayhurn reasons the adipocytes’ sheer
mass and abundance is a hypoxia-promoting factor [122].

Furthermore, a subset of parameters consistently emerged as influential
factors affecting the whole system’s PCs and the dynamics of cancer cells and
the immune system in all subjects. These parameters, which describe the ECs
promotion by VEGF ($\lambda_{VE_{N}}$), oxygen delivery by the ECs/vessels
($\lambda_{OxE_{N}}$), and the threshold value for transitioning from mild to
severe hypoxia ($Ox_{crit}^{(2)}$), are all related to angiogenesis and
oxygenation. The first parameter introduces an intervention opportunity for
controlling cancer progression and modulating the immune response. In
contrast, $\lambda_{OxE_{N}}$ depends more on the cells’ fitness landscape and
tissue content. In all three mice, increasing $\lambda_{VE_{N}}$ led to higher
cancer cell counts and greater total immune cell numbers, suggesting
angiogenesis plays a crucial role in the TME. Additionally, changing the
threshold $Ox_{crit}^{(2)}$ influenced the immune cell dynamics and cancer
progression. These findings underscore the importance of oxygen levels in the
microenvironment and suggest targeting hypoxia-related pathways could be a
promising avenue for therapy. Recall, however, that the cancer cell metabolic
plasticity is a relevant factor [17] that merits further study. The oxygen
delivery rate by ECs/vessels ($\lambda_{OxE_{N}}$) consistently impacted
cancer and immune cells in all three mice. This shows the critical role of
efficient oxygen transport in regulating the immune response within the TME.
It also highlights the importance of angiogenesis in oxygen supply, which can
affect both immune cell function and cancer cell proliferation.

We vary the three parameters discussed to study their effects on cancer and
total immune cells. Perturbing $\lambda_{VE_{N}}$ shows a time-sensitive
effect, showing the early intervention significance in managing cancer
progression. The short-lived impact on total immune cells in mouse 1 and the
longer-lasting effect in mouse 3 highlight the nuanced responses in different
scenarios. This observation aligns with existing literature cautioning against
late-stage VEGF therapy for breast cancer, emphasizing timely intervention to
yield positive outcomes [8]. Perturbing $\lambda_{OxE_{N}}$ and
$Ox_{crit}^{(2)}$ has a less significant but sustainable effect on total
immune cells, especially in mice 1 and 3 over a longer period. While directly
controlling hypoxia or oxygen delivery may be challenging, the lasting impact
suggests promising therapeutic avenues, as seen in studies targeting the HIF
pathways for treatment. [79, 97, 107].

Our results connect with several clinical trials that have demonstrated the
VEGF inhibitors efficacy, such as Bevacizumab, in breast cancer therapy [80].
These inhibitors reduce tumor vasculature, disrupt angiogenesis, and
indirectly influence oxygenation levels. Importantly, this approach can
potentially enhance immuno-therapies efficacy by alleviating hypoxia-induced
immunosuppression [55]. By targeting processes summarized as
$\lambda_{VE_{N}}$ in our model, these therapeutic approaches underscore the
potential clinical relevance of our findings. They also highlight the
importance of personalized treatment strategies that consider the specific
angiogenesis and oxygenation dynamics in breast cancer patients. Our work
paves the way for the potential in-silico exploration of the effects of such
inhibitors in future studies.

This study’s findings should be interpreted considering its limitations. One
major constraint is the lack of refined time-course data, which, if available,
could enhance model precision and prediction reliability. Exploring spatially
distributed data and comparing a spatiotemporal version of our model with such
data would be extremely interesting. Despite having a sizable system, we
acknowledge some key role players in TME are not considered. This constitutes
an opportunity for more extensions depending on the particular research needs.
Additionally, numerous factors such as the role and type of macrophages [13],
as well as the cells’ phenotypic variation under different oxygen conditions
[17] are among factors worth further considering in the context of the present
network. Lastly, we emphasize the potential metabolic changes occurring in
cancer cells, making them nearly immune to hypoxia. This is known as the
Warburg effect [67] and is a HIF pathways’ outcome [28]. Despite including the
HIF pathways’ effect on the immune system, we did not consider the pathways
leading to the Warburg effect for simplicity.

Future studies are needed to overcome the limitations highlighted above, like
seeing how the metabolic changes in cancer can affect the system outcome.
Moreover, the direct interactions between cancer cells and macrophages can be
included in the model by considering the macrophages’ polarization to pro and
anti-tumor subtypes rather than just through the molecular pathways.
Furthermore, this model lays the groundwork for developing a TME spatio-
temporal PDE model. Such a model can capture interesting phenomena such as
VEGFs’ chemotactic signaling to endothelial cells affecting the vessel network
formation and location and the oxygen penetration in the tissue influencing
the dynamics and shape of the TME.

## References

  * [1] Kasper Almholt, Anna Juncker-Jensen, Ole Didrik Lærum, Keld Danø, Morten Johnsen, Leif Røge Lund, and John Rømer.  Metastasis is strongly reduced by the matrix metalloproteinase inhibitor galardin in the mmtv-pymt transgenic breast cancer model.  Molecular cancer therapeutics, 7(9):2758–2767, 2008.
  * [2] Simon Anders and Wolfgang Huber.  Differential expression analysis for sequence count data.  Nature Precedings, pages 1–1, 2010.
  * [3] A.R. Anderson and M.A. Chaplain.  Continuous and discrete mathematical models of tumor-induced angiogenesis.  Bull Math Biol., 60:857–869, 1998.
  * [4] Shweta Aras and M Raza Zaidi.  Tameless traitors: macrophages in cancer progression and metastasis.  British journal of cancer, 117(11):1583–1591, 2017.
  * [5] Melina Arnold, Eileen Morgan, Harriet Rumgay, Allini Mafra, Deependra Singh, Mathieu Laversanne, Jerome Vignat, Julie R Gralow, Fatima Cardoso, Sabine Siesling, et al.  Current and future burden of breast cancer: Global statistics for 2020 and 2040.  The Breast, 66:15–23, 2022.
  * [6] Rachel A. Aronow, Shaya Akbarinejad, Trang Le, Sumeyye Su, and Leili Shahriyari.  TumorDecon: A digital cytometry software.  SoftwareX, 18:101072, jun 2022.
  * [7] Sherif Attalla, Tarek Taifour, Tung Bui, and William Muller.  Insights from transgenic mouse models of pymt-induced breast cancer: recapitulating human breast cancer progression in vivo.  Oncogene, 40(3):475–491, 2021.
  * [8] Gabriele Bergers and Douglas Hanahan.  Modes of resistance to anti-angiogenic therapy.  Nature Reviews Cancer, 8(8):592–603, 2008.
  * [9] Marc Beyer and Joachim L Schultze.  Regulatory t cells in cancer.  Blood, 108(3):804–811, 2006.
  * [10] Andrew F Brouwer, Marisa C Eisenberg, and Rafael Meza.  Age effects and temporal trends in hpv-related and hpv-unrelated oral cancer in the united states: a multistage carcinogenesis modeling analysis.  PloS one, 11(3):e0151098, 2016.
  * [11] Gillian Browne, Hanna Taipaleenmäki, Nicole M Bishop, Sharath C Madasu, Leslie M Shaw, Andre J Van Wijnen, Janet L Stein, Gary S Stein, and Jane B Lian.  Runx1 is associated with breast cancer progression in mmtv-pymt transgenic mice and its depletion in vitro inhibits migration and invasion.  Journal of cellular physiology, 230(10):2522–2532, 2015.
  * [12] Aparajita Budithi, Sumeyye Su, Arkadz Kirshtein, and Leili Shahriyari.  Data driven mathematical model of folfiri treatment for colon cancer.  Cancers, 13(11):2632, 2021.
  * [13] Joshua A. Bull and Helen M. Byrne.  Quantification of spatial and phenotypic heterogeneity in an agent-based model of tumour-macrophage interactions.  PLOS Computational Biology, 19(3):1–27, 03 2023.
  * [14] Ying Cai, Ruben Nogales-Cadenas, Quanwei Zhang, Jhih-Rong Lin, Wen Zhang, Kelly O’Brien, Cristina Montagna, and Zhengdong D Zhang.  Transcriptomic dynamics of breast cancer progression in the mmtv-pymt mouse model.  Bmc Genomics, 18(1):1–14, 2017.
  * [15] Lara Campana, Lidia Bosurgi, and Patrizia Rovere-Querini.  Hmgb1: a two-headed signal regulating tumor progression and immunity.  Current opinion in immunology, 20(5):518–523, 2008.
  * [16] Yihai Cao et al.  Angiogenesis modulates adipogenesis and obesity.  The Journal of clinical investigation, 117(9):2362–2368, 2007.
  * [17] Giulia L. Celora, Helen M. Byrne, Christos E. Zois, and P.G. Kevrekidis.  Phenotypic variation modulates the growth dynamics and response to radiotherapy of solid tumours under normoxia and hypoxia.  Journal of Theoretical Biology, 527:110792, 2021.
  * [18] Mark AJ Chaplain.  Mathematical modelling of angiogenesis.  Journal of neuro-oncology, 50:37–51, 2000.
  * [19] Gaoqi Chen, Kaiwen Wu, Hao Li, Demeng Xia, and Tianlin He.  Role of hypoxia in the tumor microenvironment and targeted therapy.  Frontiers in Oncology, 12, 2022.
  * [20] Juan Chen, Yanghui Wei, Weiqin Yang, Qingnan Huang, Yong Chen, Kai Zeng, and Jiawei Chen.  Il-6: the link between inflammation, immunity and breast cancer.  Frontiers in Oncology, 12:903800, 2022.
  * [21] Weizhen Chen, Wei Wei, Liya Yu, Xin Zhang, Fujing Huang, Qiping Zheng, Lingli Wang, and Cheguo Cai.  Baicalin promotes mammary gland development via steroid-like activities.  Frontiers in Cell and Developmental Biology, 9:682469, 2021.
  * [22] Yibing Chen, Yucen Song, Wei Du, Longlong Gong, Haocai Chang, and Zhengzhi Zou.  Tumor-associated macrophages: an accomplice in solid tumor progression.  Journal of biomedical science, 26(1):1–13, 2019.
  * [23] Violeta Chitu, Yee-Guide Yeung, Wenfeng Yu, Sayan Nandi, and E Richard Stanley.  Measurement of macrophage growth and differentiation.  Current protocols in immunology, 92(1):14–20, 2011.
  * [24] Jessica L Christenson, Kiel T Butterfield, Nicole S Spoelstra, John D Norris, Jatinder S Josan, Julie A Pollock, Donald P McDonnell, Benita S Katzenellenbogen, John A Katzenellenbogen, and Jennifer K Richer.  Mmtv-pymt and derived met-1 mouse mammary tumor cells as models for studying the role of the androgen receptor in triple-negative breast cancer progression.  Hormones and Cancer, 8:69–77, 2017.
  * [25] Assunta Cirella, Carlos Luri-Rey, Claudia Augusta Di Trani, Alvaro Teijeira, Irene Olivera, Elixabet Bolaños, Eduardo Castañón, Belen Palencia, Davide Brocco, Myriam Fernández-Sendin, et al.  Novel strategies exploiting interleukin-12 in cancer immunotherapy.  Pharmacology & therapeutics, 239:108189, 2022.
  * [26] Eric T Clambey, Eóin N McNamee, Joseph A Westrich, Louise E Glover, Eric L Campbell, Paul Jedlicka, Edwin F De Zoeten, John C Cambier, Kurt R Stenmark, Sean P Colgan, et al.  Hypoxia-inducible factor-1 alpha–dependent induction of foxp3 drives regulatory t-cell abundance and function during inflammatory hypoxia of the mucosa.  Proceedings of the National Academy of Sciences, 109(41):E2784–E2793, 2012.
  * [27] Oscar R Colegio, Ngoc-Quynh Chu, Alison L Szabo, Thach Chu, Anne Marie Rhebergen, Vikram Jairam, Nika Cyrus, Carolyn E Brokowski, Stephanie C Eisenbarth, Gillian M Phillips, et al.  Functional polarization of tumour-associated macrophages by tumour-derived lactic acid.  Nature, 513(7519):559–563, 2014.
  * [28] Rupert Courtnay, Darleen C Ngo, Neha Malik, Katherine Ververis, Stephanie M Tortorella, and Tom C Karagiannis.  Cancer metabolism and the warburg effect: the role of hif-1 and pi3k.  Molecular biology reports, 42:841–851, 2015.
  * [29] Ralph J DeBerardinis and Navdeep S Chandel.  We need to talk about the warburg effect.  Nature metabolism, 2(2):127–129, 2020.
  * [30] Guo-hu Di, Yang Liu, Ying Lu, Jin Liu, Chutse Wu, and Hai-Feng Duan.  Il-6 secreted from senescent mesenchymal stem cells promotes proliferation and migration of breast cancer cells.  PloS one, 9(11):e113572, 2014.
  * [31] Mirela Domijan, Paul E Brown, Boris V Shulgin, and David A Rand.  Pettsy: a computational tool for perturbation analysis of complex systems biology models.  BMC bioinformatics, 17:1–19, 2016.
  * [32] Haonan Dong, Lu Zhang, and Suling Liu.  Targeting hmgb1: an available therapeutic strategy for breast cancer therapy.  International journal of biological sciences, 18(8):3421, 2022.
  * [33] Aurélie Durgeau, Yasemin Virk, Stéphanie Corgnac, and Fathia Mami-Chouaib.  Recent advances in targeting cd8 t-cell immunity for more effective cancer immunotherapy.  Frontiers in immunology, 9:14, 2018.
  * [34] Peter Ebbesen, Kai-Uwe Eckardt, Fedor Ciampor, and Erik O Pettersen.  Linking measured intercellular oxygen concentration to human cell functions.  Acta Oncologica, 43(6):598–600, 2004.
  * [35] Tarek A El-Mihoub, Adrian A Hopgood, Lars Nolle, and Alan Battersby.  Hybrid genetic algorithms: A review.  Eng. Lett., 13(2):124–137, 2006.
  * [36] Ofer Elhanani, Raz Ben-Uri, and Leeat Keren.  Spatial profiling technologies illuminate the tumor microenvironment.  Cancer cell, 2023.
  * [37] Asieh Emami Nejad, Simin Najafgholian, Alireza Rostami, Alireza Sistani, Samaneh Shojaeifar, Mojgan Esparvarinha, Reza Nedaeinia, Shaghayegh Haghjooy Javanmard, Marjan Taherian, Mojtaba Ahmadlou, et al.  The role of hypoxia in the tumor microenvironment and development of cancer stem cell: a novel approach to developing treatment.  Cancer Cell International, 21(1):1–26, 2021.
  * [38] Napoleone Ferrara.  Vegf and the quest for tumour angiogenesis factors.  Nature Reviews Cancer, 2(10):795–803, 2002.
  * [39] Napoleone Ferrara, Kenneth J Hillan, Hans-Peter Gerber, and William Novotny.  Discovery and development of bevacizumab, an anti-vegf antibody for treating cancer.  Nature reviews Drug discovery, 3(5):391–400, 2004.
  * [40] Irene Filippi, Emilia Morena, Carlo Aldinucci, Fabio Carraro, Silvano Sozzani, and Antonella Naldini.  Short-term hypoxia enhances the migratory capability of dendritic cell through hif-1$\alpha$ and pi3k/akt pathway.  Journal of cellular physiology, 229(12):2067–2076, 2014.
  * [41] Judah Folkman.  Tumor angiogenesis: therapeutic implications.  New england journal of medicine, 285(21):1182–1186, 1971.
  * [42] Stephen B Fox, Daniele G Generali, and Adrian L Harris.  Breast tumour angiogenesis.  Breast cancer research, 9(6):1–11, 2007.
  * [43] Mercedes B Fuertes, Aalok K Kacha, Justin Kline, Seng-Ryong Woo, David M Kranz, Kenneth M Murphy, and Thomas F Gajewski.  Host type i ifn signals are required for antitumor cd8+ t cell responses through cd8$\alpha$\+ dendritic cells.  Journal of Experimental Medicine, 208(10):2005–2016, 2011.
  * [44] Yuhao Gao, Xiaosong Chen, Qing He, Ryan C Gimple, Yajin Liao, Liang Wang, Rong Wu, Qi Xie, Jeremy N Rich, Kunwei Shen, et al.  Adipocytes promote breast tumorigenesis through taz-dependent secretion of resistin.  Proceedings of the National Academy of Sciences, 117(52):33295–33304, 2020.
  * [45] Yael Gropper, Tali Feferman, Tali Shalit, Tomer-Meir Salame, Ziv Porat, and Guy Shakhar.  Culturing ctls under hypoxic conditions enhances their cytolysis and improves their anti-tumor function.  Cell reports, 20(11):2547–2555, 2017.
  * [46] Yuqi Guo, Feng Xu, TianJian Lu, Zhenfeng Duan, and Zhan Zhang.  Interleukin-6 signaling pathway in targeted therapy for cancer.  Cancer treatment reviews, 38(7):904–910, 2012.
  * [47] Carolina E Hagberg, Qian Li, Maria Kutschke, Debajit Bhowmick, Endre Kiss, Irina G Shabalina, Matthew J Harms, Olga Shilkova, Viviana Kozina, Jan Nedergaard, et al.  Flow cytometry of mouse and human adipocytes for the analysis of browning and cellular heterogeneity.  Cell reports, 24(10):2746–2756, 2018.
  * [48] Miranda A Hallett, Bin Teng, Hisashi Hasegawa, Luciana P Schwab, Tiffany N Seagroves, and Tayebeh Pourmotabbed.  Anti-matrix metalloproteinase-9 dnazyme decreases tumor growth in the mmtv-pymt mouse model of breast cancer.  Breast Cancer Research, 15(1):1–11, 2013.
  * [49] Morten Hansen and Mads Hald Andersen.  The role of dendritic cells in cancer.  In Seminars in immunopathology, volume 39, pages 307–316. Springer, 2017.
  * [50] Ning-Bo Hao, Mu-Han Lü, Ya-Han Fan, Ya-Ling Cao, Zhi-Ren Zhang, Shi-Ming Yang, et al.  Macrophages in tumor microenvironments and the progression of tumors.  Journal of Immunology Research, 2012, 2012.
  * [51] Heather A Harrington, M Maier, L Naidoo, Nathaniel Whitaker, and Panayotis G Kevrekidis.  A hybrid model for tumor-induced angiogenesis in the cornea in the presence of inhibitors.  Mathematical and Computer Modelling, 46(3-4):513–524, 2007.
  * [52] Tracy SP Heng, Michio W Painter, the Immunological Genome Project Consortium, Kutlu Elpek, Veronika Lukacs-Kornek, Nora Mauermann, Shannon J Turley, Daphne Koller, Francis S Kim, Amy J Wagers, et al.  The immunological genome project: networks of gene expression in immune cells.  Nature immunology, 9(10):1091–1094, 2008.
  * [53] John H Holland.  Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control, and artificial intelligence.  MIT press, 1992.
  * [54] Martin J Hoogduijn, Johanna C van den Beukel, Lidewij CM Wiersma, and Jooske Ijzer.  Morphology and size of stem cells from mouse and whale: observational study.  Bmj, 347, 2013.
  * [55] Rakesh K Jain.  Normalization of tumor vasculature: an emerging concept in antiangiogenic therapy.  Science, 307(5706):58–62, 2005.
  * [56] Ting Jiang, Xin-Yue Tang, Yi Mao, Yu-Qi Zhou, Jia-Jia Wang, Ruo-Mei Li, Xin-Ru Xie, Hong-Ming Zhang, Bing Fang, Ning-Juan Ouyang, et al.  Matrix mechanics regulate the polarization state of bone marrow-derived neutrophils through the jak1/stat3 signaling pathway.  Acta Biomaterialia, 168:159–173, 2023.
  * [57] Rui Kang, Qiuhong Zhang, Herbert J Zeh III, Michael T Lotze, and Daolin Tang.  Hmgb1 in cancer: good, bad, or both?  Clinical cancer research, 19(15):4046–4057, 2013.
  * [58] Dharambir Kashyap, Deeksha Pal, Riya Sharma, Vivek Kumar Garg, Neelam Goel, Deepika Koundal, Atef Zaguia, Shubham Koundal, Assaye Belay, et al.  Global increase in breast cancer incidence: risk factors and preventive measures.  BioMed research international, 2022, 2022.
  * [59] Kelly Kersten, Karin E de Visser, Martine H van Miltenburg, and Jos Jonkers.  Genetically engineered mouse models in oncology research and cancer medicine.  EMBO molecular medicine, 9(2):137–153, 2017.
  * [60] KL Knutson and ML Disis.  Il-12 enhances the generation of tumour antigen-specific th1 cd4 t cells during ex vivo expansion.  Clinical & Experimental Immunology, 135(2):322–329, 2004.
  * [61] Liselotte M Kornmann, Alma Zernecke, Daniëlle MJ Curfs, Ben JA Janssen, Christian Weber, Menno PJ de Winther, Robert S Reneman, Arnold PG Hoeks, and Koen D Reesink.  Echogenic perfluorohexane-loaded macrophages adhere in vivo to activated vascular endothelium in mice, an explorative study.  Cardiovascular Ultrasound, 13:1–8, 2015.
  * [62] Sophie Laget, Lucile Broncy, Katia Hormigos, Dalia M Dhingra, Fatima BenMohamed, Thierry Capiod, Magne Osteras, Laurent Farinelli, Stephen Jackson, and Patrizia Paterlini-Brechot.  Technical insights into highly sensitive isolation and molecular characterization of fixed and live circulating tumor cells for early detection of tumor invasion.  PloS one, 12(1):e0169427, 2017.
  * [63] Trang Le, Rachel A Aronow, Arkadz Kirshtein, and Leili Shahriyari.  A review of digital cytometry methods: estimating the relative abundance of cell types in a bulk of cells.  Briefings in Bioinformatics, 22(4):bbaa219, 07 2021.
  * [64] Trang Le, Sumeyye Su, and Leili Shahriyari.  Investigating Optimal Chemotherapy Options for Osteosarcoma Patients through a Mathematical Model.  Cells, 10(8):2009, aug 2021.
  * [65] Yoke Seng Lee and Kristen J Radford.  The role of dendritic cells in cancer.  International review of cell and molecular biology, 348:123–178, 2019.
  * [66] Lingling Li, Yulu Hu, Yunshan Xu, and Sanyi Tang.  Mathematical modeling the order of driver gene mutations in colorectal cancer.  PLOS Computational Biology, 19(6):e1011225, 2023.
  * [67] Maria V Liberti and Jason W Locasale.  The warburg effect: how does it benefit cancer cells?  Trends in biochemical sciences, 41(3):211–218, 2016.
  * [68] Juan Liu, Xiaomin Zhang, Kun Chen, Yujie Cheng, Shuxun Liu, Meng Xia, Yali Chen, Ha Zhu, Zhiqing Li, and Xuetao Cao.  Ccr7 chemokine receptor-inducible lnc-dpf3 restrains dendritic cell migration by inhibiting hif-1$\alpha$-mediated glycolysis.  Immunity, 50(3):600–615, 2019.
  * [69] Qiuyu Liu, Victoria AC Palmgren, Erik HJ Danen, and Sylvia E Le Dévédec.  Acute vs. chronic vs. intermittent hypoxia in breast cancer: a review on its application in in vitro research.  Molecular Biology Reports, 49(11):10961–10973, 2022.
  * [70] Zheng-gang Liu and Delong Jiao.  Necroptosis, tumor necrosis and tumorigenesis.  Cell stress, 4(1):1, 2020.
  * [71] Yoram Louzoun, Chuan Xue, Gregory B Lesinski, and Avner Friedman.  A mathematical model for pancreatic cancer growth and treatments.  Journal of theoretical biology, 351:74–82, 2014.
  * [72] Danny Luan, Benyamin Dadpey, Jessica Zaid, Pania E Bridge-Comer, Julia H DeLuca, Wenmin Xia, Joshua Castle, and Shannon M Reilly.  Adipocyte-secreted il-6 sensitizes macrophages to il-4 signaling.  Diabetes, 72(3):367–374, 2023.
  * [73] Xiaojing Ma, Wenjun Yan, Hua Zheng, Qinglin Du, Lixing Zhang, Yi Ban, Na Li, and Fang Wei.  Regulation of il-10 and il-12 production and function in macrophages and dendritic cells.  F1000Research, 4, 2015.
  * [74] Grace Mahlbacher, Louis T Curtis, John Lowengrub, and Hermann B Frieboes.  Mathematical modeling of tumor-associated macrophage interactions with the cancer microenvironment.  Journal for immunotherapy of cancer, 6(1):1–17, 2018.
  * [75] Ahmed M Makhlouf, Lamiaa El-Shennawy, Hesham A Elkaranshawy, et al.  Mathematical modelling for the role of cd4+ t cells in tumor-immune interactions.  Computational and mathematical methods in medicine, 2020, 2020.
  * [76] Ariel E Marciscano and Niroshana Anandasabapathy.  The role of dendritic cells in cancer and anti-tumor immunity.  In Seminars in immunology, volume 52, page 101481. Elsevier, 2021.
  * [77] Gemma Massonis, Julio R Banga, and Alejandro F Villaverde.  Autorepar: a method to obtain identifiable and observable reparameterizations of dynamic models with mechanistic insights.  International Journal of Robust and Nonlinear Control, 33(9):5039–5057, 2023.
  * [78] Gemma Massonis and Alejandro F Villaverde.  Finding and breaking lie symmetries: implications for structural identifiability and observability in biological modelling.  Symmetry, 12(3):469, 2020.
  * [79] Manuela Milani and Adrian L Harris.  Targeting tumour hypoxia in breast cancer.  European journal of cancer, 44(18):2766–2773, 2008.
  * [80] Kathy Miller, Molin Wang, Julie Gralow, Maura Dickler, Melody Cobleigh, Edith A Perez, Tamara Shenkier, David Cella, and Nancy E Davidson.  Paclitaxel plus bevacizumab versus paclitaxel alone for metastatic breast cancer.  New England journal of medicine, 357(26):2666–2676, 2007.
  * [81] Bhalchandra Mirlekar.  Tumor promoting roles of il-10, tgf-$\beta$, il-4, and il-35: Its implications in cancer immunotherapy.  SAGE open medicine, 10:20503121211069012, 2022.
  * [82] Bhalchandra Mirlekar and Yuliya Pylayeva-Gupta.  Il-12 family cytokines in cancer and immunotherapy.  Cancers, 13(2):167, 2021.
  * [83] Navid Mohammad Mirzaei, Wenrui Hao, and Leili Shahriyari.  Investigating the spatial interaction of immune cells in colon cancer.  iScience, 26:106596, 5 2023.
  * [84] Navid Mohammad Mirzaei and Leili Shahriyari.  Modeling cancer progression: An integrated workflow extending data-driven kinetic models to bio-mechanical pde models.  arXiv preprint arXiv:2311.02232, 2023.
  * [85] Navid Mohammad Mirzaei, Navid Changizi, Alireza Asadpoure, Sumeyye Su, Dilruba Sofia, Zuzana Tatarova, Ioannis K Zervantonakis, Young Hwan Chang, and Leili Shahriyari.  Investigating key cell types and molecules dynamics in pymt mice model of breast cancer through a mathematical model.  PLoS computational biology, 18(3):e1009953, 2022.
  * [86] Navid Mohammad Mirzaei, Sumeyye Su, Dilruba Sofia, Maura Hegarty, Mohamed H Abdel-Rahman, Alireza Asadpoure, Colleen M Cebulla, Young Hwan Chang, Wenrui Hao, Pamela R Jackson, et al.  A mathematical model of breast tumor progression based on immune infiltration.  Journal of Personalized Medicine, 11(10):1031, 2021.
  * [87] Navid Mohammad Mirzaei, Zuzana Tatarova, Wenrui Hao, Navid Changizi, Alireza Asadpoure, Ioannis K Zervantonakis, Yu Hu, Young Hwan Chang, and Leili Shahriyari.  A pde model of breast tumor progression in mmtv-pymt mice.  Journal of Personalized Medicine, 12(5):807, 2022.
  * [88] Severin Mühleder, Macarena Fernández-Chacón, Irene Garcia-Gonzalez, and Rui Benedito.  Endothelial sprouting, proliferation, or senescence: tipping the balance from physiology to pathology.  Cellular and Molecular Life Sciences, 78(4):1329–1354, 2021.
  * [89] Anna-katharina Müller, Ulrike A Köhler, Sébastien Trzebanski, Yaron Vinik, Harsha Mohan Raj, Jean-antoine Girault, Nir Ben-Chetrit, Antonio Maraver, Steffen Jung, and Sima Lev.  Mouse modeling dissecting macrophage–breast cancer communication uncovered roles of pyk2 in macrophage recruitment and breast tumorigenesis.  Advanced Science, 9(9):2105696, 2022.
  * [90] Barbara Muz, Pilar de la Puente, Feda Azab, and Abdel Kareem Azab.  The role of hypoxia in cancer progression, angiogenesis, metastasis, and resistance to therapy.  Hypoxia, 3:83, 2015.
  * [91] Frank Nani and Herb I Freedman.  A mathematical model of cancer treatment by immunotherapy.  Mathematical biosciences, 163(2):159–199, 2000.
  * [92] Aaron M Newman, Chloé B Steen, Chih Long Liu, Andrew J Gentles, Aadel A Chaudhuri, Florian Scherer, Michael S Khodadoust, Mohammad S Esfahani, Bogdan A Luca, David Steiner, et al.  Determining cell type abundance and expression from bulk tissues with digital cytometry.  Nature biotechnology, 37(7):773–782, 2019.
  * [93] Kristin M Nieman, Iris L Romero, Bennett Van Houten, and Ernst Lengyel.  Adipose tissue and adipocytes support tumorigenesis and metastasis.  Biochimica et Biophysica Acta (BBA)-Molecular and Cell Biology of Lipids, 1831(10):1533–1541, 2013.
  * [94] Muhammad Zaeem Noman, Meriem Hasmim, Yosra Messai, Stéphane Terry, Claudine Kieda, Bassam Janji, and Salem Chouaib.  Hypoxia: a key player in antitumor immune response. a review in the theme: cellular responses to hypoxia.  American Journal of Physiology-Cell Physiology, 309(9):C569–C579, 2015.
  * [95] Pia Nyberg, Liang Xie, and Raghu Kalluri.  Endogenous inhibitors of angiogenesis.  Cancer research, 65(10):3967–3979, 2005.
  * [96] Karyn E O’Connell, Amy M Mikkola, Aaron M Stepanek, Andyna Vernet, Christopher D Hall, Chia C Sun, Eda Yildirim, John F Staropoli, Jeannie T Lee, and Diane E Brown.  Practical murine hematopathology: a comparative review and implications for research.  Comparative medicine, 65(2):96–113, 2015.
  * [97] Barbara Onnis, Annamaria Rapisarda, and Giovanni Melillo.  Development of hif-1 inhibitors for cancer therapy.  Journal of cellular and molecular medicine, 13(9a):2780–2786, 2009.
  * [98] Karolina Palucka, Lisa M Coussens, and Joyce O’Shaughnessy.  Dendritic cells, inflammation and breast cancer.  Cancer Journal (Sudbury, Mass.), 19(6), 2013.
  * [99] Shayn M. Peirce.  Computational and mathematical modeling of angiogenesis.  Microcirculation, 15(8):739–751, 2008.
  * [100] Anand Pillarisetti, Hamid Ladjal, Antoine Ferreira, Carol Keefer, and Jaydev P Desai.  Mechanical characterization of mouse embryonic stem cells.  In 2009 Annual International Conference of the IEEE Engineering in Medicine and Biology Society, pages 1176–1179. IEEE, 2009.
  * [101] Hans Raskov, Adile Orhan, Jan Pravsgaard Christensen, and Ismail Gögenur.  Cytotoxic cd8+ t cells in cancer and cancer immunotherapy.  British journal of cancer, 124(2):359–367, 2021.
  * [102] Andreas Raue, Clemens Kreutz, Thomas Maiwald, Julie Bachmann, Marcel Schilling, Ursula Klingmüller, and Jens Timmer.  Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood.  Bioinformatics, 25(15):1923–1929, 2009.
  * [103] Andrii Rozhok and James DeGregori.  A generalized theory of age-dependent carcinogenesis.  Elife, 8:e39950, 2019.
  * [104] Mohamad Adham Salkeni and Aung Naing.  Interleukin-10 in cancer immunotherapy: from bench to bedside.  Trends in Cancer, 2023.
  * [105] Takami Sato, Mizue Terai, Yutaka Tamura, Vitali Alexeev, Michael J Mastrangelo, and Senthamil R Selvan.  Interleukin 10 in the tumor microenvironment: a target for anticancer immunotherapy.  Immunologic research, 51:170–182, 2011.
  * [106] Nicholas Schaum, Jim Karkanias, Norma F Neff, Andrew P May, Stephen R Quake, Tony Wyss-Coray, Spyros Darmanis, Joshua Batson, Olga Botvinnik, Michelle B Chen, et al.  Single-cell transcriptomics of 20 mouse organs creates a tabula muris: The tabula muris consortium.  Nature, 562(7727):367, 2018.
  * [107] Gregg L Semenza.  Targeting hif-1 for cancer therapy.  Nature reviews cancer, 3(10):721–732, 2003.
  * [108] Ganesh N Sharma, Rahul Dave, Jyotsana Sanadya, Piush Sharma, and KK22247839 Sharma.  Various types and management of breast cancer: an overview.  Journal of advanced pharmaceutical technology & research, 1(2):109, 2010.
  * [109] Norman E Sharpless and Ronald A DePinho.  The mighty mouse: genetically engineered mouse models in cancer drug development.  Nature reviews Drug discovery, 5(9):741–754, 2006.
  * [110] Elnaz Sheikhpour, Parisa Noorbakhsh, Elnaz Foroughi, Soudabeh Farahnak, Rezvan Nasiri, and Hossein Neamatzadeh.  A survey on the role of interleukin-10 in breast cancer: A narrative.  Reports of biochemistry & molecular biology, 7(1):30, 2018.
  * [111] Rebecca L Siegel, Kimberly D Miller, Nikita Sandeep Wagle, and Ahmedin Jemal.  Cancer statistics, 2023.  Ca Cancer J Clin, 73(1):17–48, 2023.
  * [112] Ludivine Soguel, Francine Durocher, André Tchernof, and Caroline Diorio.  Adiposity, breast density, and breast cancer risk: epidemiological and biological considerations.  European journal of cancer prevention, 26(6):511, 2017.
  * [113] Savas D Soysal, Alexandar Tzankov, and Simone E Muenst.  Role of the tumor microenvironment in breast cancer.  Pathobiology, 82(3-4):142–152, 2015.
  * [114] Catherine E Steding, Sung-tse Wu, Yanping Zhang, Meei-Huey Jeng, Bennett D Elzey, and Chinghai Kao.  The role of interleukin-12 on modulating myeloid-derived suppressor cells, increasing overall survival and reducing metastasis.  Immunology, 133(2):221–238, 2011.
  * [115] Ralph M Steinman and Jacques Banchereau.  Taking dendritic cells into medicine.  Nature, 449(7161):419–426, 2007.
  * [116] Lijuan Sun, Xiao Han, and Mikala Egeblad.  Isolation of mouse mammary carcinoma-derived macrophages and cancer cells for co-culture assays.  STAR protocols, 3(4):101833, 2022.
  * [117] Lijuan Sun, Tim Kees, Ana Santos Almeida, Bodu Liu, Xue-Yan He, David Ng, Xiao Han, David L Spector, Iain A McNeish, Phyllis Gimotty, et al.  Activating a collaborative innate-adaptive immune response to control metastasis.  Cancer cell, 39(10):1361–1374, 2021.
  * [118] Humayra Tasnim, G Matthew Fricke, Janie R Byrum, Justyna O Sotiris, Judy L Cannon, and Melanie E Moses.  Quantitative measurement of naïve t cell association with dendritic cells, frcs, and blood vessels in lymph nodes.  Frontiers in immunology, 9:1571, 2018.
  * [119] Zuzana Tatarova, Dylan C Blumberg, James E Korkola, Laura M Heiser, John L Muschler, Pepper J Schedin, Sebastian W Ahn, Gordon B Mills, Lisa M Coussens, Oliver Jonas, et al.  A multiplex implantable microdevice assay identifies synergistic combinations of cancer immunotherapies and conventional drugs.  Nature biotechnology, 40(12):1823–1833, 2022.
  * [120] Sheng Tong and Fan Yuan.  Numerical simulations of angiogenesis in the cornea.  Microvascular Research, 61(1):14–27, 2001.
  * [121] Jo Marie Tran Janco, Purushottam Lamichhane, Lavakumar Karyampudi, and Keith L Knutson.  Tumor-infiltrating dendritic cells in cancer pathogenesis.  The Journal of Immunology, 194(7):2985–2991, 2015.
  * [122] Paul Trayhurn.  Hypoxia and adipose tissue function and dysfunction in obesity.  Physiological reviews, 93(1):1–21, 2013.
  * [123] Alok Tripathi, Kriti Shrinet, and Arvind Kumar.  Hmgb1 protein as a novel target for cancer.  Toxicology reports, 6:253–261, 2019.
  * [124] Marta Truffi, Luca Sorrentino, and Fabio Corsi.  Fibroblasts in the tumor microenvironment.  Tumor Microenvironment: Non-Hematopoietic Cells, pages 15–29, 2020.
  * [125] Yagmur Turgay, Matthias Eibauer, Anne E Goldman, Takeshi Shimi, Maayan Khayat, Kfir Ben-Harush, Anna Dubrovsky-Gaupp, K Tanuj Sapra, Robert D Goldman, and Ohad Medalia.  The molecular architecture of lamins in somatic cells.  Nature, 543(7644):261–264, 2017.
  * [126] Peter Vaupel, Michael Höckel, and Arnulf Mayer.  Detection and characterization of tumor hypoxia using po2 histography.  Biophysical Journal, 121(19), 2007.
  * [127] Alejandro F Villaverde, Antonio Barreiro, and Antonis Papachristodoulou.  Structural identifiability of dynamic systems biology models.  PLoS computational biology, 12(10):e1005153, 2016.
  * [128] Alejandro F Villaverde, Neil D Evans, Michael J Chappell, and Julio R Banga.  Input-dependent structural identifiability of nonlinear systems.  IEEE Control Systems Letters, 3(2):272–277, 2018.
  * [129] Alejandro F Villaverde, Nikolaos Tsiantis, and Julio R Banga.  Full observability and estimation of unknown inputs, states and parameters of nonlinear biological models.  Journal of the Royal Society Interface, 16(156):20190043, 2019.
  * [130] Linda Vona-Davis and David P Rose.  Angiogenesis, adipokines and breast cancer.  Cytokine & growth factor reviews, 20(3):193–201, 2009.
  * [131] Jessica C Walrath, Jessica J Hawes, Terry Van Dyke, and Karlyne M Reilly.  Genetically engineered mouse models in cancer research.  Advances in cancer research, 106:113–164, 2010.
  * [132] Shihua Wang, Xiaodong Su, Meiqian Xu, Xian Xiao, Xiaoxia Li, Hongling Li, Armand Keating, and Robert Chunhua Zhao.  Exosomes secreted by mesenchymal stromal/stem cell-derived adipocytes promote breast cancer cell growth via activation of hippo signaling pathway.  Stem Cell Research & Therapy, 10:1–12, 2019.
  * [133] Maria Angelica Ehara Watanabe, Julie Massayo Maeda Oda, Marla Karine Amarante, and Julio Cesar Voltarelli.  Regulatory t cells and breast cancer: implications for immunopathogenesis.  Cancer and Metastasis Reviews, 29:569–579, 2010.
  * [134] Jeffrey West and Paul K Newton.  Cellular interactions constrain tumor growth.  Proceedings of the National Academy of Sciences, 116(6):1918–1923, 2019.
  * [135] Theresa L Whiteside.  What are regulatory t cells (treg) regulating in cancer and why?  In Seminars in cancer biology, volume 22, pages 327–334. Elsevier, 2012.
  * [136] Franz-Georg Wieland, Adrian L Hauber, Marcus Rosenblatt, Christian Tönsing, and Jens Timmer.  On structural and practical identifiability.  Current Opinion in Systems Biology, 25:60–69, 2021.
  * [137] Ami Yamamoto, Yin Huang, Brad A Krajina, Margaux McBirney, Andrea E Doak, Sixuan Qu, Carolyn L Wang, Michael C Haffner, and Kevin J Cheung.  Metastasis from the tumor interior and necrotic core formation are regulated by breast cancer-derived angiopoietin-like 7.  Proceedings of the National Academy of Sciences, 120(10):e2214888120, 2023.
  * [138] Ying Zhang, Yong-ming Yao, Li-feng Huang, Ning Dong, Yan Yu, and Zhi-yong Sheng.  The potential effect and mechanism of high-mobility group box 1 protein on regulatory t cell-mediated immunosuppression.  Journal of Interferon & Cytokine Research, 31(2):249–257, 2011.