Data Driven Mathematical Model of Colon Cancer Progression

Abstract

Every colon cancer has its own unique characteristics, and therefore may respond differently to identical treatments. Here, we develop a data driven mathematical model for the interaction network of key components of immune microenvironment in colon cancer. We estimate the relative abundance of each immune cell from gene expression profiles of tumors, and group patients based on their immune patterns. Then we compare the tumor sensitivity and progression in each of these groups of patients, and observe differences in the patterns of tumor growth between the groups. For instance, in tumors with a smaller density of naive macrophages than activated macrophages, a higher activation rate of macrophages leads to an increase in cancer cell density, demonstrating a negative effect of macrophages. Other tumors however, exhibit an opposite trend, showing a positive effect of macrophages in controlling tumor size. Although the results indicate that for all patients the size of the tumor is sensitive to the parameters related to macrophages, such as their activation and death rate, this research demonstrates that no single biomarker could predict the dynamics of tumors.

1. Introduction

Recent studies show that many cancers arise from sites of chronic inflammation [1,2,3,4]. Balkwill et al. [5] provide a list of inflammatory conditions that predispose an individual to cancer, in particular to colorectal cancer. Indeed, inflammatory bowel diseases like ulcerative colitis and colonic Crohn’s disease are strongly associated with colorectal cancer [6]. For example, inducing colitis to create chronic inflammation after introducing procarcinogen is an established and reliable two-step mouse model of colitis-associated cancer (CAC) [7,8,9].

Most common cancer treatments are designed to kill tumor cells. However, the way in which cells die is very important, because dying cells may release molecules that initiate an immune response. We shall refer to cells that go through the process of necrotic cell death as necrotic cells. Necrotic cells are known to release damage-associated molecular pattern (DAMP) molecules such as high mobility group box 1 (HMGB1), which triggers immune responses [10,11]. In particular, HMGB1 activates dendritic cells (DCs) [12]. There is evidence that the expressions of HMGB1 and RAGE, its receptor, are significantly higher in ulcerative colitis than in control cases [13]. HMGB1 has been observed in other cancers, as a result of treatments by radiotherapy and chemotherapy [12,14,15,16].

In colon cancer, activated CD8⁺ T cells enhance production of necrotic cells by expressing high levels of cytokines like IFN-$\gamma$ and FasL [17]. Necrotic cells and macrophages release HMGB1 to activate dendritic cells [12], which leads to activation of T-cells [18]. In addition, intestinal epithelial cells, which are in close contact with DCs, activate dendritic cells by releasing molecules like thymic stromal lymphopoietin (TSLP) [19,20]. Once activated, dendritic cells release cytokines STAT4, STAT6, and IL-4, which induce differentiation of naive T-cells into effector T cells (Th1, Th17 and Th2) [21]. CD4${}^{+}$ T-cells can also become activated by TNF-$\alpha$, which is released by $M_1$ macrophages [22]. Activated CD4${}^{+}$ T-cells release IL-2, 4, 5, 13 and 17 to activate killer cells like CD8${}^{+}$ T-cells [18,23,24]. CD4${}^{+}$ T-cells also release IFN-$\gamma$, which activates $M_1$ macrophages [25,26]. Activated macrophages and CD4${}^{+}$ effector T-cells release tumor-promoting cytokines interleukin 6 (IL-6) [27]. IL-6 promotes tumor growth by activating STAT3 in intestinal epithelial cells [28].

Knowledge of the cancer microenvironment is essential in predicting the progression of cancer. A strong correlation between in situ immune reactions in tumor regions and prognosis has been observed regardless of the local extent of the tumor and of invasion of regional lymph nodes [29]. A weak in situ immune reaction in tumor regions is associated with a poor prognosis even in patients with minimal tumor invasion (stage I). Moreover, high expression of the Th17 markers predicts a poor prognosis for patients with colorectal cancer, whereas patients with high expression of the Th1 markers have prolonged disease-free survival [30]. However, it has been observed that a high proportion of CD8${}^{+}$ T cells, effector memory T cells and CD4${}^{+}$ T cells is correlated with longer survival in colorectal cancer [31,32,33]. Moreover, in colon cancer, patients with a low level of macrophages have a deeper depth of invasion than patients with a high level of macrophages [33]. All these observations indicate the importance of the relative abundance of various immune cells, as well as their interaction networks, in the colonic tumors’ initiation and progression. Therefore, to accurately model the progress of cancer, we need to divide patients into similar cohorts based on their tumor-infiltrating immune cells and predict the progression for each group separately.

While there are many papers that use mathematical models for colon cancer progression [34,35,36,37,38,39,40,41,42,43], only a few have attempted to include immune interaction in their model. Models such as [40,41,42] define a system of ordinary differential equations (ODEs) that describe the interactions between cancerous cells and various sub-populations of immune cells (including NK cells, CD8+ T cells, lymphocytes, natural death cells and interleukins) and explore how these interactions can influence tumor growth over time. While time course data for the growth of untreated tumors are not currently easily available to verify models such as [40], other models such as [41] include simulations of treatment plans that can be compared on the population-level to results from previous clinical trials. To generate population-level simulation results while acknowledging the different responses to treatment that can arise from differing patient immune profiles, this study selects a range of parameter values to simulate 64 unique “virtual patients” for which to solve the system of ODEs describing potential treatments.

In the present paper, we develop a data driven mathematical model of colon cancer with emphasis on the role of immune cells, including T-cells, dendritic cells and macrophages. Although there are many cell types and molecules involved in colon cancer, in order to avoid too much complexity, we only model some of the key players and interaction networks that we determined by reviewing articles on colon cancer progression. The resulting mathematical model is based on the network shown in Figure 1, and it is represented by a system of ODEs within the tumor. In order to explore differences in tumor growth among patients with different immune profiles, we use cancer patients’ data to estimate the percentage of each immune cell type in their primary tumors, and use clustering to group these patients into five distinct groups of immune patterns. We then use the data within each cluster to generate five “virtual patients”, for which we can calculate patient-specific parameters to use in the mathematical model. Lastly, we examine the differences in the resulting dynamics between the 5 clusters, and look for potential biomarkers that can link details of the tumor microenvironment to the dynamics of tumor growth.

2. Materials and Methods

2.1. Mathematical Model

We develop a mathematical model for colon cancer based on the interaction network among key players in colon cancer shown in Figure 1, and the list of variables is given in

Table 1. Model’s variables. Names and descriptions of variables used in the model.

Variable	Name	Description
$T_N$	Naive T-cells
$T_h$	Helper T-cells
$T_C$	Cytotoxic cells	includes CD8+ T-cells and NK cells
$T_r$	Regulatory T-cells
$D_n$	Naive dendritic cells
D	Activated dendritic cells	antigen presenting cells
M	Macrophages
C	Cancer cells
N	Nectrotic cells
H	HMGB1
${\mu }_1$	Carcinogenic cytokines	includes effects of IL-6, IL-17, IL-21 and IL-22
${\mu }_2$	Immunosuppresive agents	includes effects of IL-10 and CCL20
$I_{\gamma }$	IFN-$\gamma$
$G_{\beta }$	TGF-$\beta$

. The model is represented by a system of differential equations for concentrations and changing in time in unit of day. For clarity, we develop a simplified model in terms of ordinary differential equations. For biochemical processes $A+B\to C$, we use the mass action law $\frac{dC}{dt}=\lambda AB$, where $\lambda$ is production rate of C [44,45]. Throughout the paper, we use the symbol $\lambda$ for production, activation or proliferation rates, and the symbol $\delta$ for decay, natural death or premature death (necrosis) rates.

2.1.1. Cytokine Approximation

In order to reduce the complexity of the system, we treat some of cytokines as independent variables and approximate the value of other cytokines through already existing variables. Additionally, we combine the cytokines that have a similar function in the interaction network (Figure 1). Therefore, we combine IL-6, IL-17, IL-21 and IL-22 and denote their sum by the variable ${\mu }_1$. We also combine IL-10 and CCL20 and denote their sum by the variable ${\mu }_2$. The cytokines treated as model variables are HMGB1, IFN-$\gamma$, TGF-$\beta$, IL-6 and IL-10. We then model the dynamics of cytokines in the following way.

HMGB1 is passively released from necrotic cells [46], or actively secreted from activated T-cells and macrophages [47,48]. Thus, we can model the dynamics of HMGB1 by the equation:

$$\frac{d\left[H\right]}{dt}={\lambda }_{HN}\left[N\right]+{\lambda }_{HM}\left[M\right]+{\lambda }_{H{T}_h}\left[T_h\right]+{\lambda }_{H{T}_C}\left[T_C\right]+{\lambda }_{H{T}_r}\left[T_r\right]-{\delta }_H\left[H\right].$$

