Mathematical modeling of hypoxia and adenosine to explore tumor escape mechanisms in DC-based immunotherapy

TL;DR

This paper uses mathematical models to study how hypoxia and adenosine help tumors escape immunotherapy, aiming to improve treatment strategies.

Contribution

The first use of mathematical models to explore hypoxia and adenosine's roles in tumor escape during DC-based immunotherapy.

Findings

Adenosine significantly impacts immunotherapy efficacy, with its suppression potentially enhancing treatment outcomes.

Mathematical models reveal how hypoxia and adenosine contribute to tumor escape mechanisms.

Optimized vaccination protocols can minimize tumor growth based on model predictions.

Abstract

Identifying and controlling tumor escape mechanisms is crucial for improving cancer treatment effectiveness. Experimental studies reveal tumor hypoxia and adenosine as significant contributors to such mechanisms. Hypoxia exacerbates adenosine levels in the tumor microenvironment. Combining inhibition of these factors with dendritic cell (DC)-based immunotherapy promises improved clinical outcomes. However, challenges include understanding dynamics, optimal vaccine dosages, and timing. Mathematical models, including agent-based, diffusion, and ordinary differential equations, address these challenges. Here, we employ these models for the first time to elucidate how hypoxia and adenosine facilitate tumor escape in DC-based immunotherapy. After parameter estimation using experimental data, we optimize vaccination protocols to minimize tumor growth. Sensitivity analysis highlights adenosine’s significant impact on immunotherapy efficacy. Its suppressive role impedes treatment success, but inhibiting adenosine could enhance therapy, as suggested by the model. Our findings shed light on hypoxia and adenosine-mediated tumor escape mechanisms, informing future treatment strategies. Additionally, identifiability analysis confirms accurate parameter determination using experimental data.

Figures

The time course of tumor volume as the parameters *µ* and *P*_*Treg*_ change. (a) Tumor volume as a result of the change in *µ*. (b) Tumor volume as a result of the change in *P*_*Treg*_.

The suggested model’s overall structure. In this model, tumor, DC, T cell, Treg, effector cell, and factors such as hypoxia and adenosine are considered. The interaction of the factors with each other is shown by the arrows and their corresponding parameter is specified.

Interactions and communication between tumor and immune system cells in the discrete grid over time period Δ*t*. The tumor cell activates its neighboring inactive DC (A1–B1). The T cell next to the active DC becomes a Treg (A3–B3) or effector cell (A5–B5) after activation. The effector cell kills the tumor cell adjacent to it (D1–D2). The Treg cell inactivates its neighboring active DC (D5–E4).

