Immunotherapy has emerged as a promising study in cancer treatment, focusing on the dynamic interaction between cancer cells and the immune system in a tumor microenvironment. Within this context, hypoxia, a common phenomenon in tumors, has been shown to have a significant impact on tumor growth dynamics and response to therapy [⁹]. Macrophages, immune system cells that can phagocytize cancer cells, play a crucial role in the immune response against cancer.

• ⁹
  Acute vs. chronic vs. intermittent hypoxia in breast cancer: A review on its application in in vitro research

  Molecular Biology Reports, 2022
  9. Liu, Q., Palmgren, V. A. C., Danen, E. H., Le-Dévédec, S. E. (2022). Acute vs. chronic vs. intermittent hypoxia in breast cancer: A review on its application in in vitro research. Molecular Biology Reports, Vol. 49, No. 11, pp. 10961–10973. DOI: 10.1007/s11033-022-07802-6.

However, their function can be influenced by the tumor environment, including oxygen levels [¹³]. There are limitations when performing in vitro clinical trials on immunotherapy and hypoxia, so it is important to rely on computational simulation [⁴]. Computational simulation has become an invaluable tool for understanding the complexity of these interactions and their therapeutic implications in immunotherapy [³]. This technology allows researchers to model and analyze various aspects of the immune system’s behavior and its response to different treatments.

• ¹³
  Hypoxia signaling in cancer: From basics to clinical practice

  Pathology and Oncology Research, 2021
  13. Sebestyén, A., Kopper, L., Dankó, T., Tímár, J. (2021). Hypoxia signaling in cancer: From basics to clinical practice. Pathology and Oncology Research, Vol. 27. DOI: 10.3389/pore.2021.1609802.
• ⁴
  PhysiCell: An open source physics-based cell simulator for 3-D multicellular systems

  PLOS Computational Biology, 2018
  4. Ghaffarizadeh, A., Heiland, R., Friedman, S. H., Mumenthaler, S. M., Macklin, P. (2018). PhysiCell: An open source physics-based cell simulator for 3-D multicellular systems. PLOS Computational Biology, Vol. 14, No. 2, pp. e1005991. DOI: 10.1371/journal.pcbi.1005991.
• ³
  BioFVM: An efficient, parallelized diffusive transport solver for 3-D biological simulations

  Bioinformatics, 2015
  3. Ghaffarizadeh, A., Friedman, S. H., Macklin, P. (2015). BioFVM: An efficient, parallelized diffusive transport solver for 3-D biological simulations. Bioinformatics, Vol. 32, No. 8, pp. 1256–1258. DOI: 10.1093/bioinformatics/btv730.

The advantages of simulation include saving money and time and the absence of contamination of in vitro [⁷] cultures. Another significant benefit is the ability to run multiple scenarios and variations in a controlled environment. In this context, this study explores how alterations in oxygen (O2) levels in the tumor microenvironment affect tumor growth in immunotherapy.

• ⁷
  Physicell studio: A graphical tool to make agent-based modeling more accessible

  Gigabyte, 2024
  7. Heiland, R., Bergman, D., Lyons, B., Waldow, G., Cass, J., Lima-da-Rocha, H., Ruscone, M., Noël, V., Macklin, P. (2024). Physicell studio: A graphical tool to make agent-based modeling more accessible. Gigabyte, Vol. 2024, pp. 1–19. DOI: 10.46471/gigabyte.128.

Using a simulation environment based on PhysiCell [³], oxygen levels were systematically altered to observe their effect on cancer cell proliferation over time, and immune response was activated on day 13 of each simulation. The results provide crucial information on how hypoxia and immunotherapy interact in the tumor microenvironment and how these findings can inform the development of more effective therapeutic strategies for cancer treatment.