(1)

IL-6 is secreted by tumor associated macrophages (TAMs) [27,49,50,51], helper T-cells [27,51,52,53] and sub-population of dendritic cells [54,55]. IL-17, IL-21 and IL-22 are produced by helper T-cells [56]. Therefore, the resulting dynamics for $\left[{\mu }_1\right]$ can be written as

$$\frac{d\left[{\mu }_1\right]}{dt}={\lambda }_{{\mu }_1{T}_h}\left[T_h\right]+{\lambda }_{{\mu }_{1}M}\left[M\right]+{\lambda }_{{\mu }_{1}D}\left[D\right]-{\delta }_{{\mu }_1}\left[{\mu }_1\right].$$

(2)

IL-10 is produced by macrophages [57,58], dendritic cells [54,59] and T-reg cells [52,56,60,61]. CCL20 is produced by macrophages [62]. Thus, the equation for $\left[{\mu }_2\right]$ is

$$\frac{d\left[{\mu }_2\right]}{dt}={\lambda }_{{\mu }_{2}M}\left[M\right]+{\lambda }_{{\mu }_{2}D}\left[D\right]+{\lambda }_{{\mu }_2{T}_r}\left[T_r\right]-{\delta }_{{\mu }_2}\left[{\mu }_2\right].$$

(3)

IFN-$\gamma$ is secreted by a sub-population of macrophages [57,63,64,65,66], helper T-cells [25,26] and cytotoxic cells [17], which results in the following equation:

$$\frac{d\left[I_{\gamma }\right]}{dt}={\lambda }_{I_{\gamma }{T}_h}\left[T_h\right]+{\lambda }_{I_{\gamma }{T}_C}\left[T_C\right]+{\lambda }_{I_{\gamma }M}\left[M\right]-{\delta }_{I_{\gamma }}\left[I_{\gamma }\right].$$

(4)

TGF-$\beta$ is produced by macrophages [57,58] and T-reg cells [52,56,60,67] leading to the equation:

$$\frac{d\left[G_{\beta }\right]}{dt}={\lambda }_{G_{\beta }M}\left[M\right]+{\lambda }_{G_{\beta }{T}_r}\left[T_r\right]-{\delta }_{G_{\beta }}\left[G_{\beta }\right].$$

(5)

Other cytokines, like IL-2, IL-4, IL-5 and IL-13, we consider to be in a quasi-equilibrium state, i.e., proportional to the concentration of cells that secrete/produce them. In particular, IL-2, IL-5 and IL-13 are produced by CD4+ T-cells [18,24,68], so we consider

$$\begin{array}{c}\hfill \left[IL\text{-}2\right]\approx \mathrm{Const}\times \left[T_h\right],\\ \hfill \left[IL\text{-}5\right]\approx \mathrm{Const}\times \left[T_h\right],\\ \hfill \left[IL\text{-}13\right]\approx \mathrm{Const}\times \left[T_h\right].\end{array}$$

IL-4 is also produced both by CD4+ T-cells [18,24,68] and dendritic cells [21], so we take

$$\left[IL\text{-}4\right]\approx \mathrm{Const}\times \left[T_h\right]+\mathrm{Const}\times \left[D\right].$$

IL-12 secreted by macrophages [57,58] and dendritic cells [49,50,54,56,69,70], thus can be approximated as

$$\left[IL\text{-}12\right]\approx \mathrm{Const}\left[M\right]+\mathrm{Const}\left[D\right],$$

while IL-23 and TNF-$\alpha$ are secreted solely by macrophages [49,50]; hence, their approximation is

$$\left[IL\text{-}23\right]\approx \mathrm{Const}\left[M\right],\left[TNF\text{-}\alpha \right]\approx \mathrm{Const}\left[M\right].$$

2.1.2. T-Cells

In this model, we differentiate four subgroups of T-cells: naive, helper, cytotoxic and regulatory.

Naive T-cells, $T_N$, are not necessarily part of tumor microenvironment, as they usually are activated within lymph nodes. However, making activation rates for other types of T-cells proportional to the density of naive cells creates a better controlled system and avoids unlimited exponential growth. Thus, we summarize the equation for the dynamics of the naive T-cells after detailing the equations of other types of T-cells. The variables $T_h$, $T_C$ and $T_r$ correspond to the concentration of activated T-helper, cytotoxic and T-reg cells, respectively.

Helper T-cells can be activated with antigen presentation by dendritic cells [18]. CD4+ T-cells can be additionally activated by IL-12, while Th17 are activated by IL-6, TNF-$\alpha$ and IL-23 [56]. Regulatory T-cells inhibit protective immune response (helper and cytotoxic T-cells) in several ways including production of immunosuppresive cytokines such as IL-10 and CCL20 as well as through contact-dependent mechanisms [56]. Additionally, we introduce the natural death rate for helper cells ${\delta }_{T_h}$. The resulting equation is

$$\frac{d\left[T_h\right]}{dt}=\left({\lambda }_{T_{h}D}\left[D\right]+{\lambda }_{T_{h}M}\left[M\right]+{\lambda }_{T_h{\mu }_1}\left[{\mu }_1\right]\right)\left[T_N\right]-\left({\delta }_{T_h{\mu }_2}\left[{\mu }_2\right]+{\delta }_{T_h{T}_r}\left[T_r\right]+{\delta }_{T_h}\right)\left[T_h\right].$$

(6)

The variable corresponding to cytotoxic cells accounts for the effects of cytotoxic T-lymphocytes (mainly CD8+ T-cells) and possibly natural killer cells. CD8+ T-cells are activated by IL-2, IL-4, IL-5 and IL-13 [18,24,68]. Cumulative effect of these cytokines can be written as

$$\left[IL\text{-}2,4,5,13\right]\approx \mathrm{Const}\times \left[T_h\right]+\mathrm{Const}\times \left[D\right].$$

Activation of natural killer cells requires IL-2 [71], which is already included. We also include inhibitory effects mediated by T-reg cells. The dynamics of $T_C$ cell group is modeled by the following equation:

$$\frac{d\left[T_C\right]}{dt}=\left({\lambda }_{T_C{T}_h}\left[T_h\right]+{\lambda }_{T_{C}D}\left[D\right]\right)\left[T_N\right]-\left({\delta }_{T_C{\mu }_2}\left[{\mu }_2\right]+{\delta }_{T_C{T}_r}\left[T_r\right]+{\delta }_{T_C}\right)\left[T_C\right].$$

(7)

Regulatory T-cells can be activated by IL-2 [56,72,73], CCL20 [62] and TGF-$\beta$ [56,67]. IL-6 suppresses T-reg differentiation and shifts it towards T-helper type [74]. The resulting dynamics can be described as follow:

$$\frac{d\left[T_r\right]}{dt}=\left({\lambda }_{T_r{T}_h}\left[T_h\right]+{\lambda }_{T_r{\mu }_2}\left[{\mu }_2\right]+{\lambda }_{T_r{G}_{\beta }}\left[G_{\beta }\right]\right)\left[T_N\right]-\left({\delta }_{T_r{\mu }_1}\left[{\mu }_1\right]+{\delta }_{T_r}\right)\left[T_r\right].$$

(8)

Combining all activation and introducing independent naive T-cell production $A_{T_N}$, we get the following equation for naive T-cells:

$$\begin{array}{cc}\hfill \frac{d\left[T_N\right]}{dt}=& A_{T_N}-\left({\lambda }_{T_{h}D}\left[D\right]+{\lambda }_{T_{h}M}\left[M\right]+{\lambda }_{T_h{\mu }_1}\left[{\mu }_1\right]\right)\left[T_N\right]\hfill \\ \hfill & -\left({\lambda }_{T_C{T}_h}\left[T_h\right]+{\lambda }_{T_{C}D}\left[D\right]\right)\left[T_N\right]\hfill \\ \hfill & -\left({\lambda }_{T_r{T}_h}\left[T_h\right]+{\lambda }_{T_r{\mu }_2}\left[{\mu }_2\right]+{\lambda }_{T_r{G}_{\beta }}\left[G_{\beta }\right]\right)\left[T_N\right]\hfill \\ \hfill & -{\delta }_{T_N}\left[T_N\right].\hfill \end{array}$$

(9)

2.1.3. Dendritic Cells

Dendritic cells become activated by HMGB1 [12]. Moreover, TSLP, which is released by cancer cells [19,20], leads to the activation of dendritic cells (green arrow in Figure 1 from cancer cells to DCs). We take TSLP in quasi-equilibrium state as