Tumor volume was predicted by the model and measured in the empirical data^[3](#CR3),[4](#CR4)^ in different groups. (a) Untreated. (b) Treatment with adenosine inhibitor. (c) Treatment with DC vaccine. (d) Treatment with adenosine inhibitor and DC vaccine. (e) Treatment with hypoxia inhibitor. (f) Treatment with hypoxia inhibitor and DC vaccine. (g) Treatment with adenosine and hypoxia inhibitor. (h) Treatment with adenosine and hypoxia inhibitor and DC vaccine.

The number of effector cells and Treg cells in different groups. (a) Untreated. (b) Treatment with adenosine inhibitor. (c) Treatment with hypoxia inhibitor. (d) Treatment with adenosine inhibitor and hypoxia inhibitor.

Comparison of the oxygen concentration and tumor cells in the untreated and siRNA + Px-478 + DC groups on day 12. (a) Oxygen concentration in the untreated group on day 12. (b) Oxygen concentration in the siRNA + Px-478 + DC group on day 12. (c) Tumor cells in the untreated group on day 12. (d) Tumor cells in the siRNA + Px-478 + DC group on day 12.

Investigating the time course of effector cell cytotoxicity rate and adenosine concentration in four groups. (a) Effector cell cytotoxicity rate, (b) the adenosine concentration.

The average number of effector cells around the tumor as a result of changing the dose of the DC vaccine.

Sensitivity value of the model parameters.

The time course of tumor volume, effector cell cytotoxicity rate, adenosine concentration, and number of Treg as a result of changes in *β* (a) Tumor volume. (b) Effector cell cytotoxicity rate. (c) Level of adenosine. (d) Number of Treg*.*

Full text

Introduction

Cancer is one of the most frequent causes of mortality globally^1,2^. Hypoxia and adenosine are two of the factors that compromise the immune system and contribute to the spread of cancer^3,4^.

Uncontrolled cell growth in the tumor microenvironment (TME) frequently exceeds the capacity of the preexisting blood capillaries to supply oxygen demand, resulting in hypoxic conditions^5^. Hypoxia has been demonstrated to protect tumors by promoting angiogenesis, reducing antitumor immune cells, and activating immunosuppressive cells like regulatory T cells (Tregs)^6^. Inhibiting the transcription factors known as hypoxia-inducing factors (HIFs), which react to reductions in the amount of available oxygen in the cellular environment, can be a useful strategy for limiting tumor growth^4,7,8^.

The ectoenzymes such as CD73 generate extracellular adenosine from ATP or ADP, which has a significant impact on the activation of immune system suppressor cells like Tregs, as well as tumor angiogenesis and metastasis^9^. So, CD73 is an attractive therapeutic target for the treatment of cancer, particularly when used in conjunction with conventional therapy (chemotherapy, radiotherapy)^3,10–12^. On the other hand, the relationship between adenosine and hypoxia has also been investigated in experimental studies and it has been shown that CD73 increases with decreasing oxygen levels^13–16^.

In dendritic cell (DC)-based immunotherapy, DCs are activated ex vivo and reinjected into the host to perform anti-cancer effector actions. DCs are antigen-presenting cells that process antigen material and present it on the cell surface to the T lymphocytes. After activation, T lymphocytes differentiate into different subtypes and migrate to destroy tumor cells or cells that have been infected by a virus^17–20^.

Mathematical models can be used to explain how cancer spreads, the immune system and the tumor interact, immune cells work, and to analyze the success of novel treatments and the combination of various therapies^21^. Examples of mathematical tools used to model living organisms are ordinary differential equations (ODEs) and agent-based models (ABMs)^22^. In ODE models, spatial and temporal dynamics cannot be simultaneously captured despite their simplicity. Additionally, these models do not consider stochasticity, although there are stochastic versions known as stochastic differential equations available. ABMs, as opposed to the conventional ODE modeling technique, provide more flexibility for individual action and are able to explain some of the behaviors observed in the tumor-immune system interaction, such as cell–cell interactions involved in the immune response, increased immune system access to tumor cells with tumor growth, and the impact of tumor antigen on the movement of immune cells within the TME^23^.

Many models in the form of ODEs and ABMs have been presented for tumor proliferation and describing the interaction between the tumor and the immune system along with various treatment methods. Carlos et al.^24^ modeled the interactions between immune system and tumor using ODE models, and used the model to test different DC vaccine injection protocols to find the amounts of DCs and injection times at which tumor treatment is significantly improved. Some studies such as^23,25^ show that DCs have nonlinear and random behavior and ABM models are more suitable than ODEs to model these behaviors. Macfarlane et al.^25^ proposed a discrete mathematical model of the interactions between immune system cells and monolayer tumor development. They investigated how cytotoxic T lymphocytes (CTLs) and DCs behave in the TME. The results indicated, the level of efficiency of T cells in destroying tumor cells is important and the immune response to cancer can be improved to some extent by raising the number of cytotoxic T cells, but an increase in the number of DCs may have the opposite effect. Pourhassanzadeh et al.^26^ introduced a novel two-dimensional random ABM to explore avascular tumor growth. In this model, the tumor is considered multilayer, and the behavior of the immune cells in the TME is investigated. Their results showed that immune cells are crucial in the prevention and treatment of cancer malignancy. Later, the authors offered a two-dimensional model of tumor growth with three layers^22^. Additionally, the concentration of oxygen, acid, and glucose in the TME was examined. The results showed that in the areas where there are tumor cells, oxygen is low and the TME becomes more acidic and leads to more aggressive tumor growth. However in this model, the immune system is not considered. Arabameri and Pourgholaminejad^23^ introduced a hybrid ABM-ODE model for studying adenosine’s role in immunotherapy. Their findings highlighted how crucial adenosine inhibition is to improve the immune system. However, in this model, the role of hypoxia and its impact on adenosine and tumor growth has not been incorporated.

To the best of the authors’ knowledge, no model has yet been proposed to quantify the relationship between hypoxia and adenosine, as well as the role of both in cancer immunotherapy. Therefore, the purpose of this study is to investigate the effect of oxygen deficiency on adenosine production and also the role of these two factors (hypoxia and adenosine) in providing tumor escape mechanisms in DC-based immunotherapy.

For the following reasons, we use a spatio-temporal model in this study: (1) the oxygen production or consumption in the TME is a function of location, that is, for example, a tumor cell consumes oxygen at a high rate, and immune system cells or normal body cells have a lower consumption rate, (2) an essential aspect of immune cell activity is how easily the immune system can access tumor cells, which depends on the geometrical pattern of the tumor’s outer layer. Therefore, in “Methods” section the proposed model is explained in detail. Then, in “Numerical simulations” section, numerical simulation of the model, finding the optimal vaccination pattern, and sensitivity identifiability analyses have been done. The findings are summarized and analyzed in “Discussion” section, along with a roadmap for future model improvements.

Methods

The overall structure of the model is depicted in Fig. 1. This diagram illustrates the interactions and connections between tumor, immune system cells, hypoxia, adenosine, and vaccines. While the model presented in this study has been calibrated using data specific to immunotherapy of breast cancer in mouse models, the interactions depicted in Fig. 1 are generally applicable to other types of solid tumors as well. We have made the assumption that there is a single population of inactive T cells, consisting of both naive Foxp3-expressing Tregs and naive CD8+ T cells. We have also assumed that this population can differentiate into both Tregs and effector cells. However, it is important to note that this is not the case in the immune system. This is done for the simplicity of the model and to reduce the number of variables. Hypoxia and adenosine increase the rate of tumor growth and decrease the cytotoxicity rate and the number of effector cells. Also, hypoxia has a positive effect on the production of adenosine. The vaccinations shown in this figure inhibit hypoxia and adenosine and increase the number of active DCs in TME.Figure 1. The suggested model’s overall structure. In this model, tumor, DC, T cell, Treg, effector cell, and factors such as hypoxia and adenosine are considered. The interaction of the factors with each other is shown by the arrows and their corresponding parameter is specified.

Modeling of tumor growth

We consider a two-dimensional discrete grid with size N × N. Time variable t is discretized with time step Δt and the space variables x and y are discretized with space step Δx = Δy = χ. We employ a Moore neighborhood structure, in which every position in the lattice is surrounded by eight neighbors. Initially, we place a tumor cell in the center of this grid and it proliferates at rate Ptumor. During division, the mother cell remains in its place and the daughter cell occupies one of the empty neighboring positions. If the tumor cell does not have an empty neighbor, it is not possible for it to divide^22,25^. As a result, only the cells comprising the outer layer of the tumor have the ability to proliferate, mirroring the conditions observed in the actual tumor environment. This reduction in effective growth rate has also been incorporated into previous studies^27^. It should be noted that we have omitted two phenomena related to the tumor for the sake of simplicity of our model: (1) the natural death of tumor cells, as the tumor experiences rapid growth in its early stages, making natural cell death negligible compared to its proliferation. (2) Metastasis, which typically does not manifest in the initial phases of tumor development. These assumptions have also been addressed in prior studies^23,25,28^.

Modeling of immune cells

Immune cell mobility can now be more easily measured thanks to improvements in imaging techniques^25^. Before becoming activated, immune system cells are free to roam around in larger areas. These cells move across a limited area after activation. It has been demonstrated in Ref.^25^ that for a more precise representation of the movement of inactive immune cells, employing the Levy walk model is preferable, whereas the movement of active immune cells is more suitably characterized by Brownian motion. Therefore, in the vein of Ref.^25^, we consider Levy motion and Brownian motion for inactive and active immune system cells, respectively. The moving step in Levy walk is randomly generated from the cumulative distribution function:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F\left(z\right)={1-z}^{-l}, z=1, 2,\dots ,$$\end{document}

where z is the number of movement steps in one direction and l is the walk exponent. After z is determined, a direction is randomly selected from empty neighbors. The cell will move in that time step and z − 1 time steps later in the selected direction. At each time step, the cell will not move if its selected path is occupied or out of the grid. Immune cells switch to Brownian motion when the tumor antigen activates them^25^. In this motion, at each time step there is an equal chance that the activated cells will move to one of the eight nearby grid sites. The cell can be moved when there is an empty space in its neighborhood and it is not outside the grid. The difference between these two motion types is the movement step. The cells in Levy walk take more steps in a certain direction^25^.

Four different types of immune cells are taken into account in our model: (I) DCs, (II) inactive T lymphocytes, (III) Treg cells, and (IV) effector cells. DCs and inactive T lymphocytes, whose numbers we denote by ND and NT, respectively, are initially inactive and randomly distributed throughout the grid, moving in a Levy walk fashion. Please note that we have not considered the proliferation and death of immune cells in this model, nor have we included the healthy cell, similar to Ref.^25^. Moreover, cells are confined within the network and prevented from exiting. Upon nearing the network boundaries, they are repositioned in a direction opposite to the edge, employing either the Levy walk or the Brownian motion based on their movement prior to reaching the boundary.

Figure 2 displays a schematic example of a cross-sectional portion of the ABM before and after a time step Δt. Inactive DCs may become active in interaction with cancer cells. If an inactive DC is next to a cancer cell (Fig. 2A1–B1), it will be activated at rate Dact and then it moves in the grid with the Brownian motion. An activated DC can activate the inactive T lymphocytes located in its neighborhood at rate Cact (Fig. 2A3–B3,A5–B5). It should be noted that immunogenicity or tolerogenicity of DCs depends on the environmental conditions during their maturation process^29^. During activation, T-lymphocytes convert to Treg at rate PTreg (Fig. 2A3–B3) or effector cells at rate 1−PTreg (Fig. 2A5–B5). An activated DC can activate several T lymphocytes. The parameter PTreg can be considered as the tolerogenicity of activated DCs. Effector cells destroy their neighboring tumor cells at rate μ (Fig. 2D1,D2). At every time step, each effector cell has the ability to eliminate a tumor cell situated in its neighborhood. If there is more than one tumor cell in the neighborhood of the effective cell, one of them is randomly selected and removed with probability µ. Subsequently, when a tumor cell is eliminated, the corresponding location is designated as an empty space. It is shown that effector cells can make multiple contacts and kill serially in a time-dependent manner^30^. Therefore, we assume that the effector cell can continue to move and kill cancer cells again. Active DCs located in the neighborhood of a Treg are deactivated at rate Dinact (Fig. 2D5–E4).Figure 2. Interactions and communication between tumor and immune system cells in the discrete grid over time period Δt. The tumor cell activates its neighboring inactive DC (A1–B1). The T cell next to the active DC becomes a Treg (A3–B3) or effector cell (A5–B5) after activation. The effector cell kills the tumor cell adjacent to it (D1–D2). The Treg cell inactivates its neighboring active DC (D5–E4).

Hypoxia and adenosine modeling

To model the relationship between oxygen and hypoxia, it is apparent that a decrease in oxygen concentration leads to a reduction in enzymatic reactions. Based on the relationships depicted in Fig. 5 of Ref.^31^, one can deduce a straightforward relationship: the hypoxia level is inversely proportional to the oxygen concentration. Therefore, in all grid sites, we considered hypoxia to be the inverse of oxygen concentration. We use a discretized diffusion equation to update the oxygen concentration at all sites of the grid in each time step. To reduce the required computation, we used the method presented in Ref.^22^. In each time step, the grid is divided into 3 × 3 non-overlapping blocks and the concentration of oxygen at the center of each block is calculated as below:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${U}_{m,n}^{t+1}={\alpha }{\times A}_{m,n}^{t}+\left(1-\alpha \right)\times {U}_{m,n}^{t}-{f}\left(n,m,t\right)+ kpx,$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${A}_{m,n}^{t}=(\sum_{j=-1}^{j=1} \sum_{k=-1}^{k=1}{U}_{m-j, n-k}^{t})/9,$$\end{document}

and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${U}_{mn}^{t}$$\end{document} represents the concentration of oxygen and* m*, n, and t denote the column number of the spatial position, row number of the spatial position, and the time step, respectively. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${f}(n,m,t)$$\end{document} is a consumer or producer function and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f\left(m,n,t\right)=\overline{f}I\left(m,n,t\right)$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{f}$$\end{document} is a constant and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I\left(m,n,t\right)$$\end{document} is an indicator function with value 1 in tumor sites and 0 elsewhere. The parameter α ϵ (0,1) represents the diffusion coefficient and parameter k is related to the effect of the vaccine on oxygen. Hypoxia inhibitor vaccine is modeled in (5). After calculating the concentration of oxygen at the center of each 3 × 3 block, the calculated value is replaced and redistributed in the whole block.

In the next step, adenosine and the effect of hypoxia on adenosine are modeled as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{dAd}{dt}= \beta +\frac{a}{U+c}-rAd-\gamma sAd.$$\end{document}

The level of adenosine production in cells can vary depending on the specific cell type and physiological conditions. Generally, there is a base level, denoted as β, for adenosine production rate^32^. However, under certain circumstances such as low oxygen levels, the production of adenosine can exceed this base level. This increase in adenosine production serves as a protective response to the decreased oxygen availability^33^. The parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a$$\end{document} adjusts the strength of adenosine production by hypoxia. As mentioned before, we have considered that the hypoxia level is inversely proportional to the oxygen concentration. We have added a near-zero value (the parameter c) to the denominator to avoid division by zero. The parameter r is the elimination rate of adenosine, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} models the adenosine inhibition by vaccine. The parameters β, a, and r are fine-tuned to closely match the experimental data on adenosine dynamics, as illustrated in Fig. 5 of Ref.^14^. The parameter γ is estimated using experimental data to guarantee that the levels of adenosine show proportional changes both before and after inhibitor administration, closely matching empirical data of Ref.^3^. Adenosine inhibitor vaccine is modeled in (6). This vaccine reduces the level of adenosine. We will now consider two simple ODEs that represent the concentration of adenosine and hypoxia inhibitor vaccines:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{dpx}{dt}=v\left(t\right)- \eta px,$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{ds}{dt}= {v}_{s}\left(t\right)-\lambda s.$$\end{document}

