ScholarWorks@UMass Amherst ScholarWorks@UMass Amherst

: 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 proﬁles 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.


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].
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.4 The Role of Inflammation in Colitis-associated Cancer  ]can be written as IL-10 is produced by macrophages [2,18], dendritic cells [24,28] and Treg cells [10,32,58,64].CCL20 is produced by macrophages [14].Thus, the equation for [µ 2 ] is IFN-γ is secreted by sub-population of macrophages [2,13,47,52,70], helper T-cells [8,46] and cytotoxic

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.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 → C, we use the mass action law dC dt = λAB, where λ is production rate of C [44,45].Throughout the paper, we use the symbol λ for production, activation or proliferation rates, and the symbol δ for decay, natural death or premature death (necrosis) rates.Immunosuppresive agents includes effects of IL-10 and CCL20 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 µ 1 .We also combine IL-10 and CCL20 and denote their sum by the variable µ 2 .The cytokines treated as model variables are HMGB1, IFN-γ, TGF-β, IL-6 and IL-10.We then model the dynamics of cytokines in the following way.

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-α 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 δ T h .The resulting equation is 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 [IL-2, 4, 5, 13] 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: Regulatory T-cells can be activated by IL-2 [56,72,73], CCL20 [62] and TGF-β [56,67].IL-6 suppresses T-reg differentiation and shifts it towards T-helper type [74].The resulting dynamics can be described as follow: Combining all activation and introducing independent naive T-cell production A T N , we get the following equation for naive T-cells:

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 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:

Macrophages
There are two main sub-types of macrophages: M1 and M2.M1 phenotype can be activated by IFN-γ, 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 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: 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 Adding the above equations, we get dM 0 dt = A M − δ M M 0 .If we assume initial conditions for M 0 to be at the equilibrium M 0 = A M /δ M , then M 0 will remain constant at all times.Then we can express naive macrophages as [M N ] = M 0 − [M] and write the resulting equation for macrophages as follows: 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 λ 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-β [56,[86][87][88], IL-12 and IFN-γ [56]; the suppressive properties of IL-12 are mediated by IFN-γ [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 [C] (1 − [C] /C 0 ), 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 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-γ 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 (α NC ) and the resulting dynamics can be written as follows:

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 ∞ , we consider non-dimensional variable X = X/X ∞ .Then, X satisfies the equation The (first order) solution sensitivity S with respect to the model parameter θ = {θ i } i=1, N is defined as a vector 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 Ω(θ) is defined as 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 , θ) = 0, where X is the equilibrium.We then consider a small perturbation to X as X(t) = X + εX 1 (t).
The linearized system becomes where ∇F (X, θ) is the Jacobian matrix of F (X, θ) with respect to X. Thus, we have X 1 (t) ≈ e ∇F(X , θ)t and the minimal eigenvalue min λ (∇F (X , θ)) determines how fast it reaches the steady-state.

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.

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, Then to obtain the sensitivity, dX * dθ , one just needs to numerically invert the matrix ∇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: where ∆θ is a small discretization parameter.

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).

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(P), 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: .
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 α dim = 1.125 × 10 5 to approximately match the average density of cancer cells across all patients to 4.5 × 10 4 cells/cm 3 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.These data are only used as in Appendix B to generate the parameter sets given in Tables A2 and A3.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.Table 2. Steady-state cell densities.We group large tumors in each cluster and calculate their average cell densities in cells/cm 3 .These data are used for parameter derivation detailed in Appendix B.  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.The dynamics with initial conditions based on other patients are presented in Appendix C. Table 3. Dimensionless initial conditions.Values of initial conditions for the dimensionless system derived from the patients with the smallest tumor size.

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.

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 Tables A1-A3 for parameter values) and initial conditions of patients with the smallest tumor (see Table 3).Figures 5 and 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-γ and µ 1 (IL-6, IL-17, IL-21, IL-22) increase in density, but TGF-β and µ 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 µ 1 , HMGB1 and IFN-γ 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 Figures A2-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.