• ³
  BioFVM: An efficient, parallelized diffusive transport solver for 3-D biological simulations

  Bioinformatics, 2015
  3. Ghaffarizadeh, A., Friedman, S. H., Macklin, P. (2015). BioFVM: An efficient, parallelized diffusive transport solver for 3-D biological simulations. Bioinformatics, Vol. 32, No. 8, pp. 1256–1258. DOI: 10.1093/bioinformatics/btv730.

Some studies investigated how hypoxia may affect the immune response and the efficacy of immunotherapy in vitro. Table 1 summarizes various studies on the effects of hypoxia on different breast cancer cell lines. The studies use different breast cancer cell lines, including MDA-MB-231, MCF-7, and MMTVPyMT, allowing comparison of the effects of hypoxia in different breast cancer models. The use of hypoxia chambers predominates, although the AnaecroGen system is also mentioned, and they vary considerably between studies, both in the duration of intermittent hypoxia cycles and in the periods of acute or chronic hypoxia.

• ⁸
  Hypoxia-inducible factor 1 mediates intermittent hypoxia-induced migration of human breast cancer MDA-MB-231 cells

  Oncology Letters, 2017
  8. Liu, L., Liu, W., Wang, L., Zhu, T., Zhong, J., Xie, N. (2017). Hypoxia-inducible factor 1 mediates intermittent hypoxia-induced migration of human breast cancer MDA-MB-231 cells. Oncology Letters, Vol. 14, No. 6, pp. 7715–7722. DOI: 10.3892/ol.2017.7223.
• ⁶
  Hypoxia is a key driver of alternative splicing in human breast cancer cells

  Scientific Reports, 2017
  6. Han, J., Li, J., Ho, J. C., Chia, G. S., Kato, H., Jha, S., Yang, H., Poellinger, L., Lee, K. L. (2017). Hypoxia is a key driver of alternative splicing in human breast cancer cells. Scientific Reports, Vol. 7, No. 1. DOI: 10.1038/s41598-017-04333-0.
• ²
  Intermittent hypoxia induces a metastatic phenotype in breast cancer

  Oncogene, 2018
  2. Chen, A., Sceneay, J., Gödde, N., Kinwel, T., Ham, S., Thompson, E. W., Humbert, P. O., Möller, A. (2018). Intermittent hypoxia induces a metastatic phenotype in breast cancer. Oncogene, Vol. 37, No. 31, pp. 4214–4225. DOI: 10.1038/s41388-018-0259-3.
• ¹
  The effect of cycling hypoxia on MCF-7 cancer stem cells and the impact of their microenvironment on angiogenesis using human umbilical vein endothelial cells (HUVECs) as a model

  PeerJ, 2019
  1. Alhawarat, F. M., Hammad, H. M., Hijjawi, M. S., Sharab, A. S., Abuarqoub, D. A., Al-Shhab, M. A., Zihlif, M. A. (2019). The effect of cycling hypoxia on MCF-7 cancer stem cells and the impact of their microenvironment on angiogenesis using human umbilical vein endothelial cells (HUVECs) as a model. PeerJ, Vol. 7, pp. e5990. DOI: 10.7717/peerj.5990.

All studies use 1 % O2 to induce hypoxia, but exposure durations and frequencies vary. The table provides a clear and comparative view of how different hypoxia conditions can influence in vitro studies of breast cancer, which is crucial for understanding the underlying mechanisms and developing potential treatments.

Performing these in vitro assays means that if the appropriate protocol is not followed, the cell cultures are sometimes no longer optimal for use in research, so it is important to use other alternatives, such as computational simulations. Table 2 summarizes studies that developed models for immunotherapy in cancer treatment.

• ¹²
  Modeling small cell lung cancer (SCLC) biology through deterministic and stochastic mathematical models

  Oncotarget, 2018
  12. Salgia, R., Mambetsariev, I., Hewelt, B., Achuthan, S., Li, H., Poroyko, V., Wang, Y., Sattler, M. (2018). Modeling small cell lung cancer (SCLC) biology through deterministic and stochastic mathematical models. Oncotarget, Vol. 9, No. 40, pp. 26226–26242. DOI: 10.18632/oncotarget.25360.
