Patient-Specific Mathematical Model of the Clear Cell Renal Cell Carcinoma Microenvironment

The interactions between cells and molecules in the tumor microenvironment can give insight into the initiation and progression of tumors and their optimal treatment options. In this paper, we developed an ordinary differential equation (ODE) mathematical model of the interaction network of key players in the clear cell renal cell carcinoma (ccRCC) microenvironment. We then performed a global gradient-based sensitivity analysis to investigate the effects of the most sensitive parameters of the model on the number of cancer cells. The results indicate that parameters related to IL-6 have high a impact on cancer cell growth, such that decreasing the level of IL-6 can remarkably slow the tumor’s growth.


Introduction
Clear cell renal cell carcinoma (ccRCC) is a type of kidney cancer that makes up 80% of kidney cancer cases [1]. Surgery is the most common treatment for ccRCC, although it is not effective when the tumor metastasizes to other parts of the body [1][2][3]. Other treatment options, such as immunotherapy and targeted therapy, are used to regulate the growth of the tumor, and radiation therapy, thermal ablation, and cryosurgery are used to kill cancer cells [1]. However, not all patients respond well to these treatments, leading to more complications [1,2].
Mathematical models can provide valuable information about the complex interactions in biological systems, such as cancer [4][5][6][7][8][9][10][11][12][13]. For example, Pillis et al. created a mathematical model on the effects of regulatory T-cells (Tregs) on ccRCC treatment via a drug called sunitinib, which resulted in an improvement in controlling cancer in 40% of patients compared to patients without treatment [14]. Since vascular endothelial growth factor (VEGF) helps in angiogenesis in tumors, another study by Sharma et al. studied the effectiveness of the VEGF receptor 2 (VEGFR2) inhibitor in renal cell carcinoma (RCC) patients and found that an inhibitor compound named SCHEMBL469307 is the most effective [15]. In addition, complex data-driven mathematical models of other cancers, such as breast cancer, osteosarcoma, and colon cancer, were developed by Mohammad Mirzaei et al., Kirshtein et al., and Le et al. to investigate the interactions between the important immune cells and cytokines involving many key cells and molecules corresponding to the particular type of cancer and found that the interactions between the cells and molecules are important to understand the dynamics of cancer growth [16][17][18].
The ccRCC microenvironment has the highest levels of immune cell infiltration among all epithelial cancer types, making it a pro-inflammatory environment [36].
Infiltrating T-cells, such as cytotoxic and helper T-cells, play essential roles in controlling the tumor by targeting antigenic tumor cells [37]; thus, treatments (such as sunitinib and pazopanib), which directly target these T-cells in ccRCC, have been dominant for a long time. However, it has been observed that tumor microenvironment heterogeneity can significantly affect the outcomes of these treatments. A study on patients in the phase III trial receiving sunitinib or pazopanib shows poor outcomes for a cluster of patients with a high immune infiltration (especially macrophages) and significantly higher PD-L1 expression on tumor cells compared to the other clusters [38]. Another study based on VEGFR TKI therapy on 53 metastatic ccRCC patients revealed undesirable outcomes in a group of patients whose tumors had an intense Th1-oriented inflammatory and suppressive immune environment, with high levels of PD-1/PD-L1 [39]. These observations emphasize the sensitivity of clinical treatment responses to tumors' immune profiles.
This paper investigates the dynamics of tumors based on their immune patterns by developing an ODE model that considers the interaction network of key players that we mentioned above (see Figure 1). We clustered the tumors based on their immune patterns and compares their dynamics. More importantly, we found the most sensitive parameters for each cluster of tumors. The parameters in this model were estimated using the steady state assumptions and data acquired from gene expression data from The Cancer Genome Atlas (TCGA), separately for each cluster of tumors. The result was a data-driven patientspecific ODE model aimed at understanding the impacts of immune patterns and their interaction network in tumor progression. Moreover, regarding its capability to extend to a treatment model, one can treat the immune cells as boundary variables and integrate this model as a compartment for larger-scale models. Furthermore, by including metabolites, such as nutrients and oxygen, one can extend the model to incorporate angiogenesis. However, at this stage, we avoid introducing more complexity to this model.

