In Silico Modeling of Immunotherapy and Stroma-Targeting Therapies in Human Colorectal Cancer

In recent years, tumor immunotherapy has become available to treat malignant tumors with several drugs already approved for solid tumors. Their main mode of action is the activation of the adaptive immune system via checkpoint inhibition (1). Other widely used therapeutic strategies aim at increasing the number and the reactivity of effector cells by adoptive cell transfer or vaccination (2, 3). Additionally, many more drugs are currently investigated in clinical trials, especially combination therapies aiming at the stromal compartment (4, 5). These approaches act on the complex interactions in the tumor microenvironment, as we could show recently in macrophage-targeted immunotherapies (6).

Still, for most patients with solid tumors, no effective immunotherapy strategy is available. Especially for microsatellite-stable colorectal cancer, the most common form of this disease, immunotherapy has been largely ineffective. With a 5-year survival of just 11% in metastatic diseases (7) and a huge disease burden (8), the therapeutic need is high. In this setting, the development of more complex interventions into the immune landscape requires detailed knowledge of the interactions of relevant players and their response to interventions. Although our understanding of these processes has advanced considerably in the last years, there is still no comprehensive systems perspective of all relevant interactions. To understand a complex system, it is not sufficient to have a detailed characterization of all its components. Instead, a complex system can show emergent behavior that does not arise from a specific component but from the interaction of different components.

Agent-based models are a powerful method to investigate the interactions in complex systems (9). An agent is the smallest unit in this model and can show different types of stochastic behavior, including interaction with other agents. Although these models simplify many aspects of reality, they have been shown to be extremely useful in a wide number of circumstances (10–12). In cancer research, these models are emerging as valuable tools to study emergent behavior in complex ecosystems (13), especially in stem-cell models of tumor growth (14, 15) and are used to study the mutational landscape of solid tumors (16, 17). Furthermore, they are increasingly used to optimize therapies, for example radiation therapy of solid tumors (18). Also, some models of immune-cell interactions with (19–24) or without tumor cells (25) have been proposed. Although these studies gave important insight into parts of the tumor-immune interaction, they did not accurately reproduce the diverse spatial patterns in human tumors and did not investigate therapeutic strategies.

In the present study, we generated a multiagent-based model from quantitative histological and other wet lab data, based on the concept of immune surveillance (26). We focused on parameters that could be morphologically measured and created a simplified yet powerful model of cellular interactions that shows emergent behavior. This agent-based model incorporates stochastic interactions between tumor cells, immune cells, and stroma and faithfully represents diverse spatial pattern observed in histological samples of human colorectal cancer tissue. Furthermore, the derived clinical predictions could be validated in an independent colorectal cancer cohort. This model was then used to systematically test the effect of different therapeutic interventions on this system and to create specific recommendations for effective immunotherapies.

All experiments were conducted in accordance with the Declaration of Helsinki, the International Ethical Guidelines for Biomedical Research Involving Human Subjects (CIOMS), the Belmont Report and the U.S. Common Rule. N = 20 human tissue samples of colorectal adenocarcinoma were used as a calibration cohort. These samples were retrieved from the pathology archive at UMM (University Medical Center Mannheim, Heidelberg University, Mannheim, Germany) after approval by the institutional ethics board (Ethics Board II at University Medical Center Mannheim, decision number 2017-806R-MA, granted to A. Marx and waiving the need for informed consent for this retrospective and fully anonymized analysis of archival samples; Supplementary Table S2).

Another set of 37 colorectal adenocarcinoma samples was used as a morphological validation cohort, composed as follows: N = 22 tissue samples were provided by the tissue bank of the National Center for Tumor diseases (NCT, Heidelberg, Germany) in accordance with the regulations of the tissue bank and the approval of the ethics committee of Heidelberg University (tissue bank decision number 2152, granted to N. Halama and J.N. Kather; informed consent was obtained from all patients as part of the NCT tissue bank protocol; Supplementary Table S3). Fifteen additional samples were used as described before (Supplementary Table S3; ref. 27). As the model was designed to be valid for all types of colorectal cancer tissues, we included primary tumor samples and colorectal cancer liver metastases in all cohorts (Supplementary Tables S2 and S3).