• ¹⁴
  Virtual clinical trials of anti-PD-1 and anti-CTLA-4 immunotherapy in advanced hepatocellular carcinoma using a quantitative systems pharmacology model

  Journal for ImmunoTherapy of Cancer, 2022
  14. Sové, R. J., Verma, B. K., Wang, H., Ho, W. J., Yarchoan, M., Popel, A. S. (2022). Virtual clinical trials of anti-PD-1 and anti-CTLA-4 immunotherapy in advanced hepatocellular carcinoma using a quantitative systems pharmacology model. Journal for ImmunoTherapy of Cancer, Vol. 10, No. 11, pp. e005414. DOI: 10.1136/jitc-2022-005414.
• ¹¹
  Modeling cancer immunoediting in tumor microenvironment with system characterization through the ising-model Hamiltonian

  Bioinformatics, 2022
  11. Rojas-Domínguez, A., Arroyo-Duarte, R., Rincón-Vieyra, F., Alvarado-Mentado, M. (2022). Modeling cancer immunoediting in tumor microenvironment with system characterization through the ising-model Hamiltonian. BMC Bioinformatics, Vol. 23, No. 1. DOI: 10.1186/s12859-022-04731-w.
• ⁵
  Modeling tumor growth using fractal calculus: Insights into tumor dynamics

  BioSystems, 2024
  5. Golmankhaneh, A. K., Tunç, S., Schlichtinger, A. M., Asanza, D. M., Golmankhaneh, A. K. (2024). Modeling tumor growth using fractal calculus: Insights into tumor dynamics. BioSystems, Vol. 235, pp. 105071. DOI: 10.1016/j.biosystems.2023.105071.
• ¹⁰
  Simulation tumor growth in heterogeneous medium based on diffusion equation

  International Journal of Modern Physics C, 2023
  10. Polyakov, M. V., Ten, V. V. (2023). Simulation tumor growth in heterogeneous medium based on diffusion equation. International Journal of Modern Physics C, Vol. 35, No. 1. DOI: 10.1142/s0129183124500104.

These studies use a variety of simulation approaches and tools. For example, Salgia (2018) and Golmankhaneh (2023) employ deterministic and stochastic models in MatLab to simulate the behavior of cancer cells in 2D environments.

Wei (2022) uses ordinary differential equations (ODEs) to model different types of immune and cancer cells also in MatLab. Rojas-Domínguez (2022) adopts an approach based on the Ising Hamiltonian model using NetLogo to study the interaction between cancer cells and T lymphocytes.

Polyakov (2023) applies diffusion equations in 2D and 3D environments using C++ and Python to simulate the growth of cancerous tumors. Although these models offer valuable insights into the interaction between cancer cells and immunotherapeutic treatments, they have certain limitations compared to PhysiCell.

Many of these models are restricted to 2D environments, which can limit the accuracy of the simulations by not considering the 3D complexity of tumors. While powerful, the simulation environments, such as MatLab and NetLogo, may not be as efficient or scalable as PhysiCell, specifically optimized for large-scale 3D multicellular simulation.

PhysiCell also provides more robust integration with experimental biological data and is highly adaptable to incorporate new biological features and mechanisms, which is crucial for accurate and up-to-date immunotherapy models.

In addition, oxygen levels can be altered in PhysiCell, which is essential for studying the tumor microenvironment and its impact on the effectiveness of treatments.

Cancer cells are cells that have undergone genetic transformations that alter their growth and division cycle, allowing them to proliferate in an uncontrolled manner. In that context, the code provided describes the behavior of hanging tumor spheroids in the simulation model, which are very common in breast cancer.

The life cycle details include a continuous growth cycle and two cell death pathways, apoptosis, and necrosis, with specific rates and durations for each process. The cell has a defined total volume with specific proportions of biomass and fluid, and its mechanical dynamics include cell-cell adhesion and repulsion.