Variables and Network
As mentioned in the introduction, ccRCC has the highest infiltration of cytotoxic cells, dendritic cells, Th1, and macrophages among all of the other epithelial cancers, and low levels of regulatory T-cells and Th2 [36]; the response to treatments are associated with the percentage of these cells in the tumors. Therefore, including Th1 and Th2 together as helper T-cells, we used these key cells in addition to cancer and necrotic cells, and the most important cytokines secreted by them as our model variables.
In addition, both cancer and necrotic cells release relatively large portions of HMGB1 that helps promote dendritic cells and helper T-cells [30,31,33,51]. Dendritic cells then promote both helper T-cells and cytotoxic cells [19,30,31,51,52]. To demonstrate the interactions among the immune cells, cytokines, cancer, and necrotic cells, we developed an interaction network in Figure 1 and provide the list of variables and how we derived them in Table 1.

Mathematical Model
We derived an ordinary differential equation (ODE) for each variable in the model using λ to denote the growth/promotion rates and δ to denote the decay/inhibition rates of the associated cells and molecules. In the case of promotion and inhibition, the second subscript refers to the promoter or inhibitor, and the first subscript is the cell or molecule subjected to promotion or inhibition.

Helper T-cells (T h ):
Helper T-cells are activated by macrophages, dendritic cells, HMGB1, and IL-12 [34,[53][54][55]. All of the T-cells in the model differentiate from naive T-cells. However, this differentiation happens mostly outside of the microenvironment. For simplicity, we added naive T-cells to the model to prevent blow-ups in the population of active T-cells without having to introduce nonlinearity in their corresponding ODEs [16,56].
Regulatory T-cells inhibit helper T-cells, and IL-10 [45,46]. Therefore, the dynamics of helper T-cells are modeled by the following equation.
Regulatory T-cells (T r ): IL-12 and IL-2 activate regulatory T-cells [28,41,43,55,63]. We only assume natural death for the regulatory T-cells since the model has no major inhibitors for this cell type. Thus, we model the dynamics of T-reg cells by Naive T-cells (T N ): As mentioned earlier, naive T-cells are not a part of the tumor microenvironment. They differentiate into other T-cell types mostly before the microenvironment immune infiltration [34,56]. We modeled naive T-cells starting with a constant growth rate and then deducted the differentiation rate of other T-cells.
Naive macrophages (M N ): Since macrophages are derived from naive macrophages [44], we used an approach similar to the naive T-cells to model naive macrophages.
Naive dendritic cells (D N ): Mature dendritic cells are derived from naive dendritic cells [68][69][70]. So, we modeled naive dendritic cells similar to naive T-cells and naive macrophages.
Cancer cells (C): Although cancer cells grow rapidly, their growth can be affected by the lack of space or nutrients. Therefore, we added a logistic model term C 0 ) to control the population of cancer cells corresponding to a carrying capacity C 0 . Cancer cells are also activated by IL-6 [49,71] and deactivated by CD8+ T-cells, NK cells [47,50,56], and interferon-γ [29,61]. Moreover, cytotoxic cell apoptosis happens when the programmed cell death proteins (PD-1 and PD-2) expressed by CD8+ T-cells attach to the programmed cell death ligands (PD-L1 and PD-L2) expressed by cancer cells [26,27]. This way, cancer cells can avoid CD8+ T-cell cytotoxicity. We combined the concentrations of PD-1 and PD-2 and denote them by [PD]; we denote PD-L1 and PD-L2 together as [PDL]. We make their concentrations proportional to the cells they are expressed by. So, we assume that the concentrations of PD-1/PD-2 and PDL-2 are proportional to cytotoxic and cancer cells, respectively.
We model the cancer cell population by the following.
Necrotic cells (N): Necrosis is a process in which the dead cells are not cleared out of the body or the system, unlike apoptosis, autophagy, etc. [30,33]. In a tumor site, necrosis happens mostly through cancer cell death; thus, we modeled the necrotic cell populations by taking a proportion of cancer cell death and a much smaller decay rate [35].
Hence, we have 15 ODEs in the system and 67 parameter values to determine to solve the system.