$$\left[TSLP\right]\approx \mathrm{Const}\times \left[C\right].$$

On the other hand, multiple factors induced by cancer cells may promote natural death of dendritic cells [75,76,77,78,79] (black arrow in Figure 1 from cancer cells to DCs). Additionally, there’s evidence that HMGB1 can reduce the maturation rate of dendritic cells [61,79]. Introducing the independent production rate of naive dendritic cells $A_{D_N}$, we get the following system for dynamics of naive ($D_N$) and activated (D) dendritic cells:

$$\begin{array}{cc}\hfill \frac{d\left[D_N\right]}{dt}=& A_{D_N}-\left({\lambda }_{DH}\left[H\right]+{\lambda }_{DC}\left[C\right]\right)\left[D_N\right]-\left({\delta }_{DH}\left[H\right]+{\delta }_D\right)\left[D_N\right],\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \frac{d\left[D\right]}{dt}=& \left({\lambda }_{DH}\left[H\right]+{\lambda }_{DC}\left[C\right]\right)\left[D_N\right]-\left({\delta }_{DH}\left[H\right]+{\delta }_{DC}\left[C\right]+{\delta }_D\right)\left[D\right].\hfill \end{array}$$

(11)

2.1.4. Macrophages

There are two main sub-types of macrophages: M1 and M2. M1 phenotype can be activated by IFN-$\gamma$, while M2 can be activated IL-4 and IL-13, which are secreted by helper T-cells [57,58]. Additionally there’s a possibility of TAM activation by IL-10 [57,80,81]. Introducing naive ($M_N$) and activated (M) TAMs, as well as production rate for naive macrophages $A_M$, we can write the following system:

$$\begin{array}{cc}\hfill \frac{d\left[M_N\right]}{dt}=& A_M-\left({\lambda }_{M{\mu }_2}\left[{\mu }_2\right]+{\lambda }_{M{I}_{\gamma }}\left[I_{\gamma }\right]+{\lambda }_{M{T}_h}\left[T_h\right]\right)\left[M_N\right]-{\delta }_M\left[M_N\right],\hfill \\ \hfill \frac{d\left[M\right]}{dt}=& \left({\lambda }_{M{\mu }_2}\left[{\mu }_2\right]+{\lambda }_{M{I}_{\gamma }}\left[I_{\gamma }\right]+{\lambda }_{M{T}_h}\left[T_h\right]\right)\left[M_N\right]-{\delta }_M\left[M\right].\hfill \end{array}$$

It has been shown that there is a diverse spectrum of TAM sub-types [82]. We therefore decided to combine all activated macrophages, including both M1 and M2 macrophages into one variable M and let the system parameters carry the differences between sub-populations and determine which effects prevail.

Next, to simplify the system, we introduce the total amount of macrophages $M_0=\left[M_N\right]+\left[M\right]$. Adding the above equations, we get $\frac{d{M}_0}{dt}=A_M-{\delta }_M{M}_0.$ If we assume initial conditions for $M_0$ to be at the equilibrium $M_0=A_M/{\delta }_M$, then $M_0$ will remain constant at all times. Then we can express naive macrophages as $\left[M_N\right]=M_0-\left[M\right]$ and write the resulting equation for macrophages as follows:

$$\frac{d\left[M\right]}{dt}=\left({\lambda }_{M{\mu }_2}\right.\left[{\mu }_2\right]+{\lambda }_{M{I}_{\gamma }}\left[I_{\gamma }\right]\left.+{\lambda }_{M{T}_h}\left[T_h\right]\right)\left(M_0-\left[M\right]\right)-{\delta }_M\left[M\right].$$

(12)

2.1.5. Cancer Cells

Cancer cells are epithelial cells with abnormally high growth and abnormally small death rate. Additional loss of apoptosis (cell death) in cancer cells is induced by IL-6 [75,77,83,84]. In addition to innate abnormally high proliferation rate ${\lambda }_C$, proliferation in cancer can be stimulated by expression of STAT3 in cancer cells, where STAT3 is activated by cytokines such as IL-6, IL-17, IL-21 and IL-22 [56,85]. On the other hand, cancer development is suppressed by TGF-$\beta$ [56,86,87,88], IL-12 and IFN-$\gamma$ [56]; the suppressive properties of IL-12 are mediated by IFN-$\gamma$ [89] (so it is not directly included in the equation). Cytotoxic T-cells also directly target cancer cells for destruction [56]. In cancer modeling, proliferation is traditionally taken to be proportional to $\left[C\right]\left(1-\left[C\right]/C_0\right)$, where $C_0$ is the total capacity [90,91]. Effect of dendritic cells on cancer cells is only modeled through intermediary agents, such as T-cells, IL-6 and IL-10. Thus, the resulting equation is

$$\frac{d\left[C\right]}{dt}=\left({\lambda }_C+{\lambda }_{C{\mu }_1}\left[{\mu }_1\right]\right)\left[C\right]\left(1-\frac{\left[C\right]}{C_0}\right)-\left({\delta }_{C{G}_{\beta }}\left[G_{\beta }\right]+{\delta }_{C{I}_{\gamma }}\left[I_{\gamma }\right]+{\delta }_{C{T}_C}\left[T_C\right]+{\delta }_C\right)\left[C\right].$$

(13)

2.1.6. Necrotic Cells

We designate cells which go through the process of necrotic cell death as necrotic cells. Since there is a limited amount of resources in the tumor microenvironment, and cells are under pressure, there are always some necrotic cells produced by the tumor. In addition, when activated cytotoxic T-cells kill colorectal cancer cells by expressing high levels of cytokines like IFN-$\gamma$ and FasL [17], a fraction of the cancer cells may go through the stage of first becoming necrotic cells. Therefore, the rate of “production” of the necrotic cells is given by the fraction of dying cancer cells (${\alpha }_{NC}$) and the resulting dynamics can be written as follows:

$$\frac{d\left[N\right]}{dt}={\alpha }_{NC}\left({\delta }_{C{G}_{\beta }}\left[G_{\beta }\right]+{\delta }_{C{I}_{\gamma }}\left[I_{\gamma }\right]+{\delta }_{C{T}_C}\left[T_C\right]+{\delta }_C\right)\left[C\right]-{\delta }_N\left[N\right].$$

(14)

2.2. Non-Dimensionalization and Sensitivity Analysis

For additional numerical stability and to eliminate scale dependence, we perform non-dimensionalization of the system. For variable X converging to a steady-state $X_{\infty }$, we consider non-dimensional variable $\overline{X}=X/X_{\infty }$. Then, $\overline{X}$ satisfies the equation

$$\frac{d\overline{X}}{dt}=F\left(\overline{X},\theta ,t\right).$$

The (first order) solution sensitivity S with respect to the model parameter $\theta ={\left\{{\theta }_i\right\}}_{i=\overline{1,N}}$ is defined as a vector

$$S_i=\frac{d\overline{X}}{d{\theta }_i},i=\overline{1,N}.$$

In general, the sensitivity vector is time dependent and varies for different solutions and parameter sets [92,93,94]. However, here we consider sensitivity at the steady-state of the equation. The sensitivity of each parameter in the neighborhood of a chosen parameter set $\Omega \left(\theta \right)$ is defined as

$$S_i={\int }_{\Omega }{S}_i\left(\theta \right)d\theta ,$$

where the integration is evaluated numerically with sparse grid points [95,96].

We choose three quantities of interest for the sensitivity analysis: amount of cancer cells C, total amount of cells, and a measure of how fast the system is converging to the steady-state. Consider general steady-state system as follows

$$F(X^{★},\theta )=0,$$

where $X^{★}$ is the equilibrium. We then consider a small perturbation to $X^{★}$ as $\overline{X}\left(t\right)=X^{★}+\epsilon X_1\left(t\right)$. The linearized system becomes

$$\frac{d{X}_1\left(t\right)}{dt}=\nabla F(X^{★},\theta )X_1\left(t\right)+O\left(\epsilon \right),$$

where $\nabla F\left(X,\theta \right)$ is the Jacobian matrix of $F\left(X,\theta \right)$ with respect to X. Thus, we have $X_1\left(t\right)\approx e^{\nabla F(X^{★},\theta )t}$ and the minimal eigenvalue $min\lambda (\nabla F(X^{★},\theta \left)\right)$ determines how fast it reaches the steady-state.

2.3. Cancer Patients’ Data