We performed histological staining for Ki67 (Dako M7240 antibody, 1:100), active caspase 3 (Abcam ab2302 antibody, 1:50), and CD3 (Leica Novocastra NCL-L-PS1 antibody, 1:50) on a Leica Bond automatic staining device using a hematoxylin–diaminobenzidine (DAB) staining protocol as described previously (6, 28). Stained whole slide tissue sections were digitized as described previously (6, 28). On histological sections, we manually identified areas homogenously occupied by tumor cells or immune cells (in the tumor or around the tumor). The fraction of Ki67-positive cells (active caspase 3–positive cells, respectively) was quantified in these regions of interest (ROI) per slide using a digital pathology approach analogous to our previously published approaches (29, 30). Definiens Tissue Studio (Definiens) was used for semiautomatic tissue segmentation and automatic cell segmentation. The fraction of positive cells was calculated as the number of positively (diaminobenzidine-positive) stained cells divided by the number of all cells in the respective ROI. On average, each ROI contained approximately 3,000 cells. Intensity thresholds for cell detection and classification were set manually and were identical for all images. The quality of cell segmentation was checked manually for each image and was found to be sufficient. We assumed that the median Ki67-positive fraction and the median active caspase 3–positive fraction approximated the probability of a given cell type to proliferate or die at one time point. All raw measurements are listed in the Supplementary Data.

For the horizontal migration experiments, Matrigel (undiluted, BD) was evenly plated at the sides of a 24-well chamber in two half-moon shapes (cat eye configuration) and the remaining third in between was filled with either pure collagen or Matrigel with CXCL9 (10 ng/mL) and CXCL10 (10 ng/mL). After gelling overnight in a humid chamber at 37°C, 5% CO₂, the well was filled with a thin layer of T-cell culture media and healthy donor T lymphocytes (CD3/CD28 activated and in an independent experiment non-activated) were placed on the right half-moon shaped third. After migration for 48 hours at 37°C, the resulting distribution in the well was documented and distances were documented and used for the multiagent model system.

Initially, our model had 22 parameters, as shown in Supplementary Table S1. These parameters were based on a clear set of assumptions. Key assumptions based on previous studies are related the tendency of immune cells to migrate towards tumor cells (31–33), lymphocyte exhaustion (23, 34), desmoplastic reaction (35, 36), and stromal permeability (37). All assumptions are listed in Table 1. Some of the model parameters had been estimated in previous studies. Specifically, this applies to the maximum proliferation capacity of non-stem cells (tumor cells and immune cells alike), which we set to 10 analogous to (14). Also, it was previously shown that the maximum number of kills a lymphocyte can deliver can be validly estimated as five (23). Of the remaining parameters, five were measured histologically: Tumor cell proliferation and apoptosis, immune cell proliferation and apoptosis, and distance to necrosis (Supplementary Fig. S1A–S1F). Thus, 12 free parameters remained and were set to biologically plausible values (Supplementary Table S1).

Of 20 tissue samples, 12 contained necrosis and the distance to necrosis from the outer tumor margin was measured at three locations, giving 36 distance values. Mean distance was 1.01 mm, standard deviation was 0.62 mm so that 90% of all necrotic areas occurred within approximately 2 mm (1.64 sigma, rounded). In the model, the occurrence of necrosis was determined by probSeedNecr, and it was more likely to be located at the tumor core, with the probability linearly increasing from 0 to 2 mm from the outer tumor margin.

Regarding the scale of the model, we measured cell density of tumor cells in 20 histological samples. We found that tumor cells and immune cells occupy an area of 222.8 μm² (median value). Although tumor cells are larger than lymphocytes, the uniform grid required that one grid cell should accommodate one cell. This assumption yielded a rectangular grid cells size of 14.9 μm. Thus, the length of 67 grid cells corresponds to 1 mm.