In addition, the cell’s motility is disabled, reflecting its inability to move autonomously. The rates of secretion and uptake of substances such as oxygen, essential for their survival and growth, are also specified. In Fig. 1, the way cancer cells interact with their environment and respond to the oncogene to proliferate uncontrollably is presented.

A macrophage is a type of white blood cell essential in the immune system, responsible for detecting, engulfing, and destroying pathogens and dead or damaged cells by a process known as phagocytosis. This cell plays a crucial role in the immune response alongside B cells (see Fig. 2). B cells are a type of lymphocyte originating in the bone marrow and playing a crucial role in the humoral immune response.

They produce and secrete antibodies, proteins that bind to antigens to neutralize or mark them for destruction by macrophages. Macrophages in simulation have a phenotype that initially disables their motility capabilities, with the speed of movement set at 1 micron per minute and a persistence time of 10 minutes. Migration is partially biased (0.5).

PhysiCell is designed to investigate the dynamics and interactions of thousands or millions of cells in three-dimensional microenvironments with environment-dependent phenotypes. It employs a physical, lattice-free approach to minimize artifacts associated with grids.

After initializing the microenvironment using BioFVM [⁵] and cells, as well as the current simulation time t=0, PhysiCell internally tracks tmech (the next time at which cell mechanics functions are executed), tcell (the next time at which cellular processes are executed) and tsave (the next simulation data output time), with an output frequency tdelta, tsave. We initially set:

• ⁵
  Modeling tumor growth using fractal calculus: Insights into tumor dynamics

  BioSystems, 2024
  5. Golmankhaneh, A. K., Tunç, S., Schlichtinger, A. M., Asanza, D. M., Golmankhaneh, A. K. (2024). Modeling tumor growth using fractal calculus: Insights into tumor dynamics. BioSystems, Vol. 235, pp. 105071. DOI: 10.1016/j.biosystems.2023.105071.

The software repeats the main program loop until the maximum simulation time is reached. To learn more about the loop, see [³].

• ³
  BioFVM: An efficient, parallelized diffusive transport solver for 3-D biological simulations

  Bioinformatics, 2015
  3. Ghaffarizadeh, A., Friedman, S. H., Macklin, P. (2015). BioFVM: An efficient, parallelized diffusive transport solver for 3-D biological simulations. Bioinformatics, Vol. 32, No. 8, pp. 1256–1258. DOI: 10.1093/bioinformatics/btv730.

The BioFVM environment for simulation of the chemical microenvironment uses a vector of reaction-diffusion partial differential equations (PDE). The modeling of the biochemical microenvironment (within a Ω computational domain and its ∂Ω boundary, discretized as a Cartesian grid to optimize computational efficiency) as a vector of reaction-diffusion PDEs applied to a vector of ρ chemical substrates as follows:

is the diffusion

decay

mass source

mass uptake

sources and uptake by cells

With null flow conditions in ∂Ω. In this context, δ(x) represents the Dirac delta function, xk indicates the position of cell number k, Wk its volume, sk the vector detailing its secretion rates, Uk the vector listing its uptake rates, and P∗ the vector of saturation densities (density levels at which cells cease secretion).

In addition, D and λ correspond to the vectors of diffusion coefficients and degradation rates, respectively. S is the supply rate, and U is the uptake function. All products between ab vectors are performed element by element (Hadamard product).

Numerically, we solve the solution in time t+Δt by a first-order operator decomposition: first, the bulk source/sink terms are solved over the entire domain (and the stored solution is overwritten), then the cell-centered source/sink terms are solved (and the solution is overwritten), and finally the reaction-diffusion terms (again, the stored solutions are overwritten) [⁵].