In recent years, several tumor deconvolution methods have been developed to estimate the relative abundance of various cell types in a tumor from its gene expression profile. A review of these methods [97] and an application of CIBERSORTx on renal cancer [98] show a great performance of CIBERSORTx model. To identify the immune profiles of colonic tumors, we applied CIBERSORTx [99] on RNA-seq gene expression profiles of primary tumors of patients with colon cancer from the Cancer Genome Atlas (TCGA) project of Colon Adenocarcinoma (COAD) downloaded from University of California Santa Cruz (UCSC) Xena web portal. There are a total of 329 patients with RSEM normalized RNA-seq data in ${log}_2$ scale. Before applying CIBERSORTx on this data set, we transformed the gene expression values to the linear space.

2.4. Numerical Methods

In order to solve the time dependent system, we employ the SciPy odeint function [100] using initial conditions based on patients with the smallest tumor area within each cluster. The sensitivity analysis of the system based on the cancer and total cell density at steady-state is obtained analytically by differentiating the steady-state equation with respect to the parameters, namely,

$$\nabla F(X^{★},\theta )\frac{d{X}^{*}}{d\theta }+\frac{\partial F(X^{★},\theta )}{\partial \theta }=0.$$

Then to obtain the sensitivity, $\frac{d{X}^{*}}{d\theta }$, one just needs to numerically invert the matrix $\nabla F$. On the other hand, it is hard to analytically obtain the sensitivity of the eigenvalue, so instead a finite-difference approach is used as follows:

$$\frac{dmin\lambda (\nabla F(X^{★},\theta \left)\right)}{d\theta }\approx \frac{min\lambda \left(\nabla F\left(X^{★},\theta +\frac{1}{2}\Delta \theta \right)\right)-min\lambda \left(\nabla F\left(X^{★},\theta -\frac{1}{2}\Delta \theta \right)\right)}{\Delta \theta },$$

where $\Delta \theta$ is a small discretization parameter.

3. Results

We derived an ODE system describing complex dynamics in the colon cancer microenvironment. Assuming non-negative values of all parameters and non-negative initial conditions, the solution of the system remains non-negative and globally bounded (see Appendix A).

3.1. Patient Data Analysis

We downloaded TCGA clinical data, which includes tumor dimension, stage, gender, vital and tumor status at last follow up, as well as gene expression profiles of primary tumors for patients with colon cancer from the Genomic Data Commons (GDC) portal. We applied CIBERSORTx B-mode on gene expression profiles to estimate the fraction of each immune cell type in each tumor. Elbow method applied on estimated cell fractions (Figure 2A) showed the existence of five distinct immune patterns. We hence performed K-means clustering with $K=5$, in order to group patients based on the immune pattern of their primary tumors. Figure 2B shows average cell fractions for the patients in each of the five clusters. In order to demonstrate the immune variation between these clusters, the average frequencies shown in this plot are of immune cells that have high discrepancy in abundance between clusters. To investigate the effect of these immune patterns on the dynamics of tumors, we model each cluster separately, and based on the steady-state assumptions (see Appendix B). We assume that tumors in each cluster might behave differently not only because of immune cell variations, but also because of variations in their parameter values. For this reason, we estimate parameter values for each cluster separately based on steady-state values derived from patient data, as described further below. The effects of variations in parameter values are investigated through sensitivity and dynamic analyses.

The deconvolution data, described in Section 2.3, only provide the ratios of immune cells in the tumor microenvironment. These data are utilized to obtain the values of variables as detailed in Appendix C. For each patient P, we define their size of tumor, size$\left(P\right)$, to be the product of the longest and the shortest dimensions of the tumor, and we assume total cell density is proportional to the size of the tumor:

$$\mathrm{Total}\_\mathrm{Cell}\_{\mathrm{Density}}_P={\alpha }_{dim}\frac{\mathrm{size}\left(P\right)}{\frac{1}{K}{\sum }_{\mathrm{all}P}\mathrm{size}\left(P\right)}.$$

Then, we take each immune cell value from deconvolution multiplied by $0.4{\alpha }_{dim}\sum \left(\mathrm{Immune}\mathrm{cell}\mathrm{ratios}\right)$ and

$$C=\frac{2}{3}\left(\mathrm{Total}\_\mathrm{Cell}\_\mathrm{Density}-\mathrm{Total}\_\mathrm{Immune}\_\mathrm{Density}\right),N=0.5C.$$

while macrophage capacity $M_0$ is derived from the data, we assume cancer capacity to be $C_0=2*C$ for both mean-based and extreme-based data. We choose the scaling factor ${\alpha }_{dim}=1.125\times {10}^5$ to approximately match the average density of cancer cells across all patients to $4.5\times {10}^4$ cells/cm³ reported in [101]. However, it is important to note that this is no more than scaling and has no effect on the dynamics of the dimensionless system.

We further investigate the clinical features of the clusters to see if there are other differences between clusters rather than ratio of immune cells. Although distributions of gender (Figure 3E,F), tumor dimension (multiplication of the longest and the shortest dimension) (Figure 3D), density of cancer cells, and ratio of cancer to immune cells (Figure 3F) show similar trends in each cluster, we observed some differences in clinical outcomes between the clusters. For example, cluster 5 has the highest percentage of alive patients and tumor-free patients, while cluster 4, which has the lowest M0 macrophages and the highest M2 macrophages, has the highest percentage of patients with tumor at the last time of follow up. In addition, a Chi-squared test shows that cluster 4 has significantly different percentage of tumor status compared to clusters 1 and 5 with p-values 0.0002 and 0.0006, respectively. In addition to these, cluster 3, which has the highest frequency of M0 macrophages, has the highest proportion of deceased patients and stages III and IV tumors. These observations suggest that these clusters represent different immune patterns that might lead to different outcomes. Importantly, no immune cell has a correlation higher than 0.6 with any other immune cells or cancer cells in each cluster. Moreover, Figure 3F demonstrates that as the size of a tumor increases, the ratio of cancer cells over the total number of cells increases.

For each cluster, we consider the mean of variables of patients with tumor size above the average of their cluster as the steady-state values of the variables for the corresponding cluster. The resulting data are given in

Table 2. Steady-state cell densities. We group large tumors in each cluster and calculate their average cell densities in cells/cm³. These data are used for parameter derivation detailed in Appendix B.

Cluster	${\mathbf{T}}_{\mathbf{N}}^{\mathbf{\infty }}$	${\mathbf{T}}_{\mathbf{h}}^{\mathbf{\infty }}$	${\mathbf{T}}_{\mathbf{C}}^{\mathbf{\infty }}$	${\mathbf{T}}_{\mathbf{r}}^{\mathbf{\infty }}$	${\mathbf{D}}_{\mathbf{N}}^{\mathbf{\infty }}$	${\mathbf{D}}^{\mathbf{\infty }}$	${\mathbf{M}}^{\mathbf{\infty }}$	${\mathbf{M}}_0$
1	1.4914 $\times {10}^4$	4.6358 $\times {10}^3$	2.5845 $\times {10}^3$	2.3891 $\times {10}^3$	3.0504 $\times {10}^2$	6.0214 $\times {10}^2$	1.1798 $\times {10}^4$	2.1004 $\times {10}^4$
2	1.1429 $\times {10}^4$	6.0411 $\times {10}^3$	5.3853 $\times {10}^3$	3.3646 $\times {10}^3$	1.0329 $\times {10}^2$	5.1299 $\times {10}^2$	8.6227 $\times {10}^3$	1.6445 $\times {10}^4$
3	9.2381 $\times {10}^3$	1.3864 $\times {10}^3$	1.1139 $\times {10}^3$	2.7910 $\times {10}^3$	1.8878 $\times {10}^{-1}$	1.8635 $\times {10}^2$	6.7972 $\times {10}^3$	3.2146 $\times {10}^4$
4	1.3878 $\times {10}^4$	2.4910 $\times {10}^3$	3.2172 $\times {10}^3$	2.2783 $\times {10}^3$	1.4196 $\times {10}^2$	6.2154 $\times {10}^2$	1.2931 $\times {10}^4$	1.5761 $\times {10}^4$
5	1.0262 $\times {10}^4$	3.7844 $\times {10}^3$	1.6853 $\times {10}^3$	2.6394 $\times {10}^3$	8.0199 $\times 10$	1.9084 $\times {10}^2$	1.1603 $\times {10}^4$	2.8198 $\times {10}^4$
	${\mathbf{C}}^{\infty }$	${\mathbf{N}}^{\infty }$	${\mathit{\mu }}_{\mathbf{1}}^{\infty }$	${\mathit{\mu }}_{\mathbf{2}}^{\infty }$	${\mathbf{H}}^{\infty }$	${\mathbf{I}}_{\mathit{\gamma }}^{\infty }$	${\mathbf{G}}_{\mathit{\beta }}^{\infty }$