Discussion
Although there are many mathematical models for cancer [34][35][36][37][38][39][40][41][42][43], 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 (Figures A2-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 (Figures 5 and 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, Figures 4 and 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, Figures 4 and 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.

Appendix A. ODE System and Analysis
Combining Equations ( 1)-( 14), we obtain the following system The system has 14 variables and 59 different parameters.The λ parameters correspond to proliferation, activation and production rates, δ parameters correspond to degradation and cell death rates and four parameters: A T N and A D N respectively are the production rates of naive T-cells and dendritic cells, M 0 and C 0 are the total density of macrophages (naive and activated together) and cancer cells maximum capacity, respectively.

Appendix A.1. Positivity
To prove that the system with positive coefficients and positive initial conditions has positive solution, let us consider the set of integrating factors, one for each variable: These integrating factors will not allow us to derive explicit solution as some of them are defined through the unknown variables.It is important to note that the factors are strictly positive and allow us to rewrite the system as

Macrophages
Let us rewrite Equation (A7) as Integrating (with implicit dependence on variables [µ 2 ], [I γ ], and Since [C] and [µ 1 ] are proven to remain positive, the right-hand side is bounded, hence [C] is bounded.
Appendix A.2.5.Interferon-γ and TGF-β Here we show the bound on [I γ ] and G β as we need them to prove the bound on [N].
Remark A1.Alternatively we could show a bound on [µ 1 ] and subsequent bound on Observe that on the right-hand sides of Equations (A13) and (A14), the positive terms are already proven to be bounded: Then combining these with Equations (A13) and (A14) we get the following differential inequalities: Integrating, we get which proves the bound.
Appendix A.2.6.Necrotic Cells Now we notice that for Equation (A9) all the components of the positive term are already proven to remain bounded which results in the differential inequality subsequently resulting after integration in the following bound:

Remaining Cytokines
With the bound on necrotic cells we have proven boundedness for all the positive components of the right-hand sides of the Equations (A10)-(A12): Thus, the following differential inequalities are valid: Integrating, we obtain Thus, [H], [µ 1 ] and [µ 2 ] are bounded for positive t.
Here we consider µ mean Further assumptions are based on maximal observable quantities for all the variable across all patients: See Appendix C for more details on patient data and specific values.We assume that most of T-helper cells are activated by antigen-presenting dendritic cells, so we take We also assume that inhibition of T-helper cells by µ 2 family of cytokines and by Treg cells are each 20 times more effective than the natural degradation: For cytotoxic T-cells, we assume that activation by T-helper cells is twice as effective as activation by Dendritic cells and same as for T-helper cells inhibition of cytotoxic T-cells by µ 2 family of cytokines and by Treg cells are each 20 times more effective than the natural degradation: Next assumption is that activation of Treg cells by T-helper cells is four times larger than activation by µ 2 family of cytokines and by TGF-β: while inhibition of Treg cells by µ 1 family of cytokines is 20 times larger than their natural degradation rate: We impose that activation of dendritic cells by HMGB1 is twice more effective than activation by cancer cells, inhibition of dendritic cells by HMGB1 is twice less effective than inhibiiton by cancer cells, and cumulative inhibition of dendritic cells by HMGB1 and cancer cells is equivalent to the natural degradation rate of dendritic cells: For macrophages, we assume that most macrophages are activated by T-helper cells, thus Next, we look at the cancer death rates.We assume TGF-β and IFN-γ are equally effective in killing cancer cells, but cytotoxic T-cells to be twice more effective, so where the nondimensional parameters can be expressed as follows: Cancer proliferation rate λ C and all the innate degradation/death rates remain unscaled.
Then the equations for doubling rate (A18) and (A19) become and the system of restrictions (A20)-(A30) in dimensionless form can be rewritten as These 29 restriction, together with 14 equations from requiring the steady-state of (A31)-(A44), and 11 given decay rates (A15)-(A17), scaling coefficients and given α NC , C 0 and M 0 are enough to derive all 63 non-dimensional coefficients of (A31)-(A44) from

Appendix C. Patient Data and Initial Conditions
Here we describe the processing of the data to be used for parameter estimation and initial conditions.The clustered deconvolution data, described in Section 2.3, and original TCGA data are used to calculate the immune variables as described in Table A4.TGFB1, TGFB1I1, TGFB2, TGFB3, TGFBI size(P) multiply_dimensions Total_Immune_Density T N , T h , T C , T r , D N , D, M 0 For variables related to immune cells, we substitute zero values with 10% of the smallest positive cell density for numerical stability.
We estimate the relative amount of cancer cells and necrotic cells as follows: we start by assuming that on average across all patients the ratio of immune cells:cancer cells:necrotic cells will be approximately 0.4:0.4:0.2 with variability between clusters based on tumor size.For patient P, we consider tumor size (size(P)) to be the product of the longest dimension and the shortest dimension.We assume total cell density at the steady-state to be proportional to this product as Total_Cell_Density P = α dim size(P) 1 K ∑ all P size(P) .
Tumor deconvolution data only provide us with ratios of immune cells relative to each other.Thus, to properly adjust the scaling, we take each immune cell value from deconvolution multiplied by 0.4α dim , and consider 0.4α dim ∑ (Immune cell ratios) as total immune density and Next, for each cluster, we divide patients into three groups according to their their tumor size: above average, below average and no data.Resulting patient numbers of each group are given in Table A5.We use the means from the group "above average" as steady-state assumptions given in Table 2.The data in the "below average" group as evaluated by (A45) may contain negative values for cancer and necrotic cells.The data in "no data" group have no values for cancer and necrotic cells.Thus, we substitute all non-positive and absent cancer values with 10% of the smallest positive cancer density value.We substitute all non-positive and absent necrotic cell values with zero.These changes violate the 0.4:0.4:0.2 ratio of immune cells:cancer cells:necrotic cells, and the updated ratio is 0.4475:0.3684:0.1841.
The steady-state assumptions (see Appendix B) are partially based on maximum values of each variable in the ODE system (A1)-(A14) across all patients, as well as mean value of variables T C , µ 1 , I γ and G β across all patients.The corresponding values are given in Table A6.In each cluster, a patient with the smallest known tumor size is used as initial condition (given in Table 3) for the dynamics computations presented in Figures 5 and 6.However, it is interesting to look at the dynamics of one cluster with initial conditions from another. Figure A1 shows the dynamics of each cluster with every initial condition from all clusters.It demonstrates that the curves originating from different initial conditions quickly converge to the same dynamics, more influenced by the parameters and not the initial conditions.Additionally, any patient in the "below average" group can be reasonably used as initial condition.The resulting dynamics is given by cluster on Figures A2-A6.The colors on the plot correspond to tumor size categories introduced on Figure 3F.These results suggest that the differences in the dynamics within the same cluster are mostly due to the variations in the initial tumor size.

T
h , TC , M T h , M , D Tr, M , D T h , TC , Tr M , N M , Tr

Figure 1 :
Figure 1: Network of cells and cytokines.Sharp arrows indicate activation or proliferation, and the blucked arrow indicates killing

Figure 1 .
Figure 1.Network of cells and cytokines.Sharp arrows indicate activation or proliferation, and the blunt arrow indicates inhibitions.

Figure 2 .
Figure 2. Immune cell fractions.(A) The fraction of immune cells in each colonic tumor.(B) The frequencies of immune cell types in each cluster of patients.Clusters were formed based on variations in 22 immune cell types, some of which were later combined and others that were not included in the model.Cell frequencies in this figure are averaged within the cluster.The vertical bars show the standard deviations.

Figure 3 .
Figure 3. Clinical features of the clusters.(A-C) The percentage of patients alive or dead at the last time of follow up (A), the stage of tumors I-IV at time of initial diagnosis (B) and tumor status (with tumors or without tumors) at the last time of follow up (C).(D-F) Tumor dimension (D), density of cancer cells (E) and ratio of cancer to total immune cells (F) in each cluster, respectively.Colors in (F) show the different tumor dimensions, grouped into six categories (cm 2 ): S1: 0-0.25, S2: 0.25-0.5,M1: 0.5-0.75,M2: 0.75-1, L1: 1-1.25, L2: >1.25.

Figure 4 .
Figure 4. Sensitivity analysis.The first, second and third columns of (A) respectively present the results of non-dimensional sensitivity of cancer cell density, total cell density and minimal eigenvalue of the Jacobian of the system at the steady-state.Minimal eigenvalue is used as a measure of how fast the system converges to the steady-state.(B) The sensitive parameters related to immune cells.Each row of plots shows the most sensitive parameters for each cluster of patients.

Figure 5 .
Figure 5. Cells' dynamics in colonic tumors.Time evolution of cells' density (cell/cm 3 ) for each cell type in the model and total cell density.Different colors represent the models derived for different clusters of patients and shaded regions represent the 10% variation in the most sensitive parameters.

Figure 6 .
Figure 6.Cytokines' dynamics in colonic tumors.Time evolution of RNA-seq expression rate of cytokines.Different colors represent the models derived from different clusters of patients and shaded regions represent the 10% variation in the most sensitive parameters.
The right-hand side function is bounded for positive [µ 2 ], [I γ ] and [T h ], and thus proves the bound on [M].Appendix A.2.4.Cancer Cells Rewriting the Equation (A8) as d values of the corresponding variable across all patients.

Figure A1 .
Figure A1.Cross-cluster dynamics.(A) Dynamics with parameters from cluster 1 and initial conditions from all clusters.(B) Dynamics with parameters from cluster 2 and initial conditions from all clusters.(C) Dynamics with parameters from cluster 3 and initial conditions from all clusters.(D) Dynamics with parameters from cluster 4 and initial conditions from all clusters.(E) Dynamics with parameters from cluster 5 and initial conditions from all clusters.(F) Initial dynamics with parameters from cluster 5 and initial conditions from all clusters on a time interval [0, 300].

Figure A2 .
Figure A2.Different initial conditions for cluster 1.Based on patients in the "below average" category.Colors correspond to tumor size categories introduced in Figure 3F.

Figure A3. Different initial conditions for cluster 2 .
Figure A3.Different initial conditions for cluster 2. Based on patients in the "below average" category.Colors correspond to tumor size categories introduced in Figure 3F.

Figure A4 .
Figure A4.Different initial conditions for cluster 3. Based on patients in the "below average" category.Colors correspond to tumor size categories introduced in Figure 3F.

Figure A5. Different initial conditions for cluster 4 .
Figure A5.Different initial conditions for cluster 4. Based on patients in the "below average" category.Colors correspond to tumor size categories introduced in Figure 3F.

Figure A6 .
Figure A6.Different initial conditions for cluster 5. Based on patients in the "below average" category.Colors correspond to tumor size categories introduced in Figure 3F.

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

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

Table A4 . Patient data correspondence with variables.
Correspondence between the system variables and the source data from TCGA and deconvolution results.

Table A5 . Distribution of patients according to their tumor size.
Evaluated relative to the average tumor size within each cluster.