In (5), px represents the normalized concentration of hypoxia inhibitor (unitless), and η is the drug elimination rate. The parameter η is estimated to match experimental data by considering pre- and post-vaccine hypoxia levels^4^. In (6), s represents the normalized concentration of adenosine inhibitor (unitless), and λ is its elimination rate. The parameters of the adenosine inhibitor drug are estimated based on experimental data, where the drug levels were measured at different time points following injection. For further information, please refer to Fig. 1 in Ref.^3^. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v\left(t\right)=40\delta (t-{t}_{v})$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${v}_{s}\left(t\right)=4200\delta (t-{t}_{vs})$$\end{document} are the vaccination functions for hypoxia and adenosine inhibitor drugs respectively. tv and tvs are the vaccination time points for hypoxia and adenosine inhibitor drugs respectively. For the sake of simplicity, we have omitted the inclusion of relationships to simulate the diffusion of these drugs within the tumor environment, assuming a uniform concentration across all grid cells at each time point.

Adenosine inhibitor drug has adhesive properties and accumulates in the tumor area^3^. Therefore, a moving average filter is considered to model that properly by:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${S}_{t}\left[n\right]=\frac{1}{M}{\sum }_{k=0}^{M-1}s[n-k].$$\end{document}

Equation (7) means that a moving average filter is applied to the output of (6). In (7), M is the number of points used for averaging (see the output of (5) and (7) in Fig. A2 in Appendices).