The events in our model (cell proliferation, migration, death) do necessarily occur with the same rate at each iteration. To account for these temporal differences, we introduced intermediary steps. We assumed that the median Ki67-positive fraction fk is equal to the probability of cell division pk in one iteration. Median tumor cell proliferation fraction (as measured in N = 20 tissue samples) was approximately 0.5 (Supplementary Fig. S1A) and a full cell cycle typically takes approximately 24 hours. Therefore, we set one full iteration of the model as 12 hours. To find out the probability of cell death in each iteration, we measured the fraction fc of tumor cells and immune cells positively stained for active caspase 3. The process of apoptosis induction to completion takes approximately 3 hours (38), just a quarter of the time step in our model. Therefore, the probability of cell death pc in each iteration was scaled appropriately. Regarding immune cell movement, we assumed that tumor-infiltrating lymphocytes (immune cell agents) migrate with an average speed of approximately 2 μm/min (2,880 μm/24 hours; ref. 32). These data were also qualitatively validated by our in vitro experiments (data not shown). In our model, this corresponds to 97 grid cells/12 hours (97 grid cells/iteration). Therefore, we introduced intermediary steps and allowed immune cells to move up to 97 times per iteration. Lastly, we scale the tumor cell killing events: It has been shown that tumor cell killing by cytotoxic T-lymphocytes is initiated within minutes (39) but takes approximately 6 hours to complete in vitro (40) and in vivo (31). Accordingly, we required that a killing event keeps a lymphocyte engaged for 6 hours and only after this period the lymphocyte can kill again (if it is not exhausted). It has been shown that tumor cell killing might require several hits by lymphocytes (40). This was not explicitly modeled in our system. Instead, we assumed that the killing probability parameter IMpkill already included multiple hits.

In our simulation runs, we let the tumors grow for a fixed number of iterations (nSteps) before performing changes of the parameters. After an additional number of steps addSteps, tumor cell number was compared with the baseline. We assessed the outcome analogously to the RECIST criteria (41) based on the number of tumor cells at T = nSteps+addSteps as compared with T = nSteps. Complete remission (CR) was equivalent to the eradication of all tumor cell, partial remission (PR) was a reduction in tumor cell number by at least 30%, progressive disease (PD) was an increase of tumor cell number by at least 20% and stable disease (SD) described all other outcomes.

All simulations were implemented in MATLAB (Mathworks) R2017a. Parts of the code were run in parallel with MATLAB's Parallel Processing Toolbox. All experiments were run on a standard workstation (Intel i7 Processor, 8 cores, 32 GB RAM, Microsoft Windows 10.1). Typically, computing speed was several hundred simulations per hour. We release all source codes for the agent-based model under an open-source license (http://dx.doi.org/10.5281/zenodo.853342).

To validate the predictions of the model, a clinical validation cohort of N = 261 colorectal adenocarcinoma (COAD) patients from the NIH (National Institutes of Health) The Cancer Genome Atlas (TCGA) collective was used (42). The data were downloaded via the TCGA Data Portal as described before (43). All TCGA samples that were included in the analysis are listed in Supplementary Table S4.

To assess the association between overall survival, stroma, lymphocytes, and their interaction, we performed a survival analysis of the TCGA collective. Cox proportional hazard models were fitted, including overall survival as a dependent variable, and lymphocytes (high/low), stroma (high/low), and the interaction of lymphocytes and stroma as fixed factors. Lymphocyte infiltration and stroma content of tumor tissue were part of the data tables available at the TCGA data portal. These variables had been manually measured by pathologists as part of the original TCGA data curation. Patients were split into low and high at the median. Due to its role as a potential confounder, TNM status was also included into the model as a fixed factor. Hazard ratios for effect estimates with corresponding 95% confidence intervals, and P values for hazard ratios and the interaction term were computed. P values smaller than 0.05 were regarded as statistically significant. The analysis was carried out using SAS v9.4 (SAS Institute). Furthermore, a Kaplan–Meier plot displaying the product-limit survival estimates alongside the number of subjects at risk for each of the four strata was created.

We propose a new agent-based model of tumor-immune cell interactions that is based on a minimal set of assumptions (Table 1) and parameters (Supplementary Table S1). In this model, we observed emergent behavior on different scales, particularly with regard to tumor tissue morphology. Generally, four types of immunological phenotypes can be distinguished histologically: hot tumors, cold tumors, immune excluded tumors (with an immune cell rim around the tumor) and tumors that have been (almost) completely eradicated by immune cells (44). Our agent-based model was able to reproduce all these spatial patterns as shown in Fig. 2A–D. We concluded that the model is in principle able to model all relevant types of immune surveillance in solid tumors.