1	9.1531 $\times {10}^4$	4.5765 $\times {10}^4$	1.6328 $\times {10}^2$	1.2987 $\times {10}^3$	8.9811 $\times {10}^3$	8.5737	1.9037 $\times {10}^4$
2	9.7064 $\times {10}^4$	4.8532 $\times {10}^4$	1.7552 $\times {10}^2$	1.3249 $\times {10}^3$	8.5279 $\times {10}^3$	10.5677	2.2275 $\times {10}^4$
3	9.0029 $\times {10}^4$	4.5014 $\times {10}^4$	1.9866 $\times {10}^2$	1.2906 $\times {10}^3$	9.5122 $\times {10}^3$	0.8287	2.5145 $\times {10}^4$
4	9.6956 $\times {10}^4$	4.8478 $\times {10}^4$	1.1410 $\times {10}^2$	3.3689 $\times {10}^2$	5.1782 $\times {10}^3$	1.2703	8.1734 $\times {10}^3$
5	8.0584 $\times {10}^4$	4.0292 $\times {10}^4$	1.2058 $\times {10}^2$	7.1551 $\times {10}^2$	7.7848 $\times {10}^3$	5.7892	2.6260 $\times {10}^4$

. These data are only used as in Appendix B to generate the parameter sets given in

Table A2. Scaling-dependent parameters. Non-dimensional values of scaling-dependent parameters (in day⁻¹, because the time was not scaled) for each cluster derived from the steady-state assumptions and patient data.

Parameter	Cluster 1	Cluster 2	Cluster 3	Cluster 4	Cluster 5
${\overline{\lambda }}_{T_{h}D}$	2.1399	2.7619	2.3163	1.7613	2.0179
${\overline{\lambda }}_{T_{h}M}$	3.9171 $\times {10}^{-2}$	4.3371 $\times {10}^{-2}$	7.8933 $\times {10}^{-2}$	3.4235 $\times {10}^{-2}$	1.1462 $\times {10}^{-1}$
${\overline{\lambda }}_{T_h{\mu }_1}$	1.1871 $\times {10}^{-2}$	1.9332 $\times {10}^{-2}$	5.0516 $\times {10}^{-2}$	6.6146 $\times {10}^{-3}$	2.6083 $\times {10}^{-2}$
${\overline{\lambda }}_{T_C{T}_h}$	3.2503	4.4295	3.6087	2.3377	3.5400
${\overline{\lambda }}_{T_{C}D}$	6.0053 $\times {10}^{-1}$	5.3505 $\times {10}^{-1}$	6.8998 $\times {10}^{-1}$	8.2970 $\times {10}^{-1}$	2.5392 $\times {10}^{-1}$
${\overline{\lambda }}_{T_r{T}_h}$	8.2898 $\times {10}^{-1}$	8.9260 $\times {10}^{-1}$	7.1048 $\times {10}^{-1}$	6.6821 $\times {10}^{-1}$	6.4634 $\times {10}^{-1}$
${\overline{\lambda }}_{T_r{\mu }_2}$	5.1120 $\times {10}^{-2}$	4.3090 $\times {10}^{-2}$	1.4558 $\times {10}^{-1}$	1.9891 $\times {10}^{-2}$	2.6898 $\times {10}^{-2}$
${\overline{\lambda }}_{T_r{G}_{\beta }}$	6.4120 $\times {10}^{-2}$	6.1990 $\times {10}^{-2}$	2.4271 $\times {10}^{-1}$	4.1296 $\times {10}^{-2}$	8.4474 $\times {10}^{-2}$
${\overline{\lambda }}_{DH}$	2.6134 $\times {10}^{-1}$	2.5410 $\times {10}^{-1}$	2.6749 $\times {10}^{-1}$	2.0531 $\times {10}^{-1}$	2.5253 $\times {10}^{-1}$
${\overline{\lambda }}_{DC}$	1.0129 $\times {10}^{-1}$	1.0998 $\times {10}^{-1}$	9.6275 $\times {10}^{-2}$	1.4619 $\times {10}^{-1}$	9.9406 $\times {10}^{-2}$
${\overline{\lambda }}_{M{\mu }_2}$	6.0389 $\times {10}^{-4}$	4.0918 $\times {10}^{-4}$	4.0357 $\times {10}^{-4}$	1.0688 $\times {10}^{-3}$	2.2488 $\times {10}^{-4}$
${\overline{\lambda }}_{M{I}_{\gamma }}$	5.4536 $\times {10}^{-4}$	4.4647 $\times {10}^{-4}$	3.5449 $\times {10}^{-5}$	5.5129 $\times {10}^{-4}$	2.4891 $\times {10}^{-4}$
${\overline{\lambda }}_{M{T}_h}$	2.4482 $\times {10}^{-2}$	2.1190 $\times {10}^{-2}$	4.9239 $\times {10}^{-3}$	8.9758 $\times {10}^{-2}$	1.3509 $\times {10}^{-2}$
${\overline{\lambda }}_{C{\mu }_1}$	3.6453 $\times {10}^{-3}$	5.6528 $\times {10}^{-3}$	1.3143 $\times {10}^{-3}$	2.5321 $\times {10}^{-3}$	2.5536 $\times {10}^{-3}$
${\overline{\alpha }}_{NC}$	1.5	1.5	1.5	1.5	1.5
${\overline{\lambda }}_{{\mu }_1{T}_h}$	9.5154 $\times {10}^{-2}$	1.5848 $\times {10}^{-1}$	5.1746 $\times {10}^{-2}$	4.8884 $\times {10}^{-2}$	8.0520 $\times {10}^{-2}$
${\overline{\lambda }}_{{\mu }_{1}M}$	9.6867 $\times {10}^{-1}$	9.0479 $\times {10}^{-1}$	1.0148	1.0150	9.8745 $\times {10}^{-1}$
${\overline{\lambda }}_{{\mu }_{1}D}$	6.1798 $\times {10}^{-3}$	6.7287 $\times {10}^{-3}$	3.4777 $\times {10}^{-3}$	6.0984 $\times {10}^{-3}$	2.0302 $\times {10}^{-3}$
${\overline{\lambda }}_{{\mu }_{2}M}$	3.6856	3.1869	3.2127	3.7737	3.7140
${\overline{\lambda }}_{{\mu }_{2}D}$	1.8810 $\times {10}^{-1}$	1.8960 $\times {10}^{-1}$	8.8081 $\times {10}^{-2}$	1.8139 $\times {10}^{-1}$	6.1088 $\times {10}^{-2}$
${\overline{\lambda }}_{{\mu }_2{T}_r}$	7.4634 $\times {10}^{-1}$	1.2435	1.3192	6.6490 $\times {10}^{-1}$	8.4487 $\times {10}^{-1}$
${\overline{\lambda }}_{HN}$	5.6645 $\times 10$	5.6824 $\times 10$	5.7494 $\times 10$	5.6720 $\times 10$	5.6504 $\times 10$
${\overline{\lambda }}_{HM}$	1.4603	1.0096	8.6817 $\times {10}^{-1}$	1.5130	1.6271
${\overline{\lambda }}_{H{T}_h}$	2.8689 $\times {10}^{-1}$	3.5366 $\times {10}^{-1}$	8.8539 $\times {10}^{-2}$	1.4573 $\times {10}^{-1}$	2.6536 $\times {10}^{-1}$
${\overline{\lambda }}_{H{T}_C}$	1.5995 $\times {10}^{-1}$	3.1527 $\times {10}^{-1}$	7.1133 $\times {10}^{-2}$	1.8821 $\times {10}^{-1}$	1.1817 $\times {10}^{-1}$
${\overline{\lambda }}_{H{T}_r}$	1.4786 $\times {10}^{-1}$	1.9697 $\times {10}^{-1}$	1.7824 $\times {10}^{-1}$	1.3328 $\times {10}^{-1}$	1.8507 $\times {10}^{-1}$
${\overline{\lambda }}_{I_{\gamma }{T}_h}$	8.8980	6.8580	6.4053	4.6181	9.8011
${\overline{\lambda }}_{I_{\gamma }{T}_C}$	1.9843 $\times 10$	2.4454 $\times 10$	2.0584 $\times 10$	2.3857 $\times 10$	1.7459 $\times 10$
${\overline{\lambda }}_{I_{\gamma }M}$	4.5290	1.9577	6.2806	4.7945	6.0098
${\overline{\lambda }}_{G_{\beta }M}$	4.1497 $\times {10}^2$	3.5894 $\times {10}^2$	3.5375 $\times {10}^2$	4.2425 $\times {10}^2$	4.0652 $\times {10}^2$
${\overline{\lambda }}_{G_{\beta }{T}_r}$	8.4033 $\times 10$	1.4006 $\times {10}^2$	1.4525 $\times {10}^2$	7.4749 $\times 10$	9.2476 $\times 10$
${\overline{\delta }}_{T_h{\mu }_2}$	4.2911 $\times {10}^{-1}$	4.3775 $\times {10}^{-1}$	4.2642 $\times {10}^{-1}$	1.1131 $\times {10}^{-1}$	2.3641 $\times {10}^{-1}$
${\overline{\delta }}_{T_h{T}_r}$	1.5309	2.1559	1.7884	1.4599	1.6912
${\overline{\delta }}_{T_C{\mu }_2}$	7.5419 $\times {10}^{-1}$	7.6938 $\times {10}^{-1}$	7.4947 $\times {10}^{-1}$	1.9563 $\times {10}^{-1}$	4.1550 $\times {10}^{-1}$
${\overline{\delta }}_{T_C{T}_r}$	2.6906	3.7891	3.1432	2.5658	2.9724
${\overline{\delta }}_{T_r{\mu }_1}$	7.1322 $\times {10}^{-1}$	7.6668 $\times {10}^{-1}$	8.6777 $\times {10}^{-1}$	4.9839 $\times {10}^{-1}$	5.2671 $\times {10}^{-1}$
${\overline{\delta }}_{DH}$	3.3577 $\times {10}^{-2}$	3.1883 $\times {10}^{-2}$	3.5563 $\times {10}^{-2}$	1.9360 $\times {10}^{-2}$	2.9105 $\times {10}^{-2}$
${\overline{\delta }}_{DC}$	5.2053 $\times {10}^{-2}$	5.5200 $\times {10}^{-2}$	5.1199 $\times {10}^{-2}$	5.5139 $\times {10}^{-2}$	4.5828 $\times {10}^{-2}$
${\overline{\delta }}_{C{G}_{\beta }}$	9.1600 $\times {10}^{-4}$	7.1011 $\times {10}^{-4}$	1.8590 $\times {10}^{-3}$	3.9506 $\times {10}^{-4}$	1.3131 $\times {10}^{-3}$
${\overline{\delta }}_{C{I}_{\gamma }}$	6.5951 $\times {10}^{-4}$	5.3859 $\times {10}^{-4}$	9.7944 $\times {10}^{-5}$	9.8156 $\times {10}^{-5}$	4.6278 $\times {10}^{-4}$
${\overline{\delta }}_{C{T}_C}$	2.1270 $\times {10}^{-3}$	2.9365 $\times {10}^{-3}$	1.4085 $\times {10}^{-3}$	2.6597 $\times {10}^{-3}$	1.4414 $\times {10}^{-3}$
${\overline{A}}_{T_N}$	1.5006	4.1268	1.2183	1.1785	1.6150
${\overline{A}}_{Dn}$	1.0264	2.1172	3.5941 $\times {10}^2$	1.8353	1.1435
${\overline{M}}_0$	1.7803	1.9072	4.7293	1.2189	2.4303
${\overline{C}}_0$	2.0	2.0	2.0	2.0	2.0
${\alpha }_{T_N{T}_h}$	3.1084 $\times {10}^{-1}$	5.2856 $\times {10}^{-1}$	1.5008 $\times {10}^{-1}$	1.7950 $\times {10}^{-1}$	3.6879 $\times {10}^{-1}$
${\alpha }_{T_N{T}_C}$	1.7330 $\times {10}^{-1}$	4.7118 $\times {10}^{-1}$	1.2057 $\times {10}^{-1}$	2.3182 $\times {10}^{-1}$	1.6423 $\times {10}^{-1}$
${\alpha }_{T_N{T}_r}$	1.6020 $\times {10}^{-1}$	2.9438 $\times {10}^{-1}$	3.0212 $\times {10}^{-1}$	1.6417 $\times {10}^{-1}$	2.5720 $\times {10}^{-1}$
${\alpha }_{D_{N}D}$	1.9739	4.9667	9.8717 $\times {10}^2$	4.3783	2.3796