Data Preparation
We used gene expression data of ccRCC from The Cancer Genome Atlas (TCGA) along with the immune classification by Su et al. [2] where the TCGA gene expression was normalized, and CIBERSORTx B-mode was used to derive the cell fractions [73]. The patient classification was conducted using unsupervised K-means clustering, resulting in four groups based on each cell type proportion, such as T-cells, B-cells, macrophages, dendritic cells, etc. [2]. Due to the lack of time course data, we represent each cluster of patients as one virtual patient as their immune patterns are similar. We consider the smallest tumor in each cluster to represent the first stage and the largest tumor to represent the last stage of progression over time. After clustering, the cell proportions from the CIBERSORTx result and the protein and molecule concentrations from the normalized gene expression data were derived from [2], according to the variable combination in the model described in Table 1. We only considered patient data with less than 0.5 p-values from the CIBERSORTx result. The cell frequencies derived from CIBERSORTx for each cluster are given in Figure 2, produced by TumorDecon software [74]. We determined the actual cell population based on the intermediate tumor size described in the TCGA clinical data of renal cell carcinoma patients [73]. We chose the ratios of immune cells to cancer cells to necrotic cells as 0.3:0.6:0.1, as was done in [18]. We assume that the epithelial cell density is 4.5 × 10 4 cell/cm 3 in the cancer microenvironment [16][17][18]75]. We also let the average cell density scale be α = 4.5 × 10 4 across all ccRCC patients. We calculate the total cell number (TCN) of each patient P by where K = 346 is the total number of patients having primary tumors from the original data. Then the total immune cell (TIC) population is calculated by, We let the cells and cytokines of the smallest tumor in each cluster be the initial conditions of the system and nondimensionalized the variables by dividing them by their steady state values. The nondimensional initial conditions are given in Table 2 and the steady state values are given in Table 3.

Parameter Estimation
We have 15 ODEs and 67 parameters, but there is not enough biological information to estimate all of the parameter values specific to ccRCC. However, we collected and derived some of the decay rates, such as δ T h , δ T c , δ T r , δ T N , δ M , δ D , δ I γ , δ H , δ IL 10 , δ I 2 , and δ IL 6 by using where δ X is the decay rate of the corresponding variable X and t X 1/2 is its half-life [16][17][18]. We took the average for the variables referring to combined quantities and then computed δ X . For instance, the half-life of IL-2 is approximately 7 min [76] and the half-life of IL-12 is almost 3.6 h [16], so we took the average of these half-lives to form the half-life of I 2 (i.e., t I 2 1/2 = 7.743 × 10 −2 days), which gives us δ I 2 = 2.238 molecules day . The following are the decay rate values that were either collected from [16,18] or calculated based on half-lives of the variables, We estimate the remaining parameters of the model by assuming the large tumors are at a steady state and utilize their data. The ODEs in Section 2.2 provide 15 algebraic equations with 67 unknowns at the steady state. Given the 11 decay rates determined above, we need to find the values for 41 unknown parameters. This system is heavily under-determined, and it is impossible to extract a unique parameter set in its current state. However, we can add some biologically feasible mathematical assumptions to remedy this issue. These assumptions simply describe a relationship between activation or inhibition of a certain cell or cytokine by other cells or molecules. This will create a system of algebraic equations at the steady state, uniquely solvable for the parameter values. See Appendix A.1 for more details on this process. We included the nondimensional parameter values in Table A3. Furthermore, given that our extra assumptions are mostly of a mathematical nature, we will assess their effects on the model dynamics by scaling them.