To compare our model with histological spatial patterns in an objective way, CD3-stained colorectal cancer tissue samples were used. This analysis entailed all 20 tissue samples from the calibration cohort (Supplementary Table S2) and 37 additional samples from an independent validation cohort from a different institution (Supplementary Table S3). We analyzed spatial features of tumor/stroma distribution and found that the spatial layout of CD3-positive lymphocytes was always part of the spectrum of cold tumors (Fig. 2A), hot tumors (Fig. 2B), immune excluded tumors (Fig. 2C), or eradicated tumors (Fig. 2D). Specifically, resulting tumor nodules showed varying degrees of fibrosis and necrosis, mirroring spatial patterns in histological samples. This is shown in Supplementary Fig. S2A–S2C for an immune-excluded non-necrotic tumor, in Supplementary Fig. S2D– S2F for a largely necrotic tumor and in Supplementary Fig. S2G–S2I for a cold tumor with stromal core. We conclude that based on a calibration cohort and a morphological validation cohort, the model reproduced spatial architecture of tumor cells, stroma, and lymphocytes sufficiently.

To better elucidate the emergent dynamics, the behavior of tumors under different environmental conditions was investigated. In particular, we investigated immunological dynamics (immune surveillance) of a typical tumor with varying tumor–stroma ratios and varying numbers of tumor-reactive lymphocytes. To this end, we simulated the growth of 50 tumors for 60 days. At this point, tumors had reached a size of close to 5,000 cells. Then, the immune cell influx was strongly increased (immune boost) and the parameter for stroma induction was varied to generate tumors with different lymphocyte and stroma contents (Fig. 3A). According to these variations, four types of host response were investigated: Low and high stroma (Stro) generation (fibrosis seeding) and low and high immune cell (Lym) number.

We observed that the groups showed a drastically different behavior: Tumors in the “Stro low, Lym low” group showed an unhindered, exponential growth (Fig. 3B). In comparison, the growth in the “Stro high, Lym low” group was slower, but still steadily rising. As can be seen in Fig. 3C, both “Lym low” groups presented largely with progressive disease (PD) states at six months (180 days) after baseline. In contrast, a higher number of lymphocytes (Lym high) led to a phenotype with restrained tumor growth, as expected. Specifically, in the subgroup “Stro high, Lym high”, tumor size was constrained to under 10,000 cells (Fig. 3B and C). Surprisingly, the “Stro low, Lym high” group showed an altogether different behavior: after an initial increase in tumor mass (Fig. 3B), inflowing immune cells regained control and completely eradicated the tumor in almost all simulation runs (Fig. 3C; Supplementary Fig. S3).

In the model, the subgroup of tumors “Stro low, Lym high” by far the most favorable outcome of all simulated tumors. Tumor eradication reproducibly occurred only in the “Stro low, Lym high” group and not in the “Stro high, Lym high” group. This suggests that a high number of lymphocytes can only successfully constrain tumor growth if there is little stroma in the tumor.

We validated this prediction by analyzing a cohort of 261 patients from the TCGA database based on publicly available records (42, 43). For all patients, a manual histopathological quantification of stroma and lymphocytes was available as well as clinical follow-up data. Patients were stratified into high and low stroma content and lymphocyte number at the median. The only subgroup with a significant overall survival benefit as compared with the other groups was “Stro low, Lym high” in comparison with “Stro high, Lym high” group, as assessed by our Cox proportional hazards model (hazard ratio 0.309 for overall survival 0.322, P = 0.0219, confidence interval 0.122 to 0.849, Fig. 4A), reflected also by Kaplan–Meier curves for all groups (Fig. 4B). “Stro” and “Lym” alone were no significant predictors of overall survival, and neither were all other subgroups as shown in Fig. 4A. The P value for the interaction between “Lym” and “Stro” was P = 0.1192. Taking into account that the interaction test is commonly subject to a very small power (45), P = 0.1192 for an interaction test, even though not statistically significant, can be deemed as an indicator for a quite prominent interaction between “Lym” and “Stro.”