and

Table A3. Scaling-independent parameters. Values of scaling-independent parameters (in day${}^{-1}$) for each cluster derived from the steady-state assumptions and patient data.

Parameter	Cluster 1	Cluster 2	Cluster 3	Cluster 4	Cluster 5
${\lambda }_C$	5.3323 $\times {10}^{-3}$	3.7596 $\times {10}^{-3}$	7.8327 $\times {10}^{-3}$	5.3535 $\times {10}^{-3}$	5.5150 $\times {10}^{-3}$
${\delta }_C$	7.8626 $\times {10}^{-4}$	5.2094 $\times {10}^{-4}$	1.2081 $\times {10}^{-3}$	7.8983 $\times {10}^{-4}$	8.1708 $\times {10}^{-4}$
${\delta }_N$	6.7332 $\times {10}^{-3}$	7.0593 $\times {10}^{-3}$	6.8602 $\times {10}^{-3}$	5.9142 $\times {10}^{-3}$	6.0514 $\times {10}^{-3}$

. In other words, we assume the largest tumors in each cluster are near their steady-state. Since we have estimated the values of all variables for all tumors from their gene expression profiles, we have the values of all variables at the steady-state. To estimate the parameters, we first set the left side of the ODEs to zero, in order to create a system of equations that does not depend on time. In order to obtain unique solutions, we need to have the same number of equations as the number of parameters. Therefore, we estimate a small number of parameters using biological studies. We also assume some relations among parameters as described in Appendix B. Hence, by these assumptions, we obtain a unique set of parameter values for each cluster, and therefore all parameters are identifiable.

In order to investigate the dynamics of tumors and solve the system of ODEs, we need the initial values of variables. We assume the smallest tumors in each cluster represent the initial conditions. Therefore, we use the values of the variables of the smallest tumors (estimated using their gene expression profiles) as the initial conditions. The relative values are given in

Table 3. Dimensionless initial conditions. Values of initial conditions for the dimensionless system derived from the patients with the smallest tumor size.

Cluster	${\mathbf{T}}_{\mathbf{N}}/{\mathbf{T}}_{\mathbf{N}}^{\mathbf{\infty }}$	${\mathbf{T}}_{\mathbf{h}}/{\mathbf{T}}_{\mathbf{h}}^{\mathbf{\infty }}$	${\mathbf{T}}_{\mathbf{C}}/{\mathbf{T}}_{\mathbf{C}}^{\mathbf{\infty }}$	${\mathbf{T}}_{\mathbf{r}}/{\mathbf{T}}_{\mathbf{r}}^{\mathbf{\infty }}$	${\mathbf{D}}_{\mathbf{N}}/{\mathbf{D}}_{\mathbf{N}}^{\mathbf{\infty }}$	$\mathbf{D}/{\mathbf{D}}^{\mathbf{\infty }}$	$\mathbf{M}/{\mathbf{M}}^{\mathbf{\infty }}$
1	0.9311	1.2492	2.4626	0.6872	1.6328	0.0003	0.6737
2	1.2302	1.3155	1.5210	0.5107	2.0461	2.7822	1.2920
3	1.1997	0.8555	1.6948 $\times {10}^{-4}$	0.6572	1.0000	1.0130 $\times {10}^{-3}$	1.4150
4	1.4471	0.1571	0.5823	0.8910	5.6827	4.2945	0.9259
5	0.6794	2.6119	1.6294	1.8819	2.3538 $\times {10}^{-3}$	0.4542	0.7749
	$\mathbf{C}/{\mathbf{C}}^{\infty }$	$\mathbf{N}/{\mathbf{N}}^{\infty }$	${\mathit{\mu }}_{\mathbf{1}}/{\mathit{\mu }}_{\mathbf{1}}^{\infty }$	${\mathit{\mu }}_{\mathbf{2}}/{\mathit{\mu }}_{\mathbf{2}}^{\infty }$	$\mathbf{H}/{\mathbf{H}}^{\infty }$	${\mathbf{I}}_{\mathit{\gamma }}/{\mathbf{I}}_{\mathit{\gamma }}^{\infty }$	${\mathbf{G}}_{\mathit{\beta }}/{\mathbf{G}}_{\mathit{\beta }}^{\infty }$
1	3.1466 $\times {10}^{-4}$	0.0	0.4971	0.5124	1.4712	3.8892	0.2549
2	2.9672 $\times {10}^{-4}$	0.0	0.7578	0.1790	0.6036	0.9385	0.5566
3	3.1991 $\times {10}^{-4}$	0.0	0.1335	0.8419	1.2566	0.0	0.6851
4	2.9706 $\times {10}^{-4}$	0.0	0.4137	5.7720	1.4630	0.0	2.5629
5	3.5741 $\times {10}^{-4}$	0.0	0.4587	2.2979	1.1835	0.4084	0.3457