Sensitivity Analysis
Since there is not enough biological evidence to estimate all of the parameter values of the ccRCC model, the limitations of unknown parameter estimations must be considered when using the result of the dynamics. We performed a global gradient-based sensitivity analysis to assess our estimations by scaling the 34 parameter assumptions in Appendix A.1. We carried out 5000 scalings for each parameter, leaving us with 34 × 5000 = 170,000 variations of parameters, which is significant but still a limited number. When performing the sensitivity analysis, we used the nondimensionalized system explained in Appendix A.2 to keep the computations stable. The sensitivity level of an ODE system, such as dX dt = F(X,θ, t), whereθ = θ 1 , θ 2 , . . . , θ N represents the parameter vector, is calculated by, where X * is the solution of the system at the steady state. In this paper, we calculated the sensitivity of cancer and total cells to all the parameters at their steady states. For each variable X * , we obtained the sensitivity vector for the system by differentiating it with respect to θ i and setting F(X * ,θ) = 0. We obtain the formula where ∇F(X,θ) −1 is the numerically approximated inverse Jacobian of F with respect to X [18]. Then, we followed the methodology for global sensitivity given in [18].

Dynamics
The dynamics of the cells and molecules over 5000 days are presented in Figure 3.
The helper T-cell population decreases at the beginning and eventually increases to reach a steady state in clusters 3 and 4. In cluster 2, it increases, and in cluster 1, it decreases to reach a steady state within the first few days. As a result, the helper T-cell populations in clusters 1 and 2 remain constant almost all of the time.
Cytotoxic cells in all clusters decrease in the first few days but eventually increase in population to reach a steady state. Furthermore, cytotoxic cells in cluster 2 start with the highest population and achieve the highest saturation level. This is followed by cluster 4, with the second highest saturation level, where the cells grow faster and reach a steady state in around 1500 days. In cluster 3, cytotoxic cells reach the lowest saturation level, leaving cluster 1 with the second lowest saturation level. Overall, the steady state populations and the time points that each cluster reaches steady states are noticeably different in each cluster.
T-reg cells in clusters 1, 2, and 4 increase quickly from their initial conditions and decrease to reach a steady state. T-reg cells in cluster 4 reach a steady state at around 1500 days and clusters 1 and 2 reach a steady state at around 3000 days. In cluster 3, T-reg cells decrease and reach a steady state faster than in any other cluster. Furthermore, the slow increase in cytotoxic cells in clusters 1 and 2 could be due to regulations by T-reg cells that slowly decrease in these clusters. Similarly, the faster growth in cluster 4 of cytotoxic cells could be related to the fast decrease of T-reg cells in cluster 4.   Naive T-cells in cluster 1 attain the highest steady state population, followed by 3, 2, and 4, respectively. In clusters 1 and 2, naive T-cells increase slowly and reach a steady state of around 2000 days, whereas clusters 3 and 4 stabilize before and after 1000 days.
Macrophages (Mφ) in clusters 1 and 3 initially increase and then soon start to decrease to reach a steady state, but in clusters 2 and 4, they decrease to reach a steady state. Naive macrophage growth in all clusters does almost the opposite of macrophages, which could be because macrophages are derived from naive macrophages. The overall trend for macrophages is that they decrease to reach a steady state in all clusters, while naive macrophages increase to reach a steady state.
The overall trend for both mature and naive dendritic cells is to decrease to reach a steady state. However, dendritic cells in cluster 3 initially increase quickly and drastically and then suddenly decrease. The same happens to naive dendritic cells in cluster 1. Moreover, cluster 1 stands out in naive dendritic cells by achieving comparatively higher steady state values than the others.
Cancer and necrotic cells exhibit exponential growth until they reach a steady state and have similar curves, as necrotic cells are produced at a rate proportional to cancer cell decay. Cancer cells in cluster 1 attain the highest steady state population. Cluster 3 stands out by growing the fastest and requiring less time to reach a steady state, which could be related to its low level of cytotoxic cells (especially since cytotoxic cells also reach a steady state at around 1000 days in this cluster). Cluster 2 has the slowest growth despite having the highest initial population. Clusters 2, 3, and 4 all achieve a similar steady state population sooner or later within 5000 days. The slow growth in cluster 2 can also be attributed to the significantly high cytotoxic levels. Overall, we can see a clear correlation between cancer progression and cytotoxicity levels.
Interferon-γ initially increases but starts to decrease to reach a steady state in all clusters. HMGB1 concentration drops at the beginning and then increases to reach a similar steady state in all clusters. Cluster 1, as with cancer and necrotic cells, achieves the highest steady state concentration in HMGB1. It is consistent with our assumptions of parameters, as HMGB1 is mainly secreted by cancer and necrotic cells in the cancer microenvironment. IL-10 concentrations in clusters 1, 3, and 4 increase initially but decrease within a few days. However, cluster 3 concentration of IL-10 increases rapidly to a higher concentration than other clusters and decreases to reach a similar steady state concentration. The cluster 1 concentration of IL-10 reaches a steady state later than other clusters (around 700 days). Cluster 2's IL-10 concentration decreases to reach the lowest steady state concentration among all clusters. IL-2 and IL-12 concentrations in clusters 1, 3, and 4 also increase at the beginning, then decrease to reach a steady state. As with IL-10, cluster 3 concentration I 2 attains its maximum concentration very fast and decreases to reach a steady state concentration. In cluster 1, it follows a similar trend as IL-10 and slowly reaches a steady state at around 2000 days. Finally, the concentration of IL-6 in clusters 1, 3, and 4 increases and then decreases to reach a steady state. However, the concentration of IL-6 in cluster 2 remains constant for the most part.
Since we assume that the concentrations of PD − 1/2 are proportional to cytotoxic cells and the concentration of PDL − 2 is proportional to cancer cells, we refrain from including their dynamics.