Adenosine and hypoxia have been shown to act as pro-tumor factors through various mechanisms. One of these mechanisms is the promotion of Treg activation^3,4^. Therefore, we will consider the following simple ODE to model the probability of Treg activation rate:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\overline{P} }_{Treg}={P}_{Treg}+\frac{{r}_{Hr}}{U+c}{+ r}_{Ar} Ad,$$\end{document}

where the production rate of Treg cells \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({\overline{P} }_{Treg}$$\end{document} ) increases with decreasing oxygen concentration or rising adenosine at rates rHr and rAr, respectively. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${P}_{Treg}$$\end{document} is the base value of Treg activation rate.

Secondly, adenosine and hypoxia have been shown to directly promote the proliferation rate of tumor cells^3,4,7,9^:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\overline{P} }_{Tumor}={P}_{Tumor}+\frac{{r}_{Ht}}{U+c}{+ r}_{At} Ad,$$\end{document}

where* r_At_ and r_Ht*_ are related to the effect of adenosine and hypoxia on the production rate of tumor. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${P}_{Tumor}$$\end{document} is the base value of tumor proliferation rate.

Thirdly, adenosine and hypoxia have been found to decrease the cytotoxicity rate of effector cells^3,4^:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\mu }=\mu - \frac{{r}_{He}}{U+c}{- r}_{Ae} Ad.$$\end{document}

Similarly, in (10) rAe and rHe are related to the effect of adenosine and hypoxia on the effector cells cytotoxicity rate. µ is the base value of effector cell cytotoxicity.

The last step is DC modeling. VDC number of active DCs was randomly distributed in the grid on the seventh day (the injection day). This number was selected such that the model output corresponds to the experimental results.

Numerical simulations

Set-up of numerical simulations

We consider a 210 × 210 grid for our model. The size of the grid is chosen so that there is enough space for the tumor to grow until the 18th day (the last day of experimental data). The number of DCs and lymphocytes in the body has been extensively studied in various references. However, it is important to note that these ratios can be altered when cancer is present, and the variations can depend on the specific characteristics of the immune system^34^. Therefore, when selecting the cell counts for this model, we consider the available data set and the desired output of the model. As a result, on the first day, we randomly distribute 100 inactive DCs and 750 inactive T lymphocytes across the entire grid.

The experimental data used in this study includes the percentage of Treg cells in various groups with and without treatment. In order to convert the unit of this data into the number of Treg cells, a coefficient was applied. This coefficient was calculated to ensure that the variable representing Treg cell count matched between the simulation and experimental data for the untreated group. The total experimental data pertaining to Treg cells was then multiplied by this coefficient.

The average diameter of various tumor cells is reported in Refs.^35–37^. Based on these studies and to reduce complexity in the simulations, we assume that each grid-site is 10 µm apart in both the x and y directions. Therefore, the size of each grid-site is 10^–4^ mm^2^. To determine the tumor volume, the number of sites occupied by the tumor is counted and then multiplied by 10^–4^ mm^2^. To convert the unit of tumor size from the number of cells to mm^3^, we multiply it by 2200. We selected this factor based on the observation that significant volumetric expansion occurs in three dimensions compared to two dimensions and one dimension, as demonstrated in Ref.^38^. This study measured and compared tumor volume when cells were arranged in one dimension (on a line), two dimensions (in a plane), or in a three-dimensional space, and a significant increase in volume in three dimensions, compared to two-dimensional and one-dimensional found. Therefore, we chose this factor so that the model and real outputs in the first step of the simulation in the untreated group match. We chose a time step of Δt = 2 min to simulate our model (small enough to capture the dynamics of immune cell movements).