. The dynamics with initial conditions based on other patients are presented in Appendix C.

3.2. Sensitivity Analysis

We perform sensitivity analysis of the non-dimensionalized system with parameters derived from patient data through steady-state assumptions. Table 2 contains the steady-state values used for each cluster, and Appendix B shows the parameter derivation and non-dimensionalization in detail. It is worth pointing out that there might be a variation in the calculated parameters due to differences between patients and possible alterations in assumptions of Appendix B. In order to account for those possible variations, sensitivity analysis is performed on dimension-less system of equations (see Section 2.2). We use cancer cells, total cell density and minimal eigenvalue of the Jacobian of the ODE system as the variables of interest in the sensitivity analysis. Minimal eigenvalue of the Jacobian serves as a measure of how fast the system converges to the steady-state. Figure 4A shows the four most sensitive parameters for each cluster. Additionally, to evaluate the effect of immune microenvironment on cancer, we look at the sensitivity of cancer cells and total cell density excluding the parameters appearing in the equations for cancer and necrotic cells. The resulting data denoted as “Immune sensitivity” are given in Figure 4B.

Across all clusters, the most sensitive parameters are cancer proliferation and death rates directly present in the cancer Equation (13). From the third column in Figure 4A, we conclude that for all clusters, increased cancer proliferation coefficients correspond to faster convergence to the steady-state, while increased cancer death rates lead to a slower convergence. When considering immune sensitivity presented on Figure 4B, in clusters 1, 2, 3 and 5, the most sensitive immune parameters are those corresponding to the activation and decay rates of macrophages, with only sensitivity levels being different between clusters and variables. In clusters 1 and 2, which include tumors with a smaller density of naive macrophages than activated macrophages, an increase in decay rate of macrophages causes a decrease in the density of cancer cells and total cell density. On the other hand, an increase in any of the activation rates for macrophages causes an increase in both quantities of interest. However, for clusters 3 and 5, which include tumors with a higher density of naive macrophages than activated macrophages, the effects are reversed. Interestingly, for cluster 3, the increase in macrophage activation rate results in both lower cancer cell density and total cell density, with latter sensitivity being noticeably smaller by absolute value. On the other hand, for cluster 5 the increase in macrophage activation rate results in lower cancer cell density, but higher total cell density. This can be explained by a significant increase in immune cell density, which for cluster 5 is even higher than the corresponding decrease in cancer cell density. All these results demonstrate that at the steady-state, tumor-associated macrophages could have different effects on different clusters of patients depending on their immune profile.

The outlying cluster 4, which consists of tumors with a significantly small density of naive macrophages compared to the other clusters, is less sensitive to the activation rates of macrophages. The most sensitive immune parameters for cancer cell density are those related to the activation and degradation of regulatory T-cells. The results indicate that increased regulatory response activation rate corresponds to an increase in the cancer cell density, while an increase in T-reg cell degradation rate results in a decrease in cancer cell density, demonstrating that for this cluster of patients regulatory T-cells have mostly negative effects. Importantly, the most sensitive parameter for the total cell density is still the decay rate of macrophages, and macrophages still have a negative effect, i.e., the faster decay of macrophages leads to the smaller tumors in the steady-state.

3.3. Dynamic of Tumor Microenvironment

We investigate the dynamics of each variable, with parameters derived for each cluster based on steady-state assumptions (see Table 2 for steady-state values and

Table A1. Prescribed parameters and their references. Innate degradation and death rates (in day⁻¹) derived or adopted from given references.

Parameter	${\mathit{\delta }}_{{\mathit{T}}_{\mathit{N}}}$	${\mathit{\delta }}_{{\mathit{T}}_{\mathit{h}}}$	${\mathit{\delta }}_{{\mathit{T}}_{\mathit{C}}}$	${\mathit{\delta }}_{{\mathit{T}}_{\mathit{r}}}$	${\mathit{\delta }}_{\mathit{D}}$	${\mathit{\delta }}_{\mathit{M}}$	${\mathit{\delta }}_{{\mathit{\mu }}_1}$	${\mathit{\delta }}_{{\mathit{\mu }}_2}$	${\mathit{\delta }}_{\mathit{H}}$	${\mathit{\delta }}_{{\mathit{I}}_{\mathit{\gamma }}}$	${\mathit{\delta }}_{{\mathit{G}}_{\mathit{\beta }}}$
Value	9.4951 $\times {10}^{-4}$	0.231	0.406	0.231	0.277	0.02	1.07	4.62	58.7	33.27	499.0
Reference	[154]	[150]	[150]	[150]	[149]	[151,158]	[153]	[155]	[152]	[157]	[156]

, Table A2, Table A3 for parameter values) and initial conditions of patients with the smallest tumor (see Table 3). Figure 5 and Figure 6 show the dynamics of cell densities and cytokines expressions respectively. We investigate the effect of variations in parameters’ values on the dynamics of the tumor by varying the most sensitive parameters by 10% in both the positive and negative directions. These variations are shown as shaded regions on each of the graphs.

For most clusters, cancer cells grow as helper T-cells, cytotoxic cells (cytotoxic T-cells and NK cells), dendritic cells and macrophages increase in density over time, while naive T-cells, regulatory T-cells and naive dendritic cells decrease in density. The increase in cytotoxic cells along with tumor progression is somewhat contradicting to the finding in [102,103] that colon primary tumor growth is associated with decreased cytotoxic T-cells density. However, there is no correlation between tumor size and cytotoxic cells in the TCGA data of colonic primary tumors. Moreover, it is important to note that in our model cancer cells’ growth is multiple times faster than the rate of change of any immune cells (Figure 5). Thus, even though cytotoxic cells density grows over time, the tumor is growing at a much faster rate. Since tumor cells activate dendritic cells which then activate cytotoxic cells, it is reasonable to see some growth of cytotoxic cells when tumor cells density increases rapidly.

Clusters 2 and 4 have the highest cancer cell density at steady-state and also the highest growth rate of cancer cells. Cluster 2’s cancer cells start out with the lowest growth rate, but at around 1800 days grow significantly faster and end up growing the fastest among all clusters. Cluster 2 has the highest density of helper T-cells and cytotoxic cells, both in the early stages of cancer development and at steady-state, as well as the highest growth rate of these cells. However, cluster 2 has rather low density and low growth rate of macrophages.

Cluster 4, having the largest density of activated macrophages and a significantly small density of naive macrophages (Figure 2B), first demonstrates average cancer growth rate, but then increases and has one of the two highest cancer cell densities at the steady-state (Figure 5). Similar to cluster 2, cluster 4 has high density of cytotoxic cells (CD8 T-cells and NK cells) initially and at steady-state. Both clusters 2 and 4 have a low growth rate of macrophages and a high density of dendritic cells, compared to other clusters. Immune cell dynamics of clusters 2 and 4 demonstrate that a high density of cytotoxic cells and dendritic cells, along with a low growth rate of macrophages, correlates with a high growth rate of cancer cells.

However, unlike cluster 2, cluster 4 has a low growth rate of cytotoxic cells and helper T-cells, and a low density of helper T-cells overall. Though both clusters 2 and 4 have a low growth rate of macrophages, cluster 4 has the highest density of macrophages among all clusters, while cluster 2 has the second lowest macrophages density. Regulatory T-cells also behave very differently between cluster 2 and cluster 4. Cluster 2 has high density and high decline rate of regulatory T-cells over time, but cluster 4 has both low density and low decline rate of this cell. These observations suggest that cell densities alone cannot predict cancer progression and there are no specific biomarkers that are sufficient to model tumor growth. Instead, a time series immune interaction network with tumor cells can be useful in modeling cancer development.

Cluster 5, with the density of activated macrophages being slightly less than naive macrophages (Figure 2B), has the lowest cancer cell density at steady-state and the lowest cancer cell growth of all clusters (Figure 5). This cluster has the lowest growth rate and density at initial condition and steady-state of naive dendritic cells, activated dendritic cells and cytotoxic cells, except for cytotoxic cells density at steady-state (second lowest). It also has the highest growth rate of macrophages among the five clusters. This observation might imply that slow tumor growth is associated with low density and growth rate of naive and activated dendritic cells, cytotoxic cells and high growth rate of macrophages.

Cluster 1, which is characterized by the second largest population of macrophages and helper T-cells, demonstrates that dendritic cells alone cannot be chosen as a marker of cancer progression, as it has the second highest dendritic cell population, but only the third highest cancer cell population at the steady-state, being surpassed by cluster 2.