Sensitivity
The sensitivity analysis reveals that the parameters directly involved in the cancer ODE are the most sensitive parameters for cancer cells and total cells. Additionally, macrophagerelated parameters, such as λ MT h , λ MIL 10 , and λ MI γ , and T-reg cell-related parameters, such as λ T r D and λ T r I 2 , and an IL-6-related parameter, λ IL 6 D , also show significant sensitivity values. The sensitivity plots for the most varying parameters on cancer cells and total cells are shown in Figures 4 and 5.   We could say that a quantity is positively sensitive to a parameter when an increase in the parameter causes the quantity to become larger. A quantity is negatively sensitive to a parameter when an increase in the parameter causes a decrease in the quantity. According to the sensitivity plot in Figure 4, cancer cells in all clusters are negatively sensitive to their decay rates and then positively sensitive to their growth rates. For all clusters, the most sensitive parameters for total cells are the decay/inhibition rates of cancer and macrophages and cancer cell growth/promotion rates.
In addition, Figure 5 illustrates more sensitive parameters for cancer cells and total cells. For sensitivity to cancer and total cells, we see a significant emphasis on macrophagerelated parameters in all clusters, especially clusters 1 and 2. Additionally, regulatory T-cell promotion and decay rates play more significant parts in the cancer sensitivity results for clusters 3 and 4. Finally, we see λ IL 6 D as the sixth sensitive parameter for cancer in cluster 1.
The sensitivity analyses of the parameters show that the direct and indirect interactions between the variables impact cancer or total cell growth in the cancer microenvironment. For instance, the inhibition of cancer cells by cytotoxic T-cells, IFN-γ, or natural decay, and its promotion by IL 6 and natural growth, are direct pathways that increase or decrease the number of cancer cells. On the other hand, we notice that cancer cells are negatively sensitive to the decay rate of macrophages ( Figure 4) and positively sensitive to their promotion rates ( Figure 5). This is because the macrophages secrete IL-6, and then IL-6 helps promote cancer. So, decay in macrophages would lead to less secretion of IL-6, leaving cancer cells with fewer resources and vice versa. Moreover, T-reg cell parameters play an important role in cancer development by negatively impacting cancer with their growth and positively affecting cancer with their decay, especially in clusters 3