• ⁵
  Modeling tumor growth using fractal calculus: Insights into tumor dynamics

  BioSystems, 2024
  5. Golmankhaneh, A. K., Tunç, S., Schlichtinger, A. M., Asanza, D. M., Golmankhaneh, A. K. (2024). Modeling tumor growth using fractal calculus: Insights into tumor dynamics. BioSystems, Vol. 235, pp. 105071. DOI: 10.1016/j.biosystems.2023.105071.

We mention the parameters for the microenvironment, which are presented in the Table 3. These parameters are associated with the development of cancer cells and can affect their tumor proliferation, specifically in the oxygenation cycle.

To facilitate the use of the code, we have provided a link on GitHub where models y1, y2, and y3 are included, along with the specific parameter modifications for each of them¹.

In this work, we define the microenvironment for developing cancer cells with the conditions specified in the section “Microenvironment”. The domain for the visualization of the simulations is in the second and third dimensions, and the configuration for simulation time is in minutes; for example, to simulate 25 days, it is necessary to multiply the number of days by 24 (the number of hours in a day), and then by 60 (minutes in an hour) that is 36 000 minutes.

We initially define a set of cancer cells that will proliferate according to the parameters mentioned in section 3.3. Also, macrophages will activate and start phagocytosis from day 13. Simulation y1 was performed under average conditions without any oxygen alteration.

For simulation y2, we applied the hypoxia effect, which is 1 % of oxygen (O2) for 25 days [⁸]. Finally, in y3, we increased the oxygen level. Fig. 3 is a general representation of the methodology, in which there are cancer cells, oncogenes for cancer cells, macrophages, and B cells. Fig. 3 is the general representation of the methodology in which there are, cancer cells, oncogene for cancer cells, macrophages, and B cells.

• ⁸
  Hypoxia-inducible factor 1 mediates intermittent hypoxia-induced migration of human breast cancer MDA-MB-231 cells

  Oncology Letters, 2017
  8. Liu, L., Liu, W., Wang, L., Zhu, T., Zhong, J., Xie, N. (2017). Hypoxia-inducible factor 1 mediates intermittent hypoxia-induced migration of human breast cancer MDA-MB-231 cells. Oncology Letters, Vol. 14, No. 6, pp. 7715–7722. DOI: 10.3892/ol.2017.7223.

The immune system activation in y1, y2, and y3 is on day 13. A specific number of macrophages perform phagocytosis starting at day 13 and up to day 25. For cell growth simulations, specific parameters that were modified for y1, y2, and y3 are defined (see Table 4).

Applying the methodology mentioned in section 4, for y1, y2, and y3 we used for the model hanging tumor spheroids, a macrophage activation on day 13, and a 25-day simulation in PhysiCell. The first simulation was performed under average normal conditions for y1. For y2, we increased the 2 % of O2. On the other hand, for y3, we decreased the 1 % of O2 according to [²].

• ²
  Intermittent hypoxia induces a metastatic phenotype in breast cancer

  Oncogene, 2018
  2. Chen, A., Sceneay, J., Gödde, N., Kinwel, T., Ham, S., Thompson, E. W., Humbert, P. O., Möller, A. (2018). Intermittent hypoxia induces a metastatic phenotype in breast cancer. Oncogene, Vol. 37, No. 31, pp. 4214–4225. DOI: 10.1038/s41388-018-0259-3.

Cancer cells represent the blue color; yellow represents an oncogene, black represents the nucleus of each cancer cell, and green represents macrophages. In Fig. 4, the result of y1 can be seen, which, in the biological context [³], is a simulation with average parameters.

• ³
  BioFVM: An efficient, parallelized diffusive transport solver for 3-D biological simulations

  Bioinformatics, 2015
  3. Ghaffarizadeh, A., Friedman, S. H., Macklin, P. (2015). BioFVM: An efficient, parallelized diffusive transport solver for 3-D biological simulations. Bioinformatics, Vol. 32, No. 8, pp. 1256–1258. DOI: 10.1093/bioinformatics/btv730.