We hold the boundaries of the grid at a fixed oxygen concentration (0.8 mMol) and consider a zero initial condition for the concentration of oxygen inside the grid. It is assumed that there are blood vessels on the borders and they diffuse oxygen into the grid (see the simulation steps of the model in Fig. A4 in Appendices). The range of some model parameters has been determined based on previous studies. By utilizing experimental data, these parameters have been adjusted to ensure that the outputs of the model closely match the experimental data. The estimated parameters are listed in Table 1. For a visual representation of a sample simulation in the lattice of the proposed model, please refer to Fig. A1 in Appendices. The code pertaining to siRNA + Px-478 + DC group is available on GitHub at: https://github.com/elaheghiyabi96/Modeling-Tumor-Escape-Mechanisms-for-DC-Based-Immunotherapy-Insights-into-the-Role-of-Hypoxia-and-A.Table 1. The parameters in the model, their description and their values estimated from experimental data.ParameterDescriptionValueUnitSource1 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} Effect of the vaccine on adenosine0.000038Δt^−1^Estim2 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a$$\end{document} Effect of hypoxia on adenosine0.9412Mol. cells^−1^–3 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r$$\end{document} Adenosine elimination rate0.07Δt*^−1^Estim4 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} Constant rate of adenosine production12Δt^−1^Estim5r_AtEffect of adenosine on tumor0.000025Δt^−1^–6rHtEffect of hypoxia on tumor0.0021Mol. cells^−1^. Δt^−2^Estim7rArAdenosine effect on Treg0.0015Cells. Δt^−1^Estim8rHrHypoxia effect on Treg0.0047Mol. cells^−2^. Δt^−2^Estim9 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document} The elimination rate of adenosine inhibitor vaccine0.002Δt^−1^Estim10DactDC activation rate0.003Cells. Δt^−1^From^25^ with change11CactT lymphocyte activation rate0.24Cells. Δt^−1^From^25^ with change12µEffector cell cytotoxicity rate0.45Cells. Δt^−1^From^25^ with change13rAeEffect of adenosine on effector cell cytotoxicity rate0.001Δt^−1^Fit to data14rHeEffect of hypoxia on effector cell cytotoxicity rate0.014Mol. cells^−2^. Δt^−2^Fit to data15DinactDCs inactivation rate by Treg0.01Cells. Δt^−1^Fit to data16 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta$$\end{document} Elimination rate of hypoxia inhibitory vaccine0.014Δt^−1^Estim17KThe effect of the vaccine on oxygen0.00004Mol. cells^−1^. Δt^−1^Estim18PTregTreg cell production rate0.3Cells. Δt^−1^Fit to data19PtumorTumor growth rate0.001Δt^−1^^25^20 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\overline{f}}$$\end{document} Rate of oxygen consumption by the tumor cell0.0023Mol. cells^−1^. Δt^−1^From^22^ with change21cTo avoid division by zero in the denominator0.01Mol. cells^−1^. Δt^−1^–22LWalk exponent1.15–^25^23NTT-lymphocyte count750––24NDInactive DC count100––25VDC_DC vaccine70––26Mnumber of points in moving average800––27Nsize of grid points in both the x and y directions210––28Χsize of each grid-site in both the x and y directions10μm^25^Δt* = 2 min.

Model assessment

We investigate the dynamics of different variables in the model in 8 groups (see Table 2). Adenosine and hypoxia inhibition have been investigated in two different studies using siRNA-loaded chitosan lactate nanoparticles (siRNA for short) and Px-478 vaccines^3,4^. In both studies, tumors were established by injection of 7 × 10^5^ 4T1 breast cancer cells subcutaneously into the right flank of female BALB/c mice (6–8 weeks old) on day 0, and tumor size (mm^3^) was assessed every 2 days after that. Also, both works used DC-based immunotherapy in combination with inhibition of adenosine or hypoxia. The DC vaccine was injected into mice on day 7, and the tumor size is available up to day 18 in both works. In addition to tumor size, some other variables such as the number of Treg cells, cytotoxicity of effector cells, or the concentration of vaccines in tumor area were also measured on specific days.Table 2. Different groups in experimental and computational investigations of this work.Group namePx-478 (Days = 6, 8, 10, 13, 15)siRNA (Days = 6, 9, 12, 15)DC (Day = 7)1Untreated ×  ×  × 2Px-478✓ ×  × 3siRNA × ✓ × 4DC ×  × ✓5Px-478 + DC✓ × ✓6siRNA + DC × ✓✓7siRNA + Px-478✓✓ × 8siRNA + Px-478 + DC✓✓✓

Since our model is stochastic, we repeat each simulation 30 times and store the results, and then report the average and 95% confidence interval of stored results. For a quantitative comparison, we report R-squared (RSQ) for each group. The RSQ is calculated to assess how well the model output and the empirical data match and it has a value between 0 and 1 (the larger it is, the closer the model output and the empirical data are). In this study, our primary objective is to investigate tumor behavior during the initial days following tumor establishment, with a specific emphasis on analyzing the impact of vaccines using experimental data. Given that extending simulations may not yield additional insights relevant to our research question, our simulations are confined to the same timeframe as the experimental data collection (up to 18 days post tumor establishment).

Figure 3 shows the volume of the tumor in different groups. In the plots on this figure, the squares show the empirical data and the dashed lines show the confidence interval. According to these plots, the groups that had combined treatment are more effective than the groups that had single treatment, and the most effective treatment is the combined treatment with all vaccines. There was no experimental data for the group’s siRNA + Px-478 and siRNA + Px-478 + DC, and conclusions were made based on the suggested model. The best RSQ corresponds to untreated group and the worst RSQ corresponds to treatment with siRNA + DC group. Although the value of RSQ has decreased in some groups, but by comparing the behavior of the model in different groups, it can be concluded the model has correctly predicted the improvement or deterioration of the results in different groups. For a more detailed comparison and analysis of the results obtained from the model and experimental data in different groups, please see Fig. A3 in the Appendixes.Figure 3. Tumor volume was predicted by the model and measured in the empirical data^3,4^ in different groups. (a) Untreated. (b) Treatment with adenosine inhibitor. (c) Treatment with DC vaccine. (d) Treatment with adenosine inhibitor and DC vaccine. (e) Treatment with hypoxia inhibitor. (f) Treatment with hypoxia inhibitor and DC vaccine. (g) Treatment with adenosine and hypoxia inhibitor. (h) Treatment with adenosine and hypoxia inhibitor and DC vaccine.