Varying Dynamics of Cancer Cells with Scaled Assumptions
Appendix A.1 shows that most of the parameters (except possibly some decay rates), including the sensitive ones, were derived from restrictive assumptions. To assess the validity of these assumptions, we scaled them using factors of 1, 0.2, and 5. These scalings created a new set of parameters with significantly different values. Then from these sets, we perturbed the most sensitive parameters illustrated in Figures 4 and 5 by 5% to create an interval of confidence for the new dynamics. After scaling all the assumptions in Appendix A.1, only 3 caused significant changes to cancer dynamics. These assumptions are as follows.
These three scalings cause significant changes in the parameters as illustrated in Figure 6. We can see that assumption (20) causes a significant deviation from the original values for parameters A D N , δ DC , λ DC , and λ DH . Moreover, scaling assumptions (21) and (22) causes parameters, such as λ IL 6 D , λ IL 6 T h , λ IL 6 M , and λ IL 6 C to change but not as drastically as the changes imposed by assumption (20). What is interesting is that none of these parameters (except for λ IL 6 D in cluster 1) are among the sensitive parameters. However, the changes caused by the scaling of the assumptions were so large that the effects were tangible. We emphasize that none of the other assumptions left such an impact after scaling. Figure 7 shows the cancer dynamics with the original parameter values next to the dynamics acquired by scaling the assumptions (20)-(22) by 0.2 and 5, respectively. The shaded regions are the regions of confidence acquired from perturbing the most sensitive parameters in Figures 4 and 5 by 5%. We see slight changes in clusters 2, 3, and 4, and a more significant change in cluster 1. Moreover, the width of the shaded regions remained the same in all cases. We saw that by changing the assumptions (20)-(22), significant deviations occurred to the parameter values mainly involved in the dendritic cells and IL 6 production. Cluster 1, which was more impacted, was the only cluster sensitive to the parameter λ IL 6 D . Even though the assumptions are mostly modeling artifacts for parameter estimations (and one must be cautious when using them), these results suggest interesting control potentials for IL 6 .
Although the parameters α T c and β C were not among the most sensitive parameters, to make sure that the assumptions on these parameters were reasonable, we scaled α T c with 0.2 and 5 with 5% variations to the most sensitive parameters to see their impacts on cancer cells (Figure 8). Since α T c is multiplied by β C in the ODE system, we did not scale it as our goal was to keep the scaling of the term that involved these parameters similar to α T c . Thus, the scaled equation for α T C becomes, As a result, the term involving α T C in the cancer cell equation becomes,