Having confirmed stroma as an important modulating factor in the in silico experimental setup, we simulated an immunotherapy together with stroma targeting therapy. A simple yet realistic way of simulating stroma-targeted therapy is to modify the permeability of stroma with regard to cell migration (46, 47). Therefore, we let tumors grow to a diameter of approximately 2 mm, which was typically reached after 120 days (Fig. 5A). Then, we simulated an “immune boost,” increasing the number of immune cells 2-fold, 4-fold, and 8-fold. Also, stroma was permeabilized with a factor of 4%, 8%, and 16%. In accordance with our previous results, we saw that nonpermeable stroma inhibited a tumor eradication (Fig. 5B).

In tumors with low stroma permeability and/or low lymphocyte numbers, we observed mostly tumor progression, characterized by four typical phenotypes: stroma acting as a physical barrier that protects the tumor (Supplementary Fig. S4A), rapid tumor outgrowth of immune control (Supplementary Fig. S4B), tumors breaking through a physical barrier (Supplementary Fig. S4C), and excessive immune cell exhaustion (Supplementary Fig. S4D).

Contrariwise, tumors with high stroma permeability and high lymphocyte numbers were successfully eradicated in 75% of all simulation runs (Fig. 5B). There was a striking duality in the effect of stromal permeabilization: In a lymphocyte-deprived environment, permeabilizing stroma had an adverse effect and led to increased tumor progression (bottom row in Fig. 5B). With highly permeable stroma, this adverse effect persisted also in 2-fold immune boosting and was only superseded by a strong immune boost of 4-fold.

From these in silico experiments with clinical validation, we conclude that stroma, arising through postinflammatory fibrosis, has a dual role in solid tumors: under usual conditions, it can mitigate tumor growth to a small degree. After a simulated immune boost, stroma provides a mechanism for immune escape. Only the permeabilization of stroma in combination with immunotherapy can lead to tumor regression in this model.

Within the concept of tumor immune surveillance, tumor cells and immune cells are engaged in an ongoing battle (26). From a conceptual perspective, a tumor disease develops if tumor cells temporarily win this fight by means of immune evasion. In line with this, immunotherapy aims at activating the T cells and thereby leading to tumor eradication in some patients. For colorectal cancer (microsatellite stable tumors), so far, no immunotherapy has shown efficacy. Driven by the clinical need to better understand the processes of immune evasion and activation, we developed a computer-based model that allowed us to investigate the dynamics governing tumor-immune cell interaction on a systems level (Fig. 1A). This model was generated utilizing quantitative histological and other wet lab data as well as data from previous studies (Supplementary Table S1).

The model was fitted to histological data in two ways: First, quantitative histological data from Ki67 and active caspase 3 immunostainings was used to estimate the proportion of actively dividing and dying cells as well as the overall cell density. In a second step, N = 37 human tumor tissue samples were used to show that the naturally occurring spatial patterns could be reproduced by our model, even on a large spatial scale (Supplementary Fig. S2A–S2I).

Using this validated in silico model, a consistent and highly relevant interaction between lymphocytes and stroma was discovered: In an immune-deprived environment, stroma restrains tumor growth and tumors with little stroma grows faster. In contrast, in an immune-cell–enriched environment, stroma inhibits tumor cell killing; thus, immune-cell rich tumors with low amounts of stroma are eradicated while immune-cell rich tumors with high amounts of stroma are not (Fig. 3). This leads to a bivariate risk-factor model for colorectal cancer that was validated in two independent validation steps: a morphological validation in n = 57 tissue samples (Fig. 2A–D, Supplementary Fig. S2A–S2I) and a clinical validation in a set of n = 261 patients (Fig. 4A and B).

As a next step, immune cell numbers and stromal permeability were gradually varied (Fig. 5B). Strikingly, a very consistent pattern was apparent: increasing the number of lymphocytes could only lead to tumor eradication if stromal permeability was also increased, in a dose-dependent manner (Fig. 6). Stromal permeability increase without lymphocyte-enrichment was detrimental and led to faster tumor progression (Fig. 6). This has important clinical consequences: increasing stromal permeability might greatly enhance the effectiveness of immunotherapy, but it can also be dangerous if the tumor microenvironment is lymphocyte deprived. Currently, several experimental therapies in the clinic aim to target the stroma. Indeed, most prominent are approaches with enzymes that aim to enhance the permeability of stroma (4, 5, 48). Our data show that these approaches can be very effective but can also have adverse effects that need to be considered before clinical development of new substances. Thus, for future clinical trials, the balance between stroma-targeting interventions and lymphocyte-targeting interventions should be investigated at the preclinical stage. In silico models such as our model could be a part of this preclinical pipeline.