Treg/effector cell ratio correlates with protective immunity in mice and humans^39^. To investigate this ratio, Fig. 4 shows the number of effector cells and Tregs in four groups. As we wanted to see the effect of only adenosine and hypoxia on effector and Treg cell numbers, the groups that received the DC vaccine were not considered here. According to Fig. 6, in the siRNA and Px-478 groups (Fig. 4b,c), the difference in the number of Treg and effector cells has decreased compared to the untreated group. Because adenosine and hypoxia have been suppressed, immune suppressor cells have reduced in siRNA and Px-478 groups. However, when the outcomes are compared, we can find that the adenosine inhibitor has enhanced the Treg/effector cell ratio more than the hypoxia inhibitor. In the Px-478 + siRNA group, the number of effector cells has increased more than the number of Treg cells (Fig. 4d). Therefore, the outputs demonstrate that combination therapy is better than other groups in improving immune system cell performance.Figure 4. The number of effector cells and Treg cells in different groups. (a) Untreated. (b) Treatment with adenosine inhibitor. (c) Treatment with hypoxia inhibitor. (d) Treatment with adenosine inhibitor and hypoxia inhibitor.

Figure 5 compares the oxygen concentration and tumor cells distribution in the untreated and the siRNA + Px-478 + DC groups on day 12. The oxygen concentration in siRNA + Px-478 + DC group has increased compared to the untreated group in the surrounding areas of the tumor, and as a result, hypoxia decreased. Also, the tumor volume has decreased in the siRNA + Px-478 + DC group.Figure 5. Comparison of the oxygen concentration and tumor cells in the untreated and siRNA + Px-478 + DC groups on day 12. (a) Oxygen concentration in the untreated group on day 12. (b) Oxygen concentration in the siRNA + Px-478 + DC group on day 12. (c) Tumor cells in the untreated group on day 12. (d) Tumor cells in the siRNA + Px-478 + DC group on day 12.

Figure 6a shows the dynamic of µ (effector cell cytotoxicity rate) in different treatment groups. The fluctuation of the output in siRNA and siRNA + Px-478 groups of this Figure is due to the periodic injection of the vaccine. As is seen in this figure, in the untreated group, tumor development and increasing hypoxia and adenosine have decreased cytotoxicity rate over time. In the Px-478 and siRNA groups, the cytotoxicity rate has increased. Hypoxia and adenosine reduce the effector cell cytotoxicity rate. Therefore, their inhibiting has prevented the reduction of this rate. In comparison to others groups, the cytotoxicity rate in the Px-478 + siRNA group (combination therapy) is higher. Therefore, inhibition of hypoxia and adenosine has improved the function of effector cells.Figure 6. Investigating the time course of effector cell cytotoxicity rate and adenosine concentration in four groups. (a) Effector cell cytotoxicity rate, (b) the adenosine concentration.

Figure 6b shows the level of adenosine in different groups. Adenosine has begun to rise more quickly in the untreated group after the 14th day, because hypoxia in the TME increases as the tumor grows. Adenosine concentration in the TME also rises when hypoxia increases. In the Px-478 group, suppression of hypoxia prevents the acute increase in adenosine. In Fig. 6b the results of the siRNA and siRNA + Px-478 groups are identical. Because according to the outputs of tumor volume and immune system cells, adenosine has a greater impact on tumor growth. As a result, the tumor volume is reduced when adenosine is inhibited. As tumor volume decreases, hypoxia also decreases and a low amount of hypoxia cannot have a substantial effect on the growth of adenosine.

The effect of changing the injection dose and time

The proposed model may be used to monitor changes in tumor volume by varying the time, number, and dose of vaccines. We have observed that the siRNA + Px-478 + DC group yields the most favorable therapeutic outcomes (see Fig. 3h). The results of altering the injection patterns within this group are presented in Table 3. In the first type of these patterns, only the dose and number of injections of the DC vaccine were altered, and all other vaccinations remained unchanged compared to the siRNA + Px-478 + DC group. The results of this type emphasize how crucial it is to choose the correct dose and number of injections. It can be seen that the sooner the vaccine is injected, the more effective it is and the less the tumor grows. An interesting result of this type is that increasing the dose and number of injections may sometimes unexpectedly lead to an increase in tumor size. Because the increase of active DCs may cause traffic in the area around the tumor, preventing effector cells from reaching the tumor cells. To investigate this effect, the average number of effector cells surrounding the tumor was measured for 18 days, and the results are displayed as a function of injections dose in Fig. 7. As seen in this figure, increasing the number of injected cells decreases the access of effective cells to tumor cells. This is also reported in other experimental^40^ and computational^25^ studies.Table 3. Investigating the effect of changes in the dose, the time of the first injection, and the time interval between injections on the final and the average volume of the tumor.Type of vaccineType of changeAverage tumor volume (mm^3^)Tumor volume on day 18 (mm^3^)(1) DC vaccineTime of initial injection (one injection only)★Day 60.00050Day 721.564.25Day 835.23102.69Injection dose (injection on day 7)30%↓25.474.230%↑16.941.850%↑21.164.35The number of injections (injection starts on day 7 and other injections with an interval of one day)2 times15.1643.164 times26.968.96 times33.9785.77The number of injections (injection starts on day 7 and other injections with an interval of two day)2 times29.4582.474 times28.3780.66 times25.0466.28(2) siRNA vaccineChanging the days of siRNA vaccine injection (other vaccines are unchanged)Everyday9.8916.47Every two days19.9344.55Every three days22.2259.1Changing the dose of siRNA vaccine (other vaccines are unchanged)30%↓41.19121.8330%↑24.4355.6350%↑16.9639.68(3) Px-478 vaccineChanging the days of Px-478 vaccine injection (other vaccines are unchanged)Everyday16.741.8Every two days20.161.20Every three days31.5769.41Changing the dose of Px-478 vaccine (other vaccines are unchanged)30%↓27.3677.2430%↑24.658.950%↑21.2858.87(4) Px-478 and siRNA vaccineChanging the dose of Px-478 and siRNA vaccines and the time of injection of these vaccines every day30%↑3.02050%↑2.6070%↑1.20Changing the dose of Px-478 and siRNA vaccines and the time of injection of these vaccines every two days30%↑14.9534.56★50%↑2.64.670%↑2.93.36Changing the dose of Px-478 and siRNA vaccines and the time of injection of these vaccines every three days30%↑16.4847.9650%↑13.5831.1370%↑9.921.53Injection of Px-478 and siRNA vaccines every other dayStart with Px-47821.4248.37★Start with siRNA8.618.43Figure 7The average number of effector cells around the tumor as a result of changing the dose of the DC vaccine.