Discussion
Cancer is a diverse and complex disease that cannot be understood with a unified approach. Each cancer patient is unique in his/her immune infiltration profile and cancer type. In this paper, we proposed a model describing the critical interactions in the ccRCC tumor microenvironment. The model was inspired by the emerging clinical model of targeted therapy efficacy within distinct ccRCC subgroups. These models consider the infiltrating immune cells, such as cytotoxic, helper, and regulatory T-cells, macrophages, and dendritic cells (see Figure 1 in [36]). Therefore, we clustered four immunologically different virtual patients and observed the development of their tumors through their tumor microenvironment interactions. One of the main challenges in this study is the shortage of biological data for estimating the parameters in the model. So we added some mathematical (yet biologically reasonable) steady state assumptions to circumvent this issue. We then assessed the validity of these assumptions and the acquired parameter values in two ways. (1) By performing a global sensitivity analysis and observing the effects of varying the most sensitive parameters on cancer dynamics, and (2) by scaling our assumptions and studying their effects on cancer dynamics.
Despite the limitations of the model, we could infer interesting pathways in ccRCC and its microenvironment development. For instance, pathways that directly or indirectly affect the cytotoxicity of the microenvironment seemed to affect the cancer dynamics. We saw that a higher population of cytotoxic cells in cluster 2 is correlated to slower cancer cell growth. This could indicate that environments with initially higher cytotoxicity have a better shot at controlling ccRCC. Furthermore, the slow cancer growth in cluster 1 at the beginning can be related to naive dendritic cell growth patterns. More naive dendritic cells result in more mature dendritic cells that promote helper T-cells, aiding in more cytotoxicity in the environment [68]. Although cytotoxic cells were higher in population cluster 4, macrophages started to deplete soon for this cohort of patients, which caused faster growth in cancer cells via helper T-cell and cytotoxic cell pathways. Similar behavior can be justified for cluster 3. Another important relationship we observed in cluster 3 was that the dynamics of dendritic cells, INF-γ, IL-10, IL-6, and IL-2 and IL-12 looked very similar. Dendritic cells secreted all of the cytokines mentioned above, indicating that they played an important role in leading the dynamics of these cytokines in cluster 3. While increased IFN-γ may inhibit cancer cells or IL-2 and IL-12 may promote cytotoxicity, the inhibition of cytotoxic cells and helper T-cells by IL-10 or the promotion of cancer cells by IL-6 may have impacted more in aiding cancer cells to grow very fast in cluster 3. Therefore, although cancer cells are not directly impacted by macrophages or dendritic cells, cancer cell interaction with the cells and cytokines that macrophages and dendritic cells promote reveal potential pathways for better understanding the ccRCC microenvironment.
We performed a global sensitivity analysis on the parameters of the model to see their impacts on cancer and total cell growth. We saw that the most sensitive parameters were the production/promotion and decay/inhibition rates directly involved in cancer ODE, which was expected. In addition, we noticed a lot of macrophage-related production and inhibition rates as the second most sensitive parameters. Specifically, the sensitivity values hinted that more macrophages lead to worse prognoses. It has been shown that macrophage polarization into anti-and pro-tumor subtypes can have different implications for cancer progression and prognosis [77,78]. The patient data and the parameter sensitivity analysis of the model reveal that the macrophages act as pro-tumors in all patients, suggesting them as targets to suppress via therapy. Additionally, cytokines and cells act as promoters or inhibitors in sensitive parameters, such as IFN-γ, CD8+ T-cells, IL 6 , IL 10 , regulatory T-cells, dendritic cells, IL 2 , and IL 12 can be potential targets. We also assessed the validity of our assumptions by scaling them. Almost all the assumptions led to no significant changes in the dynamics of cancer cells, except for three. These three assumptions, directly or indirectly, contributed to the significant production or inhibition of IL 6 . Numerous studies focus on controlling the IL-6 as a therapeutic option to control the progression of ccRCC [79][80][81][82].
As with all the other mathematical models, our model is not without limitations. Many assumptions were used here, even though biologically feasible, but they were of mere mathematical nature. Some examples of such assumptions are: (1) the estimation of initial conditions by the small tumors in each cluster and assuming the largest tumors to give us the steady state values, and (2) the specific scales used to enforce the dominant activator or inhibitor of certain cells or cytokines. We tried our best to assess the sensitivity of our model to these assumptions and justify their robustness. However, none of these attempts can replace a direct validation using time-course data, and we acknowledge this limitation of the study. Moreover, for simplicity, we ignored important factors, such as angiogenesis and metabolites. We acknowledge this as another limitation of our model. However, this model provides a basis for future computational models of the ccRCC tumor microenvironment and can be extended to include new biological and biomedical discoveries and find optimal treatment options for ccRCC patients. https://events.cancer.gov/cbiit/dtwin2020 (accessed on 12 August 2020). We also thank Ana Belen Moscoso Gonzales for reading the paper and for some editorial suggestions.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:  [83]. Thus, we considered the doubling rate to be the difference between the growth and decay rate of cancer cells, as was done in [16][17][18]. The faster doubling rate includes the rate of cancer cells' natural proliferation and proliferation by IL-6, resulting in Then, we consider the slower doubling rate to be only the natural cancer cell proliferation while the decay rates of cancer cells include the decay rates by all the inhibitors, including cancer's natural decay rate, The mean values of the variables used in Equations (A1) and (A2) are given in Table A1. Since we assume the parameter relations based on the steady states of each variable, the parameters should all be positive.
The steady states are derived from the largest tumor across all patients. See Table A2. We assume that helper T-cells are mainly activated by macrophages and dendritic cells, Moreover, for helper T-cells, we assume that the natural inhibition is smaller than inhibition by T-reg cells and IL-10, so that For cytotoxic cells, we assume that the activation or proliferation is mainly done by helper T-cells, dendritic cells, and IL-2 since IL-2 has been described as an effective media to treat cancer [28,76], and we assume that the inhibitions are mainly done by T-reg cells and IL-10.
For T-reg cells, we assume that dendritic cells mainly activate them, For macrophages, we assume proliferation scales of helper T-cells, IFN-γ, and IL-10 are similar, so we have, We further assume that the decay rates of naive macrophages are similar to the decay rates of macrophages, Similarly, we assume that the proliferation rates for dendritic cells are similar, and their inhibition rate by cancer cells is larger than their natural decay rate. Moreover, we take the decay rate of naive dendritic cells to be the same as that of dendritic cells.
We calculate the parameter α T c based on the average fraction of PD − 1 and PD − 2 among all the genes and then multiply it by the maximum of cytotoxic cells. Similarly, we calculate β C based on the average of PDL − 2 among all genes(as it was the only available ligand in the data) and then multiply it with the cancer cell maximum. Since cancer cells grow rapidly, we let and For necrotic cells, we assume that We assume that IFN-γ is secreted mostly by cytotoxic cells, Since HMGB1 is mainly produced by cancer cells and necrotic cells, we let the HMGB1 production rate by cancer and necrotic be larger, For IL-10, we assume that it is secreted on a similar scale by helper T-cells, T-reg cells, cytotoxic cells, dendritic cells, and macrophages; Since helper T-cells mainly produce IL-2, we let the rate of IL-2 secretion by helper T-cells be larger [28]. We then let other cell production rates of IL-12 be similar.
Lastly, since IL-6 is mainly produced by cancer cells, we let