The contribution of stroma reaction and inflammation in colorectal cancer has been investigated before, jointly for both factors (49) or separately (50, 51). However, to our knowledge, we report for the first time a bivariate interaction effect of both risk factors and provide an explanation for this behavior through in silico experiments. Our findings regarding fibrosis are in line with previous observations in biological studies and explain previous observations. Pancreatic cancer, for example, typically has pronounced tumor fibrosis. It has been shown that reducing fibrosis in this cancer benefits the tumor and reduces survival in a mouse model (52). This behavior also arises in our model. Furthermore, in line with our simulations, it has been shown that tumors can use the stroma to mitigate the immune response (53). Regarding immunotherapy, it has been suggested that immune response against a tumor is more efficient if the immune cells also destroy the stroma (48). Our model provides a simple explanation for this behavior.

Like all models, our model is not an exact copy of reality but simplifies some key aspects of real-world tumors. For example, killing of tumor cells by immune cells is a purely stochastic process in our model. Fine-tuning through T cell specificity, tumor immunogenicity and other co-stimulatory or inhibitory factors are not explicitly modeled but are summarized in the “killing probability”. Also, the model does not explicitly include myeloid cells or other regulating factors for lymphocyte-mediated cytotoxicity. In the future, myeloid cells—and the crucial regulation of lymphocyte activity via macrophages (6)—could be incorporated into more complicated models.

While these simplifications constrain the implications that can be drawn for certain molecular or signaling processes, it allows us to focus on the emergent features of the cellular interaction between key players in the tumor microenvironment. Most importantly, our model takes into account several quantitative histological observations and is thus partly calibrated with real-world tumors. Also, our model includes postinflammatory fibrosis, which gives rise to unexpected and informative emergent behavior. It therefore guides the possible translational developments, especially in selecting the ideal combination partners for synergistic clinical effects (Fig. 6). This provides therefore a roadmap for an iterative enhancement of immunotherapy for immunologically “cold” tumors.

No potential conflicts of interest were disclosed.

Conception and design: J.N. Kather, F. Klupp, N. Halama

Development of methodology: J.N. Kather, J. Poleszczuk, F. Klupp, N. Halama

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J.N. Kather, L. Tavernar, E. Herpel, F. Klupp, M. Schneider, A. Marx, N. Halama

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J.N. Kather, J. Krisam, P. Charoentong, N.A. Valous, C.-A. Weis, L. Tavernar, F. Leiss, A. Marx, N. Halama

Writing, review, and/or revision of the manuscript: J.N. Kather, J. Poleszczuk, M. Suarez-Carmona, J. Krisam, C.-A. Weis, F. Klupp, A. Ulrich, M. Schneider, A. Marx, D. Jäger, N. Halama

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): J.N. Kather, P. Charoentong, E. Herpel, D. Jäger, N. Halama

Study supervision: F. Klupp, D. Jäger, N. Halama

The authors would like to thank Anita Heinzelmann, Rosa Eurich and Jana Wolf (National Center for Tumor Diseases, Heidelberg, Germany), Katrin Wolk (University Medical Center Mannheim, Mannheim, Germany), and Nina Wilhelm (NCT Biobank, National Center for Tumor diseases, Heidelberg, Germany) for expert technical assistance. The authors are very grateful to Dr. Charles Neu (University Hospital Jena, Jena, Germany) for proofreading the article.

The results shown here are in part based upon data generated by the TCGA Research Network: http://cancergenome.nih.gov/.

J.N. Kather is supported by the “Heidelberg School of Oncology” (NCT-HSO) and by the “German Consortium for Translational Cancer Research” (DKTK) fellowship program. A. Marx is supported by a grant of the German Federal Ministry of Education and Research (BMBF) within the Framework of the Research Campus M²oBITE (grant 13GW0091E).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.