In the second to forth part of Table 3, different patterns for the injection of hypoxia and adenosine-inhibiting vaccines have been simulated. The outcomes demonstrate that increasing the injection volume and decreasing the time between injections can lead to improved treatment. But this improvement may be insignificant in higher doses or higher number of injections. When the vaccines are given every other day, the result is more favorable when the adenosine vaccine is injected first. Also, the results show that adenosine has a greater effect on tumor growth than hypoxia. According to the results of Table 3, only a few of the tested scenarios significantly reduced the tumor volume and are more effective than the others: (1) injection of DC vaccine on day 6, (2) injection of siRNA and Px-478 vaccines once and every two days with a 50% dose increase, (3) injection of Px-478 and siRNA vaccines every other day when the siRNA vaccine is injected first.

Global sensitivity analysis

Global sensitivity analysis calculates how sensitive the model outputs are to the model parameters. Finding the model-sensitive parameters enables us to identify the crucial and influencing agents in TME. Then, using these important factors, the model and treatment techniques may be improved and expanded. Unlike local sensitivity analysis, global sensitivity analysis involves changing every model parameter at once. Since all input parameters can be changed at once, the global sensitivity analysis is more realistic and ideal for nonlinear input–output relationships. The sensitivity analysis results are shown in Fig. 8. The output considered in this analysis is the average tumor volume during 18 days (please refer to the sensitivity analysis method provided in the Appendices). According to this figure, the most sensitive parameters are β, rAt, and r. These parameters are related to the constant rate of adenosine production, the effect of adenosine on tumor growth and, the rate of adenosine elimination, respectively. This emphasizes that in this model, adenosine has a greater impact on output than hypoxia. Also, in the preceding section, the effects of the adenosine inhibitor vaccine on tumor volume were more visible than hypoxia.Figure 8. Sensitivity value of the model parameters.

According to the sensitivity analysis, β (constant rate of adenosine production) is the most sensitive parameter of the model. To investigate how this parameter affects the variables of the model, the time evolution of the tumor, the killing rate of the effective cells, adenosine, and number of Treg for four values of β are plotted in Fig. 9. By just raising the constant rate of adenosine production, the level of adenosine, tumor volume, and number of Treg are significantly enhanced while the cytotoxicity rate is reduced. This demonstrates that adenosine plays an important role in tumor growth and immune system activity. Inhibiting adenosine or other factors that boost it, such as hypoxia, can be beneficial in treating cancer.Figure 9. The time course of tumor volume, effector cell cytotoxicity rate, adenosine concentration, and number of Treg as a result of changes in β (a) Tumor volume. (b) Effector cell cytotoxicity rate. (c) Level of adenosine. (d) Number of Treg*.*

The sensitivity of parameters PTreg and µ are also high. For further investigation, the tumor volume has been examined in connection to changes in these parameters in Fig. 10. According to the results of the figure, it is necessary to identify factors that affect the production and cytotoxicity rate of effector cells. Increasing the cytotoxicity rate of the effector cells and their number leads to a decrease in tumor volume, and according to Fig. 5, decreasing the tumor volume reduces hypoxia, and according to Fig. 6, the final level of adenosine will decrease.Figure 10. The time course of tumor volume as the parameters µ and PTreg change. (a) Tumor volume as a result of the change in µ. (b) Tumor volume as a result of the change in PTreg.

Identifiability analysis

Uncorrelated parameters may be regarded as the crucial adjusting knobs of a model^41^. Interdependence among parameters makes the model unidentifiable, which means their values cannot be determined uniquely with the information at hand. Thus, characterizing the subset(s) of identifiable parameters and how they interact with each other and model output is crucial for a model. To perform this analysis, as explained in Ref.^41^, correlation coefficients between pairs of parameters were computed for different outputs in the model (see the Identifiability Analysis Method provided in the Appendices). Tumor volume, adenosine, hypoxia, effector cell cytotoxicity rate, number of effector cells, Treg cells, and DCs were considered as model outputs. If the correlation coefficient of a pair of parameters was equal to or greater than 0.9 in all outputs^42^, the corresponding parameters were interdependent and one of them could be eliminated from the model. According to the results of this analysis (Fig. A5 in Appendices), there is no pair of parameters in the model whose correlation coefficient is greater than or equal to 0.9 in all output, so, all the parameters in the model can be identified.

Discussion

Activities to develop new cancer therapies have accelerated during the past several decades. Cancer immunotherapy is one of these developments and has significantly increased the survival chances of patients^43^. Mathematical modeling has been beneficial in improving our understanding of cancer and its various treatment methods. In this study, a spatio-temporal model for tumor growth and its interaction with the immune system has been presented. Also, hypoxia and adenosine were added to this model and their effects on TME and each other were examined. The model parameters were estimated using experimental data obtained from two distinct studies, investigating the effect of adenosine and hypoxia on DC-based immunotherapy, performed on 4T1 breast cancer-bearing mice. This model was used to investigate an effective treatment pattern. Model simulations predict that by applying this pattern, the tumor volume will decrease and the immune system will perform better than other groups.