. Nondimensionalization of Parameters and Variables
First, consider the variables at a steady state as, To avoid complexity in parameter estimations, numerically solving the ODEs, and analyzing the sensitivity of parameters, we nondimensionalize the variables and parameters based on their steady state values. We denote the nondimensional variables by [X] = [X] [ where X represents each variable and [X] = 1 when it is in a steady state. However, we do not nondimensionalize the time variable since it does not introduce any complexity to the problem. Thus, the system consisting of Equations (1)-(16) (except (9)) would be nondimensionalized, as follows: Since the dimension of the time variable remains unchanged, the decay rates in Table A3 also remain unchanged. We nondimensionalize the assumptions in Appendix A.1 as follows: In addition, we nondimensionalize independent production rates, such as A T N , A M N , and A D by , and A D [D ∞ ] , respectively. According to the above nondimensionalization of the system and assumptions, we have the nondimensional parameters in Table A3. δ T N 9.49 · 10 −4 9.49 · 10 −4 9.49 · 10 −4 9.49 · 10 −4 λ MT h 6.60 · 10 −3 1.20 · 10 −2 1.33 · 10 −2 1.56 · 10 −2 λ MI γ 6.60 · 10 −3 6.57 · 10 −3 1.09 · 10 −3 2.95 · 10 −3 λ MIL 10 6.60 · 10 −3 1.24 · 10 −3 5.37 · 10 −3 1.26 · 10 −3