For y2, to which the oncogene saturation number and threshold were modified, as well as increased oxygenation, differences are observed concerning y1, at least in the increase of cells (see Fig. 5). Finally, hypoxia was applied to y3, and the simulation was obtained as shown in Fig. 6. In Fig. 7 we show the behavior of the simulations and we can observe that concerning time in the three cases, there is an exponential trend of cancer cell growth.

The pairwise correlation coefficient was calculated for each combination obtaining: the correlation coefficient between y1 and y2=0.926, y1 and y3=0.921, finally y2 and y3=0.994. We performed one-factor ANOVA for y1, y2, and y3 and obtained that the test F=38.87 and the test value P= 3.5 × 10−12. So we performed pairwise ANOVA.

Comparison between the diameter of ductal carcinoma in situ (DCIS) and the simulated models reveals that DCIS, with a diameter of 1.511 ± 0.10 mm, differs from the models in several aspects: Model y1 shows a diameter of 1.422 ± 0.10 mm, which is slightly smaller than that of DCIS; model y2 presents a diameter of 1.6 ± 0.10 mm, which is larger, indicating an overestimation of growth; and model y3 has the smallest diameter of 1.288 ± 0.10 mm.

These discrepancies underscore the need to adjust model y3. It should be noted that this model includes hypoxia generation, suggesting that the tumor in this model needs more time to reach a diameter comparable to that of DCIS, as illustrated in Fig. 8(a).

This research stands out for its interdisciplinary approach, which combines elements of biology, computational modeling, and immunotherapy, a comprehensive understanding of the interaction between the immune system and the tumor microenvironment.

This approach not only offers innovative perspectives on cancer treatment but also underscores the importance of collaboration between different scientific fields to advance the development of more effective therapies.

The lack of significance in the combination of y1 and y3 could reflect a lower relevance of the interaction between these parameters in the simulated context, or it could suggest the need to explore other factors that may be influencing tumor dynamics.

Overall, these results provide valuable information for tuning models and designing more effective treatments, highlighting the need for a more detailed evaluation of how different combinations of parameters affect cancer progression. In biopsy images, cancer cells can be seen to overlap each other, as shown in Fig. 8 making it difficult to fully visualize all cells, both cancerous and noncancerous.

In our study, we performed 2D simulations, which limits the accuracy in visualizing certain patterns, such as the ellipsoid shape. Therefore, we propose to perform 3D simulations in future work to improve the accuracy of the representation of these patterns.

Tumor growth patterns that experience hypoxia can cause necrosis, and in some types of cancer, such as ductal carcinoma in situ (DCIS), this can lead to the formation of microcalcifications. These microcalcifications in the mammary duct are detectable by mammography and may allow early detection of cancer.

A future approach is to develop models that simulate the formation of microcalcifications from tumor growth, to improve early detection of DCIS on mammography.

Other future work is planned to extend the duration of the simulation beyond 25 days to obtain a complete picture of long-term tumor dynamics and to evaluate how variable oxygen levels sustainably influence cancer progression and treatment effectiveness.

The y2 model had the highest number of cancer cells during the 25 days, y1, and y3 had close behaviors until day 13 of immune system activation. Research has shown that sensitivity to altered oxygen levels in the microenvironment, as well as oncogene expression, have a significant impact on cancer cell proliferation within an immunotherapy model.

The high correlation between parameters y1, y2, and y3 suggests that changes in these oxygen parameters produce predictable and consistent effects on tumor growth, indicating the robustness of the models used.

Furthermore, the p-value significantly less than 0.05 supports the existence of significant differences between the means of models y1, y2, and y3, reaffirming the influence of the studied factors on tumor growth.

These findings underscore the importance of understanding how oxygen levels and oncoprotein expression affect cancer behavior. This knowledge may guide future research and contribute to developing more effective therapeutic strategies, such as manipulating oxygen levels to control tumor growth.