The simulation results of this work show that adenosine grows strongly when hypoxia increases in the TME, and inhibiting hypoxia leads to a decrease in adenosine levels. Also, experimental studies confirm the adenosine increase due to the oxygen reduction^13,16^. In addition, in these studies, it has been shown that the increase in hypoxia leads to an increase in suppressor cells. Therefore, it can be concluded that an increase in hypoxia leads to an increase in adenosine, and according to Fig. 10, the number of suppressor cells increases with an increase in adenosine. As a result, hypoxia leads to the increase of suppressor cells both directly and indirectly through adenosine. Also, using sensitivity analysis, it was confirmed that adenosine has the greatest effect on tumor growth in TME compared to other model factors. Therefore, to achieve better therapeutic results, adenosine must be inhibited, either directly (with its inhibitory vaccine) or indirectly (for example, by inhibiting hypoxia).

Our model showed that the sooner the vaccine is injected, the more effective it is and the less the tumor grows. In some studies, it has been reported that a high tumor load may favor the initial expansion of immune cells, but on the other hand, it may also enhance the suppression of tumor immunity^44^. Thus, the optimal approach to using immunotherapy in conjunction with surgical resections is a question that can be addressed using model simulation with different scenarios for initial tumor size and immune system cells.

One of the challenges in vaccine-based treatments is finding the optimal dose and interval between injections. We used the proposed model to determine these factors. The results showed that increasing the dose or the number of injections does not always improve the treatment. For the DC vaccine, it was seen that there is a certain threshold beyond which increasing the vaccine dose not only does not reduce tumor growth but even increases it. This is because with increasing the number of injected DCs, the space around the tumor is completely occupied by DCs, preventing effector cells from getting to the tumor cells. Some studies also emphasize that in order to prevent the traffic of DCs and achieve successful treatment in DC immunotherapy, the appropriate and optimal number of DCs should be used^23,25^. In addition to the number, the location of the DCs and how they move in the network is important. Because the closer the DCs are to the T lymphocytes, the more likely the T lymphocytes will be activated, and this will improve cancer treatment. As a result, the number of DC cells, their locations in the grid, and their potential to activate T lymphocytes can all affect the tumor volume. By using spatial models, we have gained insights into the behavior of DC cells that temporal models, such as ODEs, have been unable to capture. This underscores the importance of incorporating spatial dynamics into our modeling approach in order to achieve a more comprehensive understanding of cellular phenomena.

In vaccination with inhibitors of hypoxia and adenosine, tumor growth is slowed by increasing the injected dose or the number of injections or reducing the time interval between injections. Adenosine inhibitory vaccine has more effect on the final volume of the tumor than the hypoxia inhibitory vaccine. The model was used to evaluate several injection strategies and identify the near-optimal injection pattern. For the inhibition of adenosine and hypoxia, three near-optimal strategies were found using a trial-and-error manner (indicated with ★ symbols in Table 3).

The sensitivity analysis of the model showed that the first three sensitive parameters of the model are β, rAt, and r, which are all related to the adenosine dynamic. According to all the mentioned results, it can be concluded that adenosine is more involved in the tumor growth than hypoxia. Another study has also emphasized that among suppressive factors, it seems that adenosine exerts higher immunosuppressive function as its producing molecules including CD73 are expressed on both tumor and immune cells^3^. Additionally, µ and PTreg (effector cell cytotoxicity rate and Treg production rate, respectively) are also important parameters, which are also emphasized in other studies^25^. Targeting these parameters can be considered in designing treatment strategies.

We can mention several limitations for the presented model to be considered in future works. The parameters of this model are estimated for a time frame of 18 days. Measuring tumor size over a longer period of time can better show the effect of hypoxia on the tumor. Because as the tumor grows, its metabolism elevates and the need for oxygen increases, which causes hypoxia. Investigating the model in larger time frames will force us to consider other phenomena such as angiogenesis (something that is usually not seen in the early days of tumor formation). Formation of blood vessels near the tumor can improve oxygen delivery to it and have a direct effect on hypoxia dynamics. In this work, we considered the effect of hypoxia on adenosine, but the effect of adenosine on hypoxia was not investigated. The feedback that these two factors have on each other can be a determining factor in the dynamics observed in the experimental data. In the proposed model, we considered that the tumor has a single layer, and therefore, all tumor cells consume oxygen similarly. The tumor can be considered as a three layers mass (as is considered in Ref.^22^), in which the inner layer is necrotic, the middle layer is quiescent and the outer layer is proliferating. These three layers have different oxygen consumption rates, and by considering them, more accurate dynamics of oxygen concentration in the TME can be obtained. We have endeavored to estimate the maximum number of parameters using our experimental data sourced from Refs.^3,4^. Notably, the data from both^3,4^ were derived from the same mouse models and tumor type, rendering them highly compatible. For parameters not directly available from these references, reasonable initial values were gleaned from other studies and fine-tuned within the model. It is conceivable that with more extensive experimental data, all model parameters could be estimated from a single dataset, thereby reducing uncertainty within the model. Finally, expanding the model and adding details of the tumor growth, cancer-promoting factors, and its interaction with other immune cells can help our understanding of the development and progression of the tumor and the use of targeted treatments, which will lead to improving current immunotherapies.

Conclusions

Overall, in this study, we presented a spatio-temporal model to investigate the role of hypoxia, adenosine, and DC cells in cancer. The simulations showed that these factors affect the proliferation of cancer and hypoxia also affects the level of adenosine. We also concluded that adenosine has a greater effect than hypoxia in TME. Also, the simulations showed that the efficacy of DC-based immunotherapy is strongly dependent on the injection rate and the trafficking of DCs in TME. In addition, we introduced near-optimal treatment protocols and showed that the tumor volume is significantly reduced with these protocols.

Supplementary Information

Supplementary Information.