Cluster 3, being a clear outlier in the immune dynamics, has near zero density of naive dendritic cells. This alone prevents it from creating significant variations in the immune response during the cancer progression. It is interesting to note, that while almost unchecked by immune responses, this cluster initially demonstrates noticeably highest cancer growth rate, but results in the second lowest cancer density at the steady-state.

Tumor cytokines’ dynamics (Figure 6) indicate that as the tumor grows, HMGB1, IFN-$\gamma$ and ${\mu }_1$ (IL-6, IL-17, IL-21, IL-22) increase in density, but TGF-$\beta$ and ${\mu }_2$ (IL-10, CCL20) stay relatively constant. Clusters 2 and 4, which have the highest cancer cell growth rate among all clusters, show different cytokines’ behaviors throughout time. At steady-state, cluster 4 has significantly lower densities of all cytokines in our model than cluster 2, despite the fact that they have the same cancer cell density then. Cluster 4 also has a much lower growth rate of ${\mu }_1$, HMGB1 and IFN-$\gamma$ compared to cluster 2. Clusters 1 and 5 have more similar growth rates of these cytokines as cluster 2, even though they have rather different tumor growth rates from cluster 2. Thus, the density or growth of any specific cytokine is not an adequate predictor of tumor progression, and we need the full interaction network to effectively model the cancer cell growth.

Additionally within each cluster, we look at the dynamics of cancer and total cell density with different initial conditions, each derived from a different patient in that cluster. See Appendix C, and specifically Figure A2, Figure A3, Figure A4, Figure A5, Figure A6, for more details on different initial conditions and resulting dynamics. This result indicates that even within the same cluster, different initial immune profiles may cause dramatic differences in the cancer progression rate. Additionally, while the dynamics of cancer cell density remains monotone across all patients, we observe oscillatory behavior in the total cell density. This can be explained by a temporary surge of immune cell density at the early stages of cancer, which also appears to correlate with slower cancer progression rate. The only cluster which does not exhibit this oscillatory behavior is cluster 3. As mentioned before, due to lack of naive dendritic cells, cluster 3 does not show significant immune cell density variations, which are the source of oscillatory behavior for other clusters.

4. Discussion

Although there are many mathematical models for cancer [34,35,36,37,38,39,40,41,42,43,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124], one of the outstanding challenges in mathematical modeling of cancers is the existence of many unknown parameters and the limited number of data sets. For this reason, the approach of many of these mathematical models is to assume some values for parameters, use estimated parameters from other diseases, or vary the parameters and initial conditions within biologically-feasible values in order to investigate their effects on the results. Here, we choose some of the parameter values based on biological studies and the rest by utilizing patients’ gene expression data. New advances in tumor deconvolution techniques help us to utilize cancer patients’ data in order to develop a data driven mathematical model of tumor growth. Using tumor deconvolution methods, we estimate the relative abundance of various cell types from gene expression profiles of tumors. The machine learning algorithm of K-means clustering indicates the existence of five distinct groups of colon cancers based on their immune patterns. The comparison of tumor behaviors in these groups suggests that the dynamics of tumors strongly depends on their immune structure.

While it would be ideal to use time course gene expression data of colon cancer patients in our framework, the availability of these time series data sets is limited. In order to combat this limitation, clustering was used to group patients with similar immune patterns and treat each group as time course data based on the size of tumor, which means the data points with small tumor density are considered data from early stages (initial conditions) and the data points with large tumor density are considered data from late stages (steady-state values). We assume a large tumor in a cluster that is in the steady-state is an evolution of a small tumor in the same cluster, because when a small tumor in a cluster (e.g., 1) evolves to a large tumor in another cluster (e.g., 3), its dynamics will quickly converge to the dynamics of a small tumor in the cluster of the steady-state tumor (e.g., 3) (Figure A1). We then follow a common approach of mathematical biology models that use assumptions on the steady-state values of the system to estimate parameters of the model [125,126]. Note that we use patient data to estimate the steady-state values of each cluster, and we then estimate parameters based on the values at the steady-state. Due to the non-dimensionalization process, the relative dynamics of the system are independent of the data scaling and only depend on relative values for patients.

The parameter values have been estimated based on the assumption that the largest tumors are near the steady-state, because the large tumors do not have much space to grow. Figure 3F shows that the ratio of cancer cells over total immune cells increases when the size of tumors increases. This may suggest that small tumors have a broader potential for decision making and more options to evolve than large tumors. We also investigate the dynamics of all tumors with a size “below average” as initial condition. The resulting dynamics (Figure A2, Figure A3, Figure A4, Figure A5, Figure A6) show that when the size of initial tumors increases, the time to reach the steady-state decreases. This observation is also in support of the steady-state assumption.

To evaluate the effect of each parameter value on the results, we perform a comprehensive sensitivity analysis, which covers a range of parameter values. For all parameters that have not been determined using biological literature, we estimate 5 different values for each parameter in each cluster. We find the most sensitive parameters by applying a gradient based sensitivity analysis method on the dimension-less system of equations. Note, although we used completely different values for some of the parameters for each cluster, our sensitivity analysis shows that a similar set of parameters are the most sensitive parameters for all clusters. This demonstrates that while many parameters are unknown, the evolution of tumors is not sensitive to many of these parameters (Figure 4). For those sensitive parameters, we show how their variations would affect the results of the model (Figure 5 and Figure 6).

The mathematical model shows a unique pattern of tumor growth in each cluster based on their immune infiltration. For example, the model indicates that a high density of cytotoxic T-cells and dendritic cells and a low growth rate of macrophages are associated with a high growth rate of cancer cells (observed in clusters 2 and 4, Figure 4 and Figure 5), while a low density and growth rate of naive and active dendritic cells and cytotoxic T-cells and a high growth rate of macrophages correlate with slow tumor growth (observed in cluster 5, Figure 4 and Figure 5). Clinical information provided for patients in the clusters also supports the results of the dynamical model. Cluster 4, which has the highest percentage of patients with tumors at the time of last follow up (Figure 3C) has poor outcome compared to other clusters, while cluster 5, which has the highest percentages of tumor-free and alive patients (Figure 3A,C) has better outcome than the other clusters.

Importantly, our results imply that macrophages’ activation rates have different effects in different clusters. A high activation rate of macrophages leads to a high density of cancer cells in clusters 1 and 2 (Figure 4), in which there are more activated than naive macrophages. This result is in agreement with an observation of [127], which indicates a high CD206/CD68 ratio is significantly associated with poor outcome. Note, CD206 is the marker for M2 macrophages and CD68 is expressed at high levels on M0 cells [128]. However, the activation rates of macrophages are negatively correlated with tumor growth in clusters 3 and 5 (Figure 4), in which they are more naive than activated macrophages. This is consistent with the observation that a high level of macrophages is associated with a favorable outcome of colon cancer patients in [102,129]. The study in [102] also shows that a high level of regulatory T-cells is related to poor prognosis of patients, which supports our results that regulatory T-cells decrease in density as cancer cells increase in density. Another similar finding between [102] and our study is that the density of dendritic cells increases along with tumor progression (Figure 5).

There is a significant body of research analyzing statistical and mathematical relations of particular components of tumor microenvironment and the disease progression and outcome for subsequent establishment of prognostic biomarkers [130,131,132,133,134,135,136,137,138,139]. Our result demonstrates that the dynamics of cancer development cannot be captured by one specific biomarker, but can rather be characterized by complex time-dependent interactions between many components of the immune system and tumor tissue. It is important to further develop and analyze these tumor–immune cell interactions and how they affect different possibilities of treatment.

It is important to point out that, similar to many other mathematical models of biological processes, this model has some limitations that arise from the lack of time course data sets. Specifically, we make the assumption that the largest tumors in each cluster are the evolution of the smallest ones. Importantly, as we mentioned above, some of the predictions of the model match with biological observations. However, these facts do not provide a complete validation for the model, and the model should be validated on a separate time course data. We hope that this work encourages scientists to collect time course data. Although this work has some limitations, it provides important insights and an opportunity for scientists to improve the mathematical model of colonic tumors and/or validate the model if they have more data or insights.

One way forward is the design of patient-specific models [140,141,142,143]. These models can utilize the tumor immune microenvironment deconvolution and clustering methods for available patient data as detailed in this paper. New prognosis can be built based on established dynamics from patients with similar immune characteristics. To better match the dynamics of the model to real patient data, various parameter fitting algorithms can be utilized [144,145,146,147]. Another possible improvement is a transition to a partial differential equations model [148] to analyze spatial properties of tumor development as well as temporal.