Bistable dynamics of TAN-NK cells in tumor growth and control of radiotherapy-induced neutropenia in lung cancer treatment

1.   Introduction

Lung cancer is the leading cause of cancer-related mortality worldwide ($\sim$ 1.8 million (18%) deaths each year ^[1]). The most common type ($\sim 85\%$) is non-small cell lung cancer (NSCLC) with common subtypes of lung squamous cell carcinoma (LUSC) and lung adenocarcinoma (LUAD) ^[2]. The tumor microenvironment (TME) includes various cell types such as immune cells (tumor-associated neutrophils (TANs), macrophages, T cells, natural killer cells) and stromal cells, extracellular matrix (ECM), and diffusible molecules such as cytokines, chemokines, and growth factors ^[3,4,5]. Certain types of immune cells in TME were reported to support tumor growth by suppressing immune activities ^[6] and exchanging tumor promoting factors such as angiogenic factors, ECM-degrading enzymes, and growth factors. It was shown that high degrees of tumor-infiltrating lymphocytes (TILs) in the large node-negative cancer tissue are positively associated with a decrease in disease recurrence and improved survival rates in NSCLC patients ^[4,7]. In particular, some of these cells show distinct functional characteristics in tumor-promoting or tumor-suppressing roles. For example, it is well established that tumor-associated macrophages (TAMs), termed M2 TAMs ^[6], promote tumor growth by secreting factors and cytokines ^[8,9]. Even though a classical form of neutrophils called N1 TANs present anti-tumor functions in the TME, increasing evidence now suggest that other neutrophils, called N2 TANs, can mediate aggressive tumor growth throughout the course of tumor progression toward malignancy ^{[10,11,12,13,14,15]}. For example, TANs infiltrate tumor tissue and enhance tumor growth and cancer cell invasion ^[16,17], leading to poor clinical outcomes ^[16,17] in many cancers, including renal carcinoma ^[18], metastatic melanoma ^[19], and bronchoalveolar carcinoma ^[17]. Clinically, the increased neutrophil to lymphocyte ratio (NLR) is positively correlated with poor prognosis in various cancers ^[20], including NSCLC ^[21], esophageal cancer ^[22], colon cancer ^[23,24], head and neck cancer ^[25], gastric cancer ^[26], hepatocellular carcinoma ^[27], bladder cancer ^[28], cervical cancer ^[29], or advanced pancreatic cancer ^[30]. The NLR is considered a potentially crucial biomarker for lung cancer patients ^[31]. Despite the significant potential for targeting these TANs ^[32], the detailed mechanism of this critical transition between those modes is poorly understood. Therefore, a better understanding of tumor-associated inflammation may help us identify novel therapeutic strategies ^[33].

The signal activator and transducer of transcription (STAT) protein family functions in a coordinated fashion to control key cellular processes such as cell division, cell growth, cell migration, and cell death ^[34,35,36]. The STAT family consists of seven members: STAT1, STAT2, STAT3, STAT4, STAT5 (STAT5A and STAT5B), and STAT6 ^[37]. Insufficient STAT1 signaling causes a defect in the response to external control modules in diseases such as cancers. For example, the lack of STAT1 promotes the migration and spread of lung cancer cells ^{[37,38,39,40]}. STAT3 plays an important role in a variety of biological functions, including inhibition of apoptosis, cell motility, and cell growth ^[37,38,41]. The relative expression of STAT1 and STAT3 may induce two cell fates: (ⅰ) An anti-apoptosis (down-regulation of STAT1 and up-regulation of STAT3) status and (ⅱ) an apoptosis (up-regulation of STAT1 and down-regulation of STAT3) status. Apoptosis is regulated by two key proteins downstream of STAT signaling, Bcl-2 and BAX, that inhibit or promote apoptosis, respectively ^[42,43,44]. The reduced DNA binding capability of STAT3 predates dynamic changes in the expression of anti-apoptotic Bcl-2 and pro-apoptotic BAX proteins (decreased Bcl-2 level and increased BAX levels) and induction of apoptosis ^[45]. How these two phenotypes are controlled in detail is largely unknown, and many experimental and clinical results suggest significant potential for therapeutic purposes.

Transforming Growth Factor-$\beta$ (TGF-$\beta$) is a multifunctional cytokine that plays a critical role in the regulation of cell growth, differentiation, apoptosis, and immune function. TGF-$\beta$ is involved in a wide range of biological processes, including embryonic development, wound healing, tissue repair, and immune responses ^[46,47,48]. TGF-$\beta$ plays a complex role in various cancers, including lung cancer. Even though TGF-$\beta$ acts as a tumor suppressor by inhibiting the growth and spread of cancer cells in an early stage of cancer progression, TGF-$\beta$ can promote the growth and spread of cancer cells and facilitate tumor invasion and metastasis in later stages of the disease ^[49,50,51]. In lung cancer, TGF-$\beta$ signaling pathways are often disrupted, which can result in the loss of the tumor suppressive effects of TGF-$\beta$ and the acquisition of pro-tumorigenic functions ^[46,48]. TGF-$\beta$ has also been implicated in promoting immune evasion in lung cancer. It can suppress the activity of immune cells, including T cells and natural killer (NK) cells, leading to reduced anti-tumor immune responses and metastasis of lung cancer ^[52,53]. TGF-$\beta$ was also shown to skew neutrophil differentiation toward N2 TANs ^{[15,54,55,56]}, and inhibition of TGF-$\beta$ can reverse this effect ^[57,58]. Therefore, the development of TGF-$\beta$ inhibitors as a potential therapeutic strategy in lung cancer is an active area of research.

Interferon gamma (IFN-$\gamma$) is a cytokine that plays an important role in the immune response to cancer. In lung cancer, IFN-$\gamma$ has been shown to have both anti-tumor and pro-tumor effects, depending on the stage of the disease and the context of its action ^[59,60]. In lung cancer, IFN-$\gamma$ can stimulate anti-tumor immune responses by activating immune cells, such as T cells and natural killer (NK) cells, which can directly kill cancer cells ^[61,62]. IFN-$\gamma$ also enhances the expression of major histocompatibility complex (MHC) molecules on the surface of cancer cells, making them more visible to the immune system. Researchers ^[10,54,63] suggest that IFN-$\gamma$ can stimulate the production and activation of N1 TANs in the TME, leading to enhanced anti-tumor immunity. In addition, IFN-$\gamma$ has been shown to up-regulate the expression of chemokines that recruit N1 TANs to the tumor site ^[10,54,63]. Studies have also shown a positive correlation between IFN-$\gamma$ and activation of STAT1 in lung cancer, leading to improved outcomes and increased sensitivity to chemotherapy and immunotherapy, and recruitment of N1 TANs ^[39,64,65].

NK cells are a type of lymphocyte that play an important role in the immune response to cancer. NK cells have been shown to have both anti-tumor and pro-tumor effects, depending on the stage of the disease and the context of their action in the tumor ^[52,66,67]. In lung cancer, NK cells are thought to play a protective role by recognizing and eliminating cancer cells. Studies have shown that higher levels of NK cell infiltration in lung tumors are associated with improved patient survival rates ^[68,69,70]. There is a positive correlation between NK cells and cytokines, particularly IFN-$\gamma$, which is a key activator of NK cell function ^[61,62]. IFN-$\gamma$ can induce the expression of receptors on NK cells that recognize and kill cancer cells, as well as enhance their production of cytokines that stimulate other immune cells to target cancer ^[59,60]. TGF-$\beta$, on the other hand, has been shown to suppress the activity of NK cells, which can contribute to immune evasion and promote the aggressive growth and spread of cancer cells ^[71,72]. Thus, the relative balance between the effects of IFN-$\gamma$ and TGF-$\beta$ on NK cell activity can be an important factor in determining the outcome of the immune response to lung cancer. Exogenous NK cells administered to a patient as a form of immunotherapy are being studied as a potential treatment option for various cancers, including lung cancer ^{[73,74,75,76]}. In oncolytic virus therapy, injection of exogenous NK cells or depletion of NK cells in the TME were shown to improve tumor control compared to the control case where NK cells can decrease anti-tumor efficacy of OV-mediated cancer cell killing ^[74,76].

Radiation therapy (RT), one of the most common types of cancer treatment, uses high-energy X-rays to destroy cancer cells. While various types of RT such as 3D conformal radiation therapy, image guided radiation therapy (IGRT), intensity modulated radiation therapy (IMRT), and volumetric modulated arc therapy (VMAT) are used treat cancers, these may also affect normal healthy cells with different sensitivity ^[77,78]. In particular, ionizing radiation (IR) kills cancer cells by destroying double-stranded DNA of tumor cells ^[79,80]. Studies indicate a significant effect of IR on the dynamics of immune cells such as such as T cells, macrophages, dendritic cells, and NK cells ^[79]. In particular, IR induces the secretion of IFN-$\gamma$ and tumor necrosis factor alpha (TNF-$\alpha$) by NK cells through the ATM and NK$\kappa$B pathway, resulting in anti-tumorigenic and apoptosis phenotypes ^[79]. IR-induced damage to tissues and cells is also shown to promote secretion of TGF-$\beta$ by fibroblasts and N2 TANs ^[81,82,83], leading to up-regulation of N2 TANs ^[15].

The complex network of interactions between tumor and the immune system include NK cells, N2/N1 TANs, STAT signaling, and immune cytokines in response to radiation therapy (Figure 1). In this work, we employed the framework of the mathematical model to investigate the key role of NK cells and radiation in regulation of TAN transitions and the STAT signaling pathway in the TME. The mathematical model consists of a system of ordinary differential equations (ODEs) and partial differential equations (PDEs) in the extended model (Section 3.7) involving the following variables at time $t$:

The mathematical model in this work is used (ⅰ) to investigate how NK cells regulate the critical balance between N1 and N2 TANs via TGF-$\beta$ and immune molecules (IFN-$\gamma$ and TGF-$\beta$ inhibitor such as Galunisertib, which blocks TGF-$\beta$ activity) in a complex interaction network, (ⅱ) to investigate how this perturbed N2-to-N1 ratio affect anti-tumor efficacy, (ⅲ) to investigate the role of intracellular signaling (STAT signaling) in regulation of the NK-controlling immune-tumor dynamics, (ⅳ) to identify the optimal schedule in a combination therapy (NK cells + IFN-$\gamma$ + TGF-$\beta$ inhibitor), and (ⅴ) to investigate how radiation therapy affect the overall complex nonlinear dynamics and to find the optimal dose schedule for the best outcomes in a clinic.

2.   Materials and methods

We developed a mathematical model based on a schematic diagram of key regulation of interaction among tumor cells, N1 TANs, N2 TANs, and NK cells in TME and intracellular STAT signaling within a tumor cell (Figure 1). The NK-controlled phenotypic balance between STAT signaling and N1/N2 TANs in the presence and absence of ionizing radiation regulates enhancement or suppression of tumor growth in response to TGF-$\beta$, IFN-$\gamma$, and TGF-$\beta$ inhibitor.

2.1.   Densities of N1 and N2 TANs: ($N_1(t), N_2(t)$)

To describe the time evolution of densities of N1 and N2 TANs, we took into account the transition from N1 TANs to N2 TANs via IL-6 and TGF-$\beta$ ^[10,84], induction of N1 TANs in the presence of IFN-$\gamma$ ^[10], natural decay processes of TANs, and the mutual inhibition between N1 and N2 TANs. The governing equations for densities of N1 and N2 TANs were modified from the model by Kim et al. ^[85]:

where $\lambda_6$ and $\lambda_{G}$ are the source of N2 TANs from IL-6 and TGF-$\beta$, respectively, from the TME, $\lambda_{S_1}$ is the source of N1 TANs from IFN-$\gamma$ from the TME, $\alpha$ is the inhibition strength of N2 TANs by N1 TANs with scaling parameters ($k_1, k_2$), $\beta$ is the inhibition strength of N1 TANs by N2 TANs with scaling parameters ($k_3, k_4$), and $\mu_{N2}$ and $\mu_{N1}$ are the decay rates of N2 TANs and N1 TANs, respectively. In this framework, the relative balance between IFN-$\gamma$ and TGF-$\beta$ determines either N1- or N2-dominant TME, which in turn either suppresses or promotes the tumor growth.

2.2.   Density of tumor population ($C(t)$)

Different kinds of mathematical models of tumor growth were introduced in literature to reproduce empirical time series data of tumor size ^[86]: Exponential growth ^[87], logistic growth ^[88], Gompertz growth ^[89], and nonlinear models ^[90]. In particular, logistic growth patterns in the presence and absence of supportive and/or inhibitive factors were observed in the experimental setting ^{[39,74,91,92,93,94]} and used in the various mathematical models (ordinary or partial differential equations, stochastic, or multi-scale models) for successful fitting to experimental data ^{[8,74,85,91,92,93,94,95,96,97,98,99,100,101,102,103,104]}.

The dynamics of the tumor cell density ($C(t)$) is governed by a production term ($P_{C}$) and death processes ($D_{C}$) as follows:

We take the logistic growth model for proliferation of tumor cells,

where $r$ is the proliferation rate of tumor cells and $C_0$ is the carrying capacity of the tumor in a given TME ($r, C_0 \in {\mathbb R}^+$).

We assume that tumor cells are killed by immune cells and IFN-$\beta$-induced apoptosis. In the model, N1 TANs and NK cells (both endogenous and exogenous NK cells, when injected) eliminate the tumor cells from the TME. Bcl-2 and BAX are key players in regulating the apoptotic death of cancer cells in the presence of various biochemical stimuli such as IFN-$\beta$ in the complex TME. As discussed below, the down-regulation of the gate keeper, Bcl-2 ($B < th_B$) ^[42], and up-regulation of BAX ($X > th_X$) ^[43,44] in response to high IFN-$\beta$ levels lead to cell death of tumor cells ^{[43,45,105,106]}. These lead to the following cell death terms

Here, $\delta_1$ and $\delta_3$ are the killing rate of tumor cells by N1 TANs and NK cells, respectively, $\delta_2$ is the elimination rate of tumor cells by the IFN-$\beta$-mediated apoptosis (STAT-Bcl-2-BAX signaling pathway), and $I[\cdot]$ is an indicator function (giving 1 or 0 based on the intracellular conditions on the tumor cell $\{(B, X) \in {\mathbb R}^2 :B < th_B, X > th_X\}$).

Combining these proliferation/death terms lead to the governing equation for the cancer cell density:

2.3.   Intracellular signaling network ($S_1(t), S_3(t), B(t), X(t)$)

Among several STAT family members ^[37], the relative balance between STAT1 and STAT3 plays a key role in regulation of various biological functions such as apoptotic status of given cell ^[37,38,41]. For instance, down-regulation of STAT1 and up-regulation of STAT3 lead to a stronger anti-apoptotic mode by maintaining the on-switch of the gate keeper, Bcl-2 ^[37,38,41], and can increase the invasive potential and aggressive spreading of lung cancer cells ^{[37,38,39,40]}. The mutual antagonism between STAT1 and STAT3 can reverse the gate keeping of the cell-death program. IFN-$\gamma$, the key inducer of the $N_2 \to N_1$ TAN transition ^[10,11,15], induces the apoptosis in malignant cancer cells ^[105] by upregulating STAT1 and down-regulating the STAT3. For instance, the reduced DNA binding capability of STAT3 predates dynamic changes in the expression of anti-apoptotic Bcl-2 and pro-apoptotic BAX proteins, inducing decreased Bcl-2 levels and increased BAX levels followed by apoptosis ^[45]. Thus, we assume that (ⅰ) IFN-$\gamma$ induce STAT1 at a rate $\lambda_{S_2}$ but inhibit STAT3 at a rate $\lambda_{S_3}$ with scaling parameters ($\lambda_1, K_1$), (ⅱ) all intracellular molecules, STAT1, STAT3, Bcl-2, and BAX, undergo natural decay at rates $\mu_{S1}, \mu_{S3}, \mu_{B}$ and $\mu_{X}$, respectively, and (ⅲ) mutual antagonism between STAT1 and STAT3 with inhibition parameters $\gamma, \delta$ and kinetic parameters ($a_1, a_2, a_3, a_4$). These biological observations (Figure 2A) and assumptions lead to the following dynamic systems of concentrations of STAT1 and STAT3 molecules within cancer cells (cf ^[107]):

Bcl-2 and BAX are downstream molecules whose molecular balance either inhibits or induce apoptosis ^[42,43,44]. While STAT1 induces apoptosis by inhibiting Bcl-2 (thus upregulation of BAX), STAT3 blocks the programmed cell death protocols by keeping Bcl-2 level intact. Thus, the STAT1-STAT3 signaling in Eqs (2.7) and (2.8) is connected to the following final steps of the apoptosis network ^[107]:

2.4.   Concentration of TGF-$\beta$ ($G(t)$) and TGF-$\beta$ inhibitors ($L(t)$)

TGF-$\beta$ is secreted by cancer cells and induces the phenotypic transition from N1 to N2 TANs ^[10,84]. Thus, TGF-$\beta$ inhibitors such as Galunisertib (LY2157299) can block this critical TGF-$\beta$-induced switch ^{[108,109,110]}, thus indirectly inducing N1-favorable TME conditions. Thus, the governing equations of concentrations of TGF-$\beta$ and its inhibitor are given by

where $\lambda_C$ is the secretion rate of TGF-$\beta$ by tumor cells, $\gamma_L$ is the degradation rate of TGF-$\beta$ by TGF-$\beta$ inhibitor, $\mu_G$ and $\mu_L$ are the decay rate of TGF-$\beta$ and TGF-$\beta$ inhibitor, respectively, and $u_L(t)$ is the external injection rate of the TGF-$\beta$ inhibitor.

2.5.   Concentration of IFN-$\gamma$ ($S(t)$)

IFN-$\gamma$ is secreted by stromal cells and immune cells, such as NK cells ^{[111,112,113,114,115,116]}, and can induce the phenotypic transition from the N2-mode to N1-mode ^[10,84] in the TME. IFN-$\gamma$ therapy alone or in a combination with other therapeutic drugs ^[60] has been reported to be effective in reducing tumor size ^{[117,118,119,120]} and in promoting the anti-tumor immune response ^[121]. In the model, we take into account the source, secretion by NK cells (both endogenous and exogenous NK cells), and natural decay. Thus, we get the governing equation:

where $u_S(t)$ is the time-dependent injection rate of IFN-$\gamma$, $\lambda_K$ and $\lambda_{K'}$ are the secretion rate of IFN-$\gamma$ by endogenous and exogenous NK cells, respectively, and $\mu_S$ is the decay rate of IFN-$\gamma$.

2.6.   Concentration of endogenous and exogenous NK cells ($K_{end}(t), K_{ex}(t)$)

NK cells are a type of immune cell that plays a critical role in the first line of defense against pathogenes and cancerous cells ^{[112,113,114]}. These immune cells are able to recognize and destroy infected or abnormal cells without the need for prior activation or antigen-specific recognition. NK cells have several mechanisms of killing, including the release of cytotoxic granules and the engagement of death receptors on target cells. In addition to their direct cytotoxic effects, NK cells also play a key role in regulating the immune response by producing cytokines such as IFN-$\gamma$ and TNF-$\alpha$. These cytokines not only stimulate other immune cells to mount an effective immune response but also have direct anti-tumor effects by inhibiting tumor growth and promoting the activation of other immune cells. However, these activities of NK cells may be effectively suppressed by other factors. For instance, TGF-$\beta$ can inhibit NK cell activation ^{[122,123,124]} by inhibiting T-bet and IFN-$\gamma$, and suppresses NK metabolism by downregulating mTORC1, SREBP, and c-MYC ^[125]. Therefore, injection of exogenous NK cells, cells that are generated outside of the body and then administrated to a patient as a form of immunotherapy, may be an alternative way of overcoming suppressed immunity with safety and efficacy in a typical TME ^[126]. For instance, Kim et al. and others showed that the injection of exo-NK cells significantly increased anti-tumor efficacy, survival time in experiments, and in silico models ^[74,76]. We assume that (ⅰ) NK cells are supplied from the TME at a rate $\lambda_{NK}$, (ⅱ) NK cell activities are suppressed by TGF-$\beta$ ^[125] with an inhibition strength $\theta$ and scaling parameters in a Hill-type function ($k_5, k_6$), (ⅲ) endogenous and exogenous NK cells go through the death processes at rates $\mu_K, \mu_{K'}$, respectively, and (ⅳ) exogenous NK cells are injected at a time-dependent rate $u_K(t)$. Based on these assumptions, the governing equations of densities of endogenous and exogenous NK cells, respectively, are given by

Non-dimensionalization of the dimensional version of governing equations (2.1)–(2.15) is shown in Appendix A. Matlab (Mathworks) was used for computation of the mathematical model and corresponding optimal control problems. See Tables 1 and 2 for parameter values of the mathematical model.

3.   Results

3.1.   Local dynamics of N1 and N2 TANs

We first investigate the fundamental dynamics of N2-N1 systems in Eqs (2.1) and (2.2) in response to two key constant signals, NK cells ($K_{end}$) and TGF-$\beta$ ($G$). In order to analyze the system structure at steady states, we can solve $N_2, N_1$ as a function of the NK cells ($K_{end}$) and TGF-$\beta\; (G)$, respectively.

We observe qualitatively different structures of the bifurcation curves of $N_2, N_1$ in response to the fixed NK signal ($K_{end}$) in the presence of low ($G = 0.1$; Figure 3A), intermediate ($G = 0.5$, Figure 3B), and high ($G = 0.9$; Figure 3C) TGF-$\beta$ levels. The equilibrium point at the lower and upper branches in Figure 3A are stable, while the steady state in the middle branch is unstable. See the details in Supporting Information S3 File. This naturally creates a bi-stability window where either N1 TANs or N2 TANs can dominate the tumor microenvironment depending on the initial configuration of the N1/N2 populations. For instance, we have a bi-stability window $W_K = [0.74, 1.35]$ in the NK cell spectrum when $G = 0.5$ in Figure 3B. This bi-stability window is shifted to the left (lower sector of NK cell signal) in response to low TGF-$\beta$ signal (Figure 3A) while the bistable region completely disappears in response to the high TGF-$\beta$ level (Figure 3C). Figure 3D–F illustrate the bifurcation curves of the densities of $N_2, N_1$ TANs in response to the fixed TGF-$\beta$ signal ($G$) in the presence of the low ($K_{end} = 0.1$; Figure 3D), intermediate ($K_{end} = 1.0$; Figure 3E), and high ($K_{end} = 2.0$; Figure 3F) NK cell activities. We observe the opposite effect on the up- or down-regulation of N1 (and N2) TANs for various TGF-$\beta$ levels. A bi-stability window emerges when NK cell activities are low (Figure 3D) but moves to the middle range of the TGF-$\beta$ level ($W_G = [0.32, 0.53]$, Figure 3E) as $K_{end}$ increases ($K_{end} = 0.1 \to 1$) and is eliminated in response to the high level of NK cell activities (Figure 3F). In both cases, the onset and size of the bi-stability windows ($W_K, W_G$) depend on the input signal ($K_{end}, G$) and parameter values. By taking the critical thresholds, $th_{N_1}( = 1.3)$ and $th_{N_2}( = 1.8)$, in the structure of the bifurcation curves in response to both NK cell signal and TGF-$\beta$ in Figure 3A–F, we can define the anti-tumorigenic ($\mathbb{P}_a$) and tumorigenic ($\mathbb{P}_t$) phases as follows:

When the TGF-$\beta$ signal is low ($G = 0.1$), the low NK cell activities lead to the co-existence of anti-tumorigenic ($\mathbb{P}_a$) and tumorigenic ($\mathbb{P}_t$) states, but the increased NK cell signals can induce $\mathbb{P}_a$ state (Figure 3A). For the intermediate level of TGF-$\beta$ ($G = 0.5$, Figure 3B), various stimuli from NK cells induce the nonlinear immune responses. For instance, a low level of NK cells (red asterisk in Figure 3B) induces $\mathbb{P}_t$-mode in given TME. On the other hand, intermediate levels of NK signals (blue asterisk in Figure 3B) can generate the bi-stable state ($\mathbb{P}_a$ or $\mathbb{P}_a$), resulting in either N1- or N2-dominant TME based on trajectories of solutions. As $K_{end}$ is increased further, passing the right knee ($K_{end} = 1.35$ of the N1 bifurcation curve (blue curve), the high level of NK cell activities enhance the immune balance toward the $\mathbb{P}_a$-state where tumor growth can be suppressed. Conversely, as the NK cell signal is slowly decreased from a high level (blue asterisk in Figure 3B), the system tends to maintain the $\mathbb{P}_a$-mode along the upper branch until it jumps down to the lower branch at the left knee point ($K_{end} = 0.74$), returning to the opposite mode ($\mathbb{P}_t$). Thus, the direction (increase or decrease) of the immune activities of NK cells can play a role in forming the tumorigenic or anti-tumorigenic mode of TANs in a given TME. When the TGF-$\beta$ signal is high ($G = 0.9$; Figure 3C), the immune system adapts to the N2-favorable mode $\mathbb{P}_t$ with a great potential of fast growing tumor regardless of the strength of the NK activities. We observe a different effect of TGF-$\beta$ on determining the TAN-induced characteristics of TME in Figure 3D–F. For the intermediate level of NK cell activities ($K_{end} = 1$; Figure 3E), the TAN system converges toward the $\mathbb{P}_t$ (or $\mathbb{P}_a$) in response to high (or low) TGF-$\beta$ stimuli while the intermediate TGF-$\beta$ level can lead to a mixed state ($\mathbb{P}_t$ or $\mathbb{P}_a$) with the bi-stability window $W_G$. However, when NK activities are suppressed ($K_{end} = 0.1$; Figure 3D), the system can adapt only to the bistable mode or $\mathbb{P}_t$. In this case, the TAN system does not return to the suppression mode ($\mathbb{P}_a$) even when the TGF-$\beta$ signal is decreased to zero from a high level with TGF-$\beta$ inhibitor treatment. On the other hand, when NK cell treatment is enhanced ($K_{end} = 1\to 2$; Figure 3E $\to$ Figure 3F), the bi-stability window $W_G$ disappears and the system gets the higher transition point from $\mathbb{P}_a$ to $\mathbb{P}_t$ as $G$ increases, resulting in the decreased N2-to-N1 ratio ^[85] and increased anti-tumorigenic potential. These are consistent with experimental observations ^{[10,11,12,13,14]}. Figure 3G shows the geometric illustration of the single mode ($\mathbb{P}_a$ (blue panel), $\mathbb{P}_t$ (pink panel)), and bi-stable mode ($\mathbb{P}_a$ or $\mathbb{P}_t$) corresponding to ($K_{end} =$ 0.1, 1.0, 2.0 in the $N_1-N2-K_{end}$ cube. The curves represent stable steady states ($N_1^s, N_2^s, K_{end}^s$) in the $N_1-N_2-K_{end}$ cube with the color coding (blue $\in \mathbb{P}_a$; red $\in \mathbb{P}_t$; green $\in \mathbb{P}_a \cup \mathbb{P}_t$). This implies that fluctuating NK cell activities can lead to active phenotypic switches of TANs between $\mathbb{P}_a$ and $\mathbb{P}_t$ with a nonlinear perturbation due to bi-stability in the middle. In a similar fashion, Figure 3H illustrates the structure of the tumorigenic or anti-tumorigenic regions and stable branches of equilibrium in the $N_1-N_2-G$ cube. Thus, changes in the concentration of TGF-$\beta$ can lead to different routes among different mono- or bistable-states of TANs. The phenotypic representation in response to various NK cell ($K_{end}$) and TGF-$\beta$ ($G$) stimuli is shown in Figure 3I. For a low, intermediate, and high level of NK cell signals, the system can be exposed to various environments ($\mathbb{P}_a$, $\mathbb{P}_t$, or mixed) as $G$ is increased. For instance, the system transits from $\mathbb{P}_a$ to a mixed state, and to $\mathbb{P}_t$ for the intermediate $K_{end}$ level. In a similar fashion, the system goes through $\mathbb{P}_t \to W_G \to \mathbb{P}_a$ transition for intermediate $G$ as $K_{end}$ increases. Overall, the system adapts to the bi-stable system for lower $K_{end}, G$ values and mono-stable system as either $K_{end}$ or $G$ is high. One can also observe that the relative balance between TGF-$\beta$ and NK cells can determine whether the TME is pro- or anti-tumorigenic. This also suggests that the combined therapy of TGF-$\beta$ inhibitor and exo-NK cells can effectively suppress the tumor growth.

3.2.   Local dynamics of intracellular signaling

In Figure 4A,B, we show the sensitivity of $K_{end}$-dependent Bcl-2 ($B$) and BAX ($X$) to changes in the inhibition strength parameter $\gamma$ of STAT1 by STAT3 ($\gamma = 0.5, 1, 1.5, 5$), respectively. As $\gamma$ is increased, both the bifurcation curve $B(K_{end})$ and bi-stability window ($D_K$) shift to the right, leading to a decrease in the size of the window ($|D_K|$). Here, we define the apoptotic ($\mathbb{T}_{a}$) and anti-apoptotic ($\mathbb{T}_{t}$) phases based on the threshold values ($th_B, th_X$):

As $\gamma$ is decreased, the probability of the phenotypic switching to the apoptosis phase ($\mathbb{T}_{a}$) is increased. For example, the Bcl-2 level is already in the $\mathbb{T}_{a}$-phase ($B < th_B, X > th_X$) when $K_{end} = 1$ in the case of the lower $\gamma$ ($\gamma = 0.5$; green solid curves in Figure 4A,B) while it is still in the anti-apoptosis phase ($\mathbb{T}_t$; $B > th_B, X < th_X$) in the base case ($\gamma = 1.5$; solid red curves in Figure 4A,B). On the other hand, the system adapts to a tumor-favoring situation even in the presence of the high level of NK cell activities as $\gamma$ is increased. For example, the system will be in the $\mathbb{T}_{t}$-phase for a high NK signal ($K_{end} = 3$) when $\gamma$ is high ($\gamma = 5$; black dotted curves in Figure 4A,B) while it should be in the $\mathbb{T}_{a}$-phase in the base case ($\gamma = 1.5$). Figure 4C shows the phenotypic transition from a $\mathbb{T}_{a}$-phase to a mixed phase ($\mathbb{T}_{a}+\mathbb{T}_{t}$), and to a $\mathbb{T}_t$-phase under various NK cells conditions as $\gamma$ increases. While the system adapts to the $\mathbb{T}_{a}$-dominant TME regardless of strength of NK cell activities for small $\gamma$ ($0 \le \gamma < \gamma^{MIN}$), it forms the $\mathbb{T}_{t}$-dominant TME regardless of NK cell conditions for larger $\gamma$ ($\gamma > \gamma^{MAX}$). For the intermediate levels of $\gamma$ ($\gamma^{MIN} < \gamma < \gamma^{MAX}$), the increased NK cell activity leads to the mode switches toward the apoptotic death of cancer cells: $\mathbb{T}_{t} \to (\mathbb{T}_{a} \cup \mathbb{T}_{t}) \to \mathbb{T}_{a}$. In Figure 4D,E, we show the sensitivity of $K_{end}$-dependent Bcl-2 ($B$) and BAX ($X$) to changes in the inhibition strength parameter of STAT3 by STAT1 ($\delta = 0.5, 1, 1.5, 6$), respectively. As $\delta$ is decreased, both the bifurcation curves $(B(K_{end}), X(K_{end}))$ and the bi-stability window ($D_K$) shifts to the right with a decrease in the size of the window ($|D_K|$).

An increase in $\delta$ leads to higher probability of the phenotypic switching to the apoptosis phase ($\mathbb{T}_{a}$). For example, the level of Bcl-2 is already in the $\mathbb{T}_{a}$-phase ($B < th_B, X > th_X$) when $K_{end} = 1$ in the case of the higher $\delta$ ($\delta = 6$; green solid curves in Figure 4D,E) while it is still in the anti-apoptosis phase ($\mathbb{T}_t$; $B > th_B, X < th_X$) in the base case ($\delta = 1$; solid red curves in Figure 4D,E). On the other hand, the system favors a tumor-promoting TME even in the presence of high NK cell activities when $\delta$ is small. For example, the system would be in the $\mathbb{T}_{t}$-phase for a high NK cells signal ($K_{end} = 3$) when $\delta$ is small ($\delta = 0.5$; black solid curves in Figure 4D,E) while it should be in the $\mathbb{T}_{a}$-phase in the base case ($\delta = 1$). Figure 4F shows the cell fate distribution of cancer cells for various NK cell levels and $\delta$'s. For a fixed NK cell activity, an increase in $\delta$ leads to a phenotypic transition from a $\mathbb{T}_{t}$-phase to a mixed phase ($\mathbb{T}_{a}+\mathbb{T}_{t}$), and to a $\mathbb{T}_a$-phase. While the system favors the $\mathbb{T}_{t}$-dominant TME regardless of strength of NK cell activities for small $\delta$ ($0 \le \delta < \delta^{MIN}$), it forms the $\mathbb{T}_{a}$-dominant TME regardless of NK cell conditions for larger $\delta$ ($\delta > \delta^{MAX}$). For the intermediate levels of $\delta$ ($\delta^{MIN} < \delta < \delta^{MAX}$), the increased NK cell activity leads to two types of mode switches toward the apoptotic death of cancer cells: (ⅰ) $\mathbb{T}_{t} \to (\mathbb{T}_{a} \cup \mathbb{T}_{t}) \to \mathbb{T}_{a}$ ($\delta^{MIN} < \delta \le \delta^\dagger$) (ⅱ) $(\mathbb{T}_{a} \cup \mathbb{T}_{t}) \to \mathbb{T}_{a}$ ($\delta^{\dagger} < \delta < \delta^{MAX}$).

3.3.   Sensitivity analysis

In the model developed in this paper, there are several parameters for which no experimental data are available, and these parameters may affect the simulation results. We selected all parameters in the model ($\lambda_6, \lambda_G, k_1, k_2, \alpha, \mu_{N2}, \lambda_{S1}, k_3, k_4, \beta, \mu_{N1}, \lambda_{S2}, a_1, a_2, \gamma, \mu_{S1}, \lambda_1, K_1, \lambda_{S3}, a_3, a_4, \delta, \mu_{S3}, \lambda_2, a_5$, $a_6, \omega, \lambda_3, \mu_B, \lambda_4, a_7, a_8, \zeta, \mu_X, \lambda_K, \mu_S, \lambda_C, \mu_G, \lambda_{NK}, k_5, k_6, \theta, \mu_K, D, \alpha_{RT}, \beta_{RT}, \epsilon_{N2}, \epsilon_{N1}, \lambda_{RK}, \lambda_{RG}, r, \delta_1$, $\delta_2, \delta_3$) for sensitivity analysis. We investigated the sensitivity of the major variables (N1 TANs, N2 TANs, STAT1, STAT3, Bcl-2, BAX, IFN-$\gamma$, NK cells, TGF-$\beta$, and Cancer cells) to these parameters at various time points. Specifically, we computed partial rank correlation coefficients (PRCCs) for $N_1, N_2, S_1, S_3, B, X, S, K_{end}, G,$ and $C$ at times $t = 30$ days in the model. We established a range for each of these parameters and divided each range into 10,000 intervals of uniform length. For each of these 54 parameters of interest, a PRCC value was calculated. A PRCC value is a real number between -1 and 1, with the sign indicating whether an increase in the parameter value either decreases (-) or increases (+) each variable at a given time. The sensitivity analysis described below was carried out using the method (General Latin Hypercube Sampling (LHS) scheme and Partial Rank Correlation Coefficient (PRCC)) ^[127].

Figure 5 shows the PRCC values of all variables at $t = 30 \; day$. Most of eleven parameters ($\lambda_6, \lambda_G, k_1, k_2, \alpha, \mu_{N2}, \lambda_{S1}, k_3, k_4, \beta, \mu_{N1}$) in the N1/N2 TANs system demonstrated sensitivity exclusively to either N1 or N2 TANs. Specifically, the N2 TANs activities are positively correlated with $\lambda_G, k_1, k_2$, $\beta, \mu_{N1}$, indicating that increases in these parameters promote the N2 phenotype. Conversely, the N2 TANs activities are negatively correlated to $\alpha, \mu_{N2}, \lambda_{S1}, k_3, k_4$, highlighting their suppressive effect on the N2 phenotype. This analysis underscores the differential impact of each parameter on the balance between N2 and N1 TANs. On the other hand, 23 parameters involved in STAT signaling pathway show specific correlations with the players (STAT1 ($S_1$), STAT3 ($S_3$), BAX ($B$), and BAX ($X$)) in the signaling network. For instance, the parameters $a_1, a_2, \mu_{S3}, \mu_B, a_7,$ and $a_8$ exhibited positive correlations with the BAX level, indicating that an increase in these parameters enhances the apoptotic phenotype. Conversely, parameters $\mu_{S1}, a_3, a_4, \lambda_3, \zeta,$ and $\mu_X$ showed negative correlations, suggesting that the larger value of these parameters suppresses the apoptosis of cancer cells. This analysis highlights the key factors influencing apoptotic regulation within the STAT signaling pathway. The rest of parameters are involved in the interactions between cancer cells, radiation therapy, and the immune system. This analysis evaluates how variations in these parameters influence key system dynamics, highlighting the relative impact of each parameter on tumor progression, immune response, and the efficacy of radiation therapy. The N1 TAN activities are positively correlated with $\lambda_K$ and $\lambda_{RK}$ but negatively correlated with $\mu_K$ and $\epsilon_{NK}$. The N2 TAN activities, located on the opposite side in the TAN spectrum, show the opposite effect from changes in these parameters. On the other hand, the IFN-$\gamma$ level is positively correlated with $\lambda_K, \lambda_{RK}, D,$ and $\alpha_{RT}$ but negatively correlated with $\mu_K$. In a similar fashion, the STAT1 level and BAX expression are positively correlated with $\lambda_K, \lambda_{RK}, D,$ and $\alpha_{RT}$ but negatively correlated with $\mu_K$ due to the feed-forward structure of the intracellular signaling network. However, the antagonist STAT3 and Bcl2 in the direct downstream shows the opposite effect: a negative correlation with $\lambda_K, \lambda_{RK}, D,$ and $\alpha_{RT}$ and a positive correlation with $\mu_K$ due to mutual antagonism between STAT1 and STAT3 upstream of the signaling. NK cell activities are negatively correlated with $\mu_K$ but positively correlated with radiation-induction parameters $D, \alpha_{RT}, \beta_{RT},$ and $\lambda_{RK}$. TGF-$\beta$ is positively correlated with radiation dose packages $D, \alpha_{RT}, \beta_{RT},$ and $\lambda_{RG}$, indicating strong impact from radiation, but do not show particular negative correlations with other parameters.

To ensure a thorough evaluation of model robustness, we have conducted sensitivity analysis across a biologically relevant range of parameters. The specific parameter ranges used in sensitivity analysis, including minimum, maximum, and baseline values, are summarized in Table 3.

3.4.   Therapeutic strategies

We now investigate the therapeutic effect of TGF-$\beta$ inhibitor, IFN-$\gamma$, and exo-NK cells on tumor growth in this section.

These anti-tumor agents were injected once every 5 days by prescribing periodic functions $u_L(t), u_S(t),$ and $u_K(t)$ in Eqs (2.12), (2.13), and (2.15), respectively, as follows (Figure 6):

Here, we set $h_L = 1, L_s = 0.5, \tau_L ( = t_{j+1}-t_j) = 5$ (pink boxes in Figure 6A1), $h_S = 1, S_s = 0.4, \tau_S ( = t_{j+1}-t_j) = 5$ (blue boxes in Figure 6B1), $h_K = 1, K_s = 0.08,$ and $\tau_K ( = t_{j+1}-t_j) = 5$ (green boxes in Figure 6C1). The periodic injection of these therapeutic agents lead to fluctuation in these concentrations within TME (Figure 6A1–C1). While the TGF-$\beta$ level (blue curve in Figure 6A1) is dramatically reduced by its inhibitor, it is not strong enough to transit the N2 TANs to N1 TANs due to the short duration of injection (Figure 6A3), leading to inefficient suppression of tumor growth (Figure 6A4). However, the elevated IFN-$\gamma$ level (red curve in Figure 6B1) forms a N1-dominant TME (high N1 TAN and low N2 TAN in Figure 6B3), resulting in efficient cancer cell killing (6B4). Note that while this up-regulated IFN-$\gamma$ slowly decreases the Bcl2 level, instant elevation of IFN-$\gamma$ is not strong enough to upregulate the BAX level during this process (Figure 6B2). When NK cells are injected into the system, direct cancer cell killing from increased numbers of NK cells (Figure 6C1), indirect cell killing via intracellular mechanisms (Figure 6C2), and N1-mediated killing from increased IFN-$\gamma$ (Figure 6C3) can effectively reduce the tumor size (Figure 6C4).

From now on, we investigate anti-tumor efficacy by monitoring relative TME states than periodic injections of anti-cancer drugs. We first investigate how controlled IFN-$\gamma$ injection can affect tumor growth by maintaining either an N1-dominant TME or apoptotic status in cancer cells. Figure 7A–C show time courses of the density of N1 and N2 TANs, key molecules in apoptosis pathway (Bcl2 and BAX), and tumor population, respectively, in response to IFN-$\gamma$ injection (blue boxes) at $t = 0, 2, 20.8,$ and 27 days.

While the injection duration was fixed, the injection time was determined so that N1-dominant TME is maintained ($N_1 > th_{N_1}$; blue curve in Figure 7A), leading to inhibition of tumor growth (Figure 7C; reduction of tumor size by 22%). In this case, the system does not adapt to apoptosis status (i.e., still in the state of down-regulation of BAX and up-regulation of Bcl-2; Figure 7B). Figure 7D–F show the time courses of TAN density (N2 and N1 TANs), concentrations of Bcl-2 and BAX, and tumor population, respectively, in the presence of apoptosis-controlled IFN-$\gamma$ therapy. To prevent a rapid drop in BAX levels and maintain the apoptosis status (up-regulation of BAX and down-regulation of Bcl-2; Figure 7E) via STAT signaling, frequent IFN-$\gamma$ injections are needed in the beginning. The next injection time was determined so that the Bcl-2 level stays below the threshold, maintaining the $\mathbb{T}_{a}$-mode. In this strategy, frequent injection of IFN-$\gamma$ at the early stage also induces the N1-dominant TME (Figure 7D). Thus, the combined effect ($\mathbb{T}_{a}$+$\mathbb{P}_a$) leads to a significant reduction in the tumor population (63% reduction compared to the control; Figure 7F).

In Figure 8, we consider the combination treatment of IFN-$\gamma$ and TGF-$\beta$ inhibitors. Galunisertib (LY2157299), a type of TGF-$\beta$ inhibitor, is typically administered for 14 days, followed by a 14-day rest period in clinics ^[109,110]. Following the protocol, TGF-$\beta$ inhibitor was injected for 14 days, beginning at day $T_\kappa$, in addition to injection of IFN-$\gamma$, four times at $t = 0, 2, 20.8, 27$ days as in Figure 7A–C. Figure 8A,B illustrate the injection rate of IFN-$\gamma$ (blue box; $u_S = 0.25$) and TGF-$\beta$ inhibitor (red box; $u_L = 0.05$) when $T_\kappa = 0$ and $T_\kappa = 12$, respectively.

Injection of the TGF-$\beta$ inhibitor in the early stage ($T_\kappa = 0$; red solid curve in Figure 8C) outperforms the case of injection in later time ($T_\kappa = 12$; diamond solid curve in Figure 8C) in suppression of tumor growth due to extensive N1-mediated cancer cell killing at later times and extended duration of the anti-tumorigenic mode (Figure 8D). However, the system maintains the anti-apoptosis phase in both cases (Figure 8E). Figure 8F shows the nonlinear anti-tumor efficacy as a function of $T_\kappa$ ($0 \le T_\kappa \le 16$). The maximum anti-tumor efficacy is obtained at an intermediate injection time of TGF-$\beta$ inhibitor ($T_\kappa = 12$; 20% reduction in the tumor population relative to the worst case ($T_\kappa = 0$)).

The third strategy depicted in Figure 9 involves the use of NK cells. Figure 9A–C show the time courses of TAN density (N1 and N2 TANs), concentrations of Bcl-2 and BAX, and tumor population, respectively, when exo-NK cells were injected to maintain the N1-dominant TME. Two injections at $t = 0, 2$ days (green boxes in Figure 9A) were enough to upregulate N1 TANs (blue curve in Figure 9A) and down-regulate N2 TANs (red curve in Figure 9A) while the system is under $\mathbb{T}_{t}$-mode ($B > th_B, X < th_X$) for most of the time (Figure 9B). On the other hand, in the case of apoptosis-controlled NK therapy, additional injections of exo-NK cells (green boxes in Figure 9D) were needed in order to maintain the apoptosis status ($B < th_B, X > th_X$; Figure 9E) via STAT signaling while maintaining the $\mathbb{P}_{a}$-mode ($N_1 > th_{N_1}, N_2 < th_{N_2}$; Figure 9D). The tumor size in the apoptosis-controlled strategy (Figure 9F; 98% reduction compared to control) is better controlled than the case of TAN-controlled strategy (Figure 9C; 66% reduction compared to control) despite higher administration costs.

NK cells were shown to play a significant role in modulating the tumor microenvironment, i.e., either enhancing or decreasing anti-tumor efficacy in a combination therapy depending on the nature of interference of anti-tumor and immunosuppressive effects ^[74]. We consider immune-based combination therapy (exo-NK + TGF-$\beta$ inhibitor) in Figure 10 to test the role of NK cells in regulating immune-tumor dynamics and cancer cell killing in the presence of TGF-$\beta$ inhibitor treatment. Here, a fixed amount of TGF-$\beta$ inhibitor was continuously injected over the time interval $[T_{\kappa}, T_{\kappa}+\tau]$ with various infusion times $T_{\kappa}$ ($0 \le T_\kappa \le 16$, $\tau > 0$; pink boxes in Figure 10A,B) while exogenous NK cells were injected twice in the early stage as in Figure 9A (green boxes in Figure 10A,B). Figure 10C–E show time courses of cancer cell densities and corresponding trajectories of solutions ($N_2(t), N_1(t)$) and ($B(t), X(t)$) in the $N_2-N_1$ and $B-X$ plane, respectively, in control and NK+ cases ($T_\kappa = 0$ (solid, worst case), and 12 (diamond, best) in Figure 10C). While two injections of exogenous NK cells in these two cases were enough to switch N2 TANs to N1 TANs in the TME (Figure 10D), these kept cancer cells in the anti-apoptosis of the cancer cells (Figure 10E). Figure 10F summarizes the normalized tumor population at final time in response to the combination therapy with various infusion times ($T_\kappa$; $0 \le T_\kappa \le 16$). The maximum anti-cancer efficacy is obtained when $T_\kappa = 12$ while the efficacy is decreased as $T_\kappa$ is further decreased.

We finally consider a combination of two drugs (IFN-$\gamma$ and TGF-$\beta$ inhibitor) and exo-NK cells in Figure 11. While IFN-$\gamma$ (blue box) and NK cells (green box) are injected 4 times and 2 times as before, various infusion times $T_{\kappa}$ of TGF-$\beta$ inhibitor (pick box; $0 \le T_\kappa \le 16$) were selected for testing (Figure 11A,B). Figure 11C shows the time courses of tumor populations when $T_{\kappa} =$16, 2 in response to triple combination treatment with the late (Figure 11A) and early (Figure 11B) injection of TGF-$\beta$ inhibitor. In these cases, the combination treatment leads to the phenotypic transition from N2 TANs to N1 TANs (Figure 11D). Unlike sustained anti-apoptotic status in the NK+TGF-$\beta$ inhibitor combination in Figure 10, these triple treatments also activate apoptosis signaling pathways within the cancer cells (Figure 11E). Figure 11F shows the tumor population at final time in control (red solid) and combination therapy with various $T_\kappa$ values ($0 \le T_\kappa \le 16$). In Figure 11G, we summarize the normalized tumor population at the final time points in response to three combination treatments (IFN-$\gamma$+TGF-$\beta$ inhibitor; NK+TGF-$\beta$ inhibitor; IFN-$\gamma$+NK+TGF-$\beta$ inhibitor).

The normalized tumor populations corresponding to different injection combinations of IFN-$\gamma$, TGF-$\beta$ inhibitor, and NK cells are shown in the circular form in Figure 12A–C where the center represents the 0 population and the far edge along the radial direction indicates the 1 value. Thus, the length of the radial line is the normalized tumor population in $[0, 1]$ corresponding to each injection scheme. The different schemes were categorized by three groups based on the initial agent of IFN-$\gamma$ ($S$ in (A)), NK cells ($K$ in (B)), and TGF-$\beta$ inhibitor ($L$ in (C)). Figure 12D–F show the temporal profiles of the $\mathbb{P}_a$- and $\mathbb{T}_a$-modes for various combinations of the combination therapies in Figure 12A–C. Simulation results suggest that combination 'LLSSKK' is the least effective in killing tumor cells while the 'KSSKLL' scheme is the most effective in increasing anti-tumor efficacy (Figure 12B,C,J). In the best scenario 'KSSKLL', we observe a wide range of the apoptotic status (blue region in Figure 12I) in response to the initial injection of NK cells and IFN-$\gamma$ at earlier time points (Figure 12G) and the sustained anti-tumorigenic mode (blue region in Figure 12H) in response to TGF-$\beta$ inhibitor treatment at later times (Figure 12G). On the contrary, in the worst scheme ('LLSSKK'), earlier two infusions of TGF-$\beta$ inhibitor (red bars, Figure 12K) are not as effective in switching TANs to the $\mathbb{P}_a$-mode (blue region in Figure 12L) and the system successfully transits to the $\mathbb{P}_a$-mode only after IFN-$\gamma$ injection after $t = 10$ days. Two injections of IFN-$\gamma$ in the middle temporal spot and additional injection of exogenous NK cells at later times ($t > 20\; days$, Figure 12K) is not strong enough to induce the apoptotic status ($\mathbb{T}_a$) later times (Figure 12M). These two factors, delay in the systemic transition to $\mathbb{P}_a$-mode and lack of $\mathbb{T}_a$-mode, result in relatively weak anti-tumor efficacy (blue curve, Figure 12J). In general, early injections of IFN-$\gamma$ tend to form the wide range of $\mathbb{P}_a$ status but the narrow or none of $\mathbb{T}_a$-mode (Figure 12D). These opposite effects lead to the intermediate level of cancer cell killing (left panel in Figure 12N). On the other hand, early injections of TGF-$\beta$ inhibitor induce the narrower ranges of $\mathbb{P}_a$ in middle and small range of $\mathbb{T}_a$-modes at later times (Figure 12F), leading to lower anti-tumor efficacy (right panel in Figure 12N). Finally, initial infusion of NK cells induce well-shaped ranges of $\mathbb{P}_a$- and mixed ($\mathbb{P}_a \cup \mathbb{T}_a$)-status, leading to effective, combined cancer cell killing (middle panel in Figure 12N).

3.5.   Radiotherapy

In this section, we investigate how the complex interactions between IR and immune cells in the TME (Figure 1) regulate tumor growth and develop a new anti-cancer strategy by combining the TGF-$\beta$ inhibitor to RT therapy. Here, we propose a mathematical model of tumor growth via tumor-NK-TAN interactions and cancer cell survival dynamics in the presence of radiotherapy (Figure 1). The linear-quadratic (LQ) model, one of the standard tools in radiation biology, has been used to calculate cell survival in response to delivered radiation dose ^[128,129]. In this framework, the survival probability $[RT]$ of cells is calculated as follow ^[130]:

where the biologically effective dose (BED), $BED = \alpha_{RT}D+\beta_{RT}D^2$, is obtained from cell-specific radiosensitivity parameters ($\alpha_{RT}$ and $\beta_{RT}$) and the physical dose delivered per fraction ($D$). With the known BED for healthy cells and cancer cells, the killing rate of the cells by radiation is given by

where $\mu_0$ is a basic killing rate due to radiation. We use the well-known relation $\alpha_{RT}/\beta_{RT} = 10Gy$ ^[81,93,131] in this study. Based on experimental observation, we also assume that the effect of radiation on N1 and N2 TANs is lower than one for cancer cells. Thus, we set $\alpha_{RT} = 0.03 \; Gy^{-1}, \beta_{RT} = 0.003 \; Gy^{-2}$.

In order to take into account the effect of the radiotherapy, we now use the modified version of Eqs (2.1), (2.2), (2.6), (2.14), and (2.11), for the densities of N1/N2 TANs, cancer cells, and NK cells, and concentration of TGF-$\beta$ as follows:

where $\varepsilon_{N2}, \varepsilon_{N1}$ are a scaling factor due to the low sensitivity of IR damage ($\varepsilon_{N2}, \varepsilon_{N1} \ll 1$), $\lambda_{RK}, \lambda_{RG}$ are production rate of NK cells and TGF-$\beta$ by IR, respectively ^[81]. We consider a treatment strategy using IR and TGF-$\beta$ inhibitor through the Eqs (3.10)–(3.12). We set $\varepsilon_{N2} = \varepsilon_{N1} = 0.5$. We also set $\lambda_{RK} = 0.05 \; g/cm^2$ and $\lambda_{RG} = 0.01 \; g/cm^2$.

In Figure 13, we investigate the effect of radiation on regulation of changes in the N1/N2 TAN composition in TME and overall anti-tumor efficacy. In this model, we set $\alpha_{RT} = 0.03$ and $\beta_{RT} = 0.003$ following ^[81,93,131]. Radiation therapy can cause neutropenia (low neutrophil count) as it damages the rapidly dividing cells in the bone marrow, where white blood cells are produced. The severity and duration of radiation-induced neutropenia depend on the dose and location of the radiation therapy, as well as other individual factors such as age, overall health, and the presence of other medical conditions ^[132,133]. In our model, the radiation-mediated cell killing from different doses of radiation (Figure 13A, $D = 0, 1, 3, 5$) lowers the overall population and the landscapes of the levels of N1 and N2 TANs (Figure 13B,C) relative to control ($D = 0$), illustrating consistent decrease in neutrophils in TME as shown in experiments. A decrease in TAN populations may interfere with strong anti-cancer activities by N1 TANs in the presence of NK cells or decrease the tumor-promoting effect by N2 TANs. As $D$ increases, the bi-stable region disappears (Figure 13B,C). Radiotherapy is typically administered to cancer patients five times per week, with treatment schedules lasting from 3 to 9 weeks ^{[134,135,136]}. In our model, we simulate a total of 4 weeks of radiation therapy, with treatment administered five times per week. Figure 13D–E show the normalized tumor population and corresponding time courses for various radiation dose levels ($D = 1$ (orange), 2 (yellow), and 3 (purple) relative to control ($D = 0$; blue curves). Radiation with strong dose is quite effective in killing cancer cells ($D\ge2$) but it can induce neutropenia ^[132,133]. During this radiation period, the system also stays in the tumorigenic mode even with low levels of radiation (Figure 13F).

The optimal control theory ^[137] was utilized to find optimal infusion strategies of anti-cancer drugs that minimize the tumor population and minimize administrative costs in mathematical oncology. Optimal control theory has been adapted to develop effective administration schedules of anti-cancer drugs in various types of cancer progression ^{[138,139,140,141,142]}, including glioblastoma ^{[76,143,144,145,146,147]} and lung cancers ^[107,148]. In particular, an optimal control theory was used to block aggressive glioma cell infiltration by maintaining miR-451 levels ^[147] and regulating the cell cycle ^[146], or to optimize the cancer cell killing via controlling a nonlinear immune-tumor interaction such as NK cells ^[76] and traditional drugs such as bortezomib with OVs via the necroptosis pathways ^[143]. Here, we use an optimal control theory to reduce the total cancer population $\int C(t) dt$ and minimize the side effects and administrative costs $u_{RT}$ associated with radiation therapy by obtaining radiation schedule of $u_D$. The control $u_D$ represents dosage of ionizing radiation. This control is also assumed to be bounded. That is, the control set is defined as

where $u_{D}^{max}$ is the maximum dosage of ionizing radiation, which is administered in a given treatment period, and $t_s$ and $t_f$ are the start and end time of treatment with optimal control applied, respectively. The cost (side effects + administrative costs) associated with radiation therapy is calculated by

where $u_0$ is a proportional constant. An objective function is to find the optimal dose strategy (temporal profiles and dose) of radiation over time for the minimal tumor population and minimal costs, leading to an objective function as follows:

where $C_1$ and $C_2$ are the weight constant for each control. Both linear and quadratic terms in the integrand are used in costs associated with radiation in order to regularize the radiation dose in the system. Linear controls are more difficult to analyze compared to quadratic controls despite their higher degree of biological relevance. Numerical solutions of the control problems (3.14) and (3.15) was obtained using the forward-backward sweep method with shooting methods ^[137].

Since radiation therapy is performed five times a week ^{[134,135,136]} in a clinical setting, we applied the optimal control to adjust the radiation dose to achieve the goal of minimizing the tumor size within a time cycle (1 week). Figure 14A shows the time courses of radiation doses without (dashed blue) and with (solid pink) the application of optimal control theory, respectively. We can reduce the tumor size at the final time by optimal control (Figure 14C) with lower radiation dose (Figure 14B). Thus, optimally controlled schedule can lead to better clinical outcomes compared to regular radiation schedule with fixed dose and duration. The radiation therapy can kill tumor cells ($> 50\%$) relative to control (w/o treatment) as shown in experiments with the non small cell lung cancer cells implanted in C57Bl6 mice ^[149].

Radiotherapy (RT) is a common treatment option for cancer patients, and it is well known that it can stimulate a robust immune response, providing a strong basis for the combination of RT and immunotherapy (iRT) ^[150]. While radiation kills cancer cells and stimulates NK cell activities, it also stimulates TGF-$\beta$ secretion. Furthermore, the radiation dose is usually set low in order to minimize the biophysical and biochemical strain on the patient's body. Here, we consider the iRT combination therapy (RT + TGF-$\beta$ inhibitor) with low dose radiation and optimal control theory. TGF-$\beta$ is a protein that plays a role in cell growth and division and is also involved in cancer development ^[48,80]. Preclinical studies have shown that combining radiotherapy with TGF-$\beta$ inhibitors can enhance the anti-tumor efficacy of radiotherapy by increasing tumor cell killing and reducing the growth, spread of tumors, and angiogenesis ^[151,152]. TGF-$\beta$ inhibitors have also been shown to decrease the adverse effects of radiotherapy on normal tissues, such as reducing inflammation and fibrosis ^[153,154]. While iRT is not yet a standard treatment option for cancer, it shows a great potential of eradicating tumor cells. iRT was also suggested to promote abscopal effect, i.e., radiation-mediated anti-tumor immune response with the regression of non-irradiated tumors at distant sites ^[150]. To investigate the anti-tumor efficacy of the combination therapy, two injection schedules of the TGF-$\beta$ inhibitor were selected in addition to regular RT (Figure 14): $T_{\kappa} = 16$ (Figure 15A), $6$ (Figure 15B). Suppression of the TGF-$\beta$ level (red curves in Figure 15C) in response to the TGF-$\beta$ inhibitor induces the transition from the tumorigenic to anti-tumorigenic mode (Figure 15D) in both cases. Ionizing radiation in the combination therapy further upregulates NK cell activities (blue curves in Figure 15C), which increases anti-tumor efficacy. Thus, this combination therapy leads to effective cancer cell killing compared to control (RT$^-$TGF-$\beta$ inhibitor$^-$) (Figure 15D). Figure 15F shows normalized tumor populations for various $T_\kappa$'s ($T_\kappa = 1, \dots, 16$), implying the importance of scheduling of the TGF-$\beta$ inhibitor. For instance, the tumor size when $T_\kappa = 6$ is reduced by about 40% compared to the worst case ($T_\kappa = 16$) and NK cell activities are also higher when $T_\kappa = 6$.

3.6.   RT-induced neutropenia

In order to treat stage 3 or 4 cancer patients, a large amount of radiation is necessary ^{[155,156,157]} but this can cause a serious side effect, such as neutropenia. Neutropenia, characterized by a deficiency of neutrophils (comprising N1 TANs + N2 TANs) in the bloodstream, is one of frequent complications of extended-field radiotherapy. This condition often leads to treatment disruptions and elevates the susceptibility to infections or bleeding ^[158,159]. For patients undergoing simultaneous chemotherapy and radiotherapy targeting the mediastinum, the administration of CSF support is not a favorable option due to heightened risks of complications and mortality ^[160,161]. However, in scenarios where chemotherapy is not part of the treatment regimen, CSF support could be contemplated for patients undergoing RT alone if there is substantial and prolonged neutropenia, particularly when RT involves a considerable volume of active bone marrow. The spontaneous resolution of neutropenia following RT alone can be protracted, often taking around 3–4 weeks after exposure to an intermediate dose of total body radiotherapy ^[162]. Hence, the utilization of a medication known as neupogen (G-CSF) is warranted to address neutropenia. Here, we develop a new combined RT treatment strategy for a high-grade tumor where the patient is treated with a high radiation dose and neutropenia is minimized. We set $C(0) = 0.5$ to reflect the later stage tumor. In this framework, neutropenia is defined by

where $th_N = 3$. Figure 16A shows the temporal profiles of $N_i$ for various radiation dose ($D = 0, 1, 2, 3, 4, 5$). When the radiation dose is low ($D =$0 or 1 in Figure 16A), the system does not induce neutropenia during the treatment ($N_i(t) > 1, \forall t\in [0, 30]$) though the anti-tumor efficacy is also low ($D = 1$, red dashed curve in Figure 16B). A higher radiation dose ($D \ge 2$) is needed to eradicate tumor cells for stage 3 or 4 cancer patients, as shown in Figure 16B. However, these higher radiation dose strategies lead to neutropenia ($D \ge 2$; yellow boxes in Figure 16A). To reduce the duration of radiation-induced neutropenia, specific agents such as neupogen ^[163] have been used. For the following simulations, we set $D = 2$ for high dose RT. In Figure 16C–H, we investigate the effect of neupogen on controlling neutropenia during RT. We assume that (ⅰ) neupogen is injected during time interval $[T_\kappa, T_\kappa+\tau]$ at a fixed injection rate $u_N = 0.1733$, and (ⅱ) neupogen is administered for 14 days (i.e., $\tau = 14\; days$), similar to the TGF-$\beta$ inhibitor treatment (Figure 15), to control neutropenia following ^[164]. By modifying Eq (3.8) for N1 TAN and introducing a new variable, neupogen concentration ($N$), we have the following dynamics:

Here, $\lambda_N$ is the source of N1 TANs from bone marrow stem cell from neupogen, $\mu_N$ is the decay rate of neupogen, the indicator function $I_{[a, b]}$ gives 1 when $t \in [a, b]$, or zero otherwise.

Figure 16C–E show the temporal changes in the TAN (N2 (red), N1 (blue)) densities in control (RT), and combination (RT+neupogen) therapy with $T_\kappa = 0$ (D), 10 (E). The neutropenia index ($N_i$) stays above the threshold value for most of the time when neupogen is injected after 10 days ($T_\kappa = 10$, red curve in Figure 16F), leading to weak neutropenia. On the other hand, the index drops down below 1 with a significant duration in response to the early injection of neupogen ($T_\kappa = 0$, blue curve in Figure 16F), leading to cyclic neutropenia (yellow boxes in Figure 16F). Temporal discrete patterns of neutropenia in response to neupogen treatment of sequential injection time ($T_\kappa = 0, 1, \dots, 16$) are shown in Figure 16G. Earlier injection of neupogen tends to generate more frequent, extensive neutropenia, including wider ranges of neutropenia at later times. Administration of neupogen at later times induces narrow strips of neutropenia but do not form these later times. Figure 16H summarizes the cumulative duration of RT-induced neutropenia (CDRN) for various $T_\kappa$'s ($T_{\kappa} = 0, 1, \dots, 16$). Neupogen infusion at 0 day (blue) and 10 days (red) gives rise to the least and most favorable outcomes, respectively. In the best scenario ($T_{\kappa} = 10$, red box in Figure 16H), the cumulative duration of neutropenia was reduced by 90% compared to the worst case scenario ($T_{\kappa} = 0$, blue box in Figure 16H). Therefore, our results suggest that even with the use of high-dose radiation in stage 3–4 cancer patients, the proper administration of neupogen can mitigate neutropenia and effectively eliminate cancer cells.

3.7.   Spatial dynamics in the heterogenous tumor microenvironment

In this section, we extend the ODE model to the PDE model in order to capture the spatio-temporal evolution of cancer cells in the presence of NK cells and TANs in a heterogenous tumor microenvironment. The mass balance equation for the main variable $w({\mathbf x}, t)$ at space ${\mathbf x}$ and time ($t$) is given by

where $J_w$ is the flux, and $P_w$ is the net production rate of the variable $w$. The flux $J_w$ is due to the cell motility ($C, \; N_1, \; N_2, \; K_{end},$ and $K_{ex}$) or diffusion process of diffusible molecules ($G, S$), such as

where $D_w({\mathbf x)}$ is the space-dependent diffusion coefficient of $w$. By taking the reaction part of the original governing equations (2.6)–(2.13) for $P_w$ for each variable, we get the following governing equations in the heterogenous TME:

The list of parameters of diffusion coefficients are given in Table 4. Our computation domain is $\Omega = [0, 1]$ in a dimensionless version, and we impose Neumann (no-flux) boundary conditions for all relevant variables ($w = C, \; N_1, \; N_2, \; K_{end}, \; K_{ex}, G, S$):

We first investigate the role of spatial distribution of NK cells and N1/N2 TANs in the regulation of tumor growth on a spatial domain $\Omega$ with two subdomains, $\Omega_1$ and $\Omega_2$: $\Omega_1 = [0, 0.5]$, $\Omega_2 = [0.5, 1]$, $\Omega = \Omega_1 \cup \Omega_2$. We set the initial conditions as follows: (ⅰ) The cancer cells are positioned at the center (Figure 17A-1), $C(x, 0) = e^{-100*(x - 0.5)^2}$; (ⅱ) N1 TANs are immersed within the tumor on the left ($\Omega_1$) whereas N2 TANs already reside in the tumor region on the right ($\Omega_2$ (Figure 17B-1), $N_1(x, 0) = 5.0\times e^{-200*(x - 2/6)^2}+5.0 \times e^{-200*(x - 5/6)^2}$, $\; N_2(x, 0) = 5.0 \times e^{-200*(x - 1/6)^2)}+5.0\times e^{-200*(x - 4/6)^2}$, and (ⅲ) NK cells are positioned near the boundary around BV (Figure 17C-1). Initial conditions of cytokines (TGF-$\beta$ and IFN-$\gamma$; Figure 17D-1) and intracellular variables in the apoptotic signaling network (Bcl-2 and BAX) are set to be zero in the domain. Figure 17A–E shows spatial profiles of the densities of key players (cancer cells, N1/N2 TANs, NK cells, TGF-$\beta$, IFN-$\gamma$, Bcl-2 and BAX) at various time points ($t = 0, 7, 14, 60 \; days$). The relatively high proportion of N2 TANs in the tumor region in $\Omega_2$ (Figure 17I) promotes tumor growth in the early stage compared to the case on the left domain where tumor growth is suppressed from N1 TANs.

This accelerated cancer growth leads to an elevated level of TGF-$\beta$ (Figure 17D,G), which in turn suppresses NK cells activity (Figure 17C) and level of IFN-$\gamma$ (Figure 17D,F), which leads to up-regulation of Bcl-2 and down-regulation of BAX within the tumor (Figure 17E). Figure 17J shows time courses of tumor populations in $\Omega_1$ (red) and $\Omega_2$ (blue), respectively. Here, the population of tumor cells was calculated by integrating the tumor cell density over the spatial domains: $\hat{C}_1(t) = \int_{\Omega_1} C(x, t) \; dx$ (blue), $\hat{C}_2(t) = \int_{\Omega_2} C(x, t) \; dx$ (red). Other populations of other cells and levels of cytokines were calculated in a similar fashion: IFN-$\gamma$ ($\hat{S}_1(t) = \int_{\Omega_1} S(x, t) \; dx$, $\hat{S}_2(t) = \int_{\Omega_2} S(x, t) \; dx$ in Figure 17F)); TGF-$\beta$ ($\hat{G}_1(t) = \int_{\Omega_1} G(x, t) \; dx$, $\hat{G}_2(t) = \int_{\Omega_2} G(x, t) \; dx$ in Figure 17G); N1 TANs ($\hat{N}_{1, 1}(t) = \int_{\Omega_1} N_1(x, t) \; dx$, ($\hat{N}_{1, 2}(t) = \int_{\Omega_2} N_1(x, t) \; dx$ in Figure 17H); N2 TANs ($\hat{N}_{2, 1}(t) = \int_{\Omega_1} N_2(x, t) \; dx$, $\hat{N}_{2, 2}(t) = \int_{\Omega_2} N_2(x, t) \; dx$ in Figure 17H). Thus immune conditions such as the N1-N2 ratio create either the immunosuppressive or immune-active tumor microenvironment (Figure 17B,G). leading to different tumor growth rates (Figure 17J).

We now investigate the role of physical microenvironment in the spatial dynamics by considering two subdomains, $\Omega_-$ and $\Omega_+$: $\Omega_- = [0, 0.5]$, $\Omega_+ = [0.5, 1]$, $\Omega = \Omega_- \cup \Omega_+$. In order to reflect the different tissue compositions in $\Omega_-$ and $\Omega_+$, the diffusion coefficienst of all diffusible variables in $\Omega_-$ were set to be 100-fold smaller than those in $\Omega_+$, i.e., $D_w \vert_{\Omega_-} = \frac{D_w \vert_{\Omega_+}}{100}$, $w = C, \; N_1, \; N_2, \; K_{end}, \; K_{ex}, G, S$. The Initial conditions are prescribed as in Figure 17 except TANs. In this analysis, we assume that N2 TANs are immersed within the tumor as in a typical tumor microenvironment. Populations of cells and levels of cytokines on each subdomain were calculated as before, for instance, $\hat{C}_-(t) = \int_{\Omega_-} C(x, t) \; dx$, $\hat{C}_+(t) = \int_{\Omega_+} C(x, t) \; dx$ for IFN-$\gamma$ for tumor population and $\hat{S}_-(t) = \int_{\Omega_-} S(x, t) \; dx$, $\hat{S}_+(t) = \int_{\Omega_+} S(x, t) \; dx$ for IFN-$\gamma$ in $\Omega_-$ and $\Omega_+$, respectively. Figure 18A illustrates the different spatial profiles of cancer cells in $\Omega_-$ and $\Omega_+$ under different landscapes of tissue composition. The initial growth rate of the tumor in $\Omega_-$ is lower due to the smaller diffusion rate (Figure 18A-2). This essentially leads to a sequential chain of events in $\Omega_-$: Lower secretion rate of TGF-$\beta$ (Figure 18D; Figure 18F), lower transition to N2 TANs (Figure 18B; Figure 18G), lower intratumoral infiltration rate of NK cells (Figure 18C; Figure 18H), and lower apoptosis rate of cancer cells (Figure 18E). These feedback cycles between the tumor and immune cells via cytokines results in the lower level of tumor growth in $\Omega_-$ and faster proliferation of cancer cells in $\Omega_+$ in the later stage of cancer progression (Figure 18I).

In Figure 19, we investigate the effect of the combination therapy (NK cells + IFN-$\gamma$) on tumor growth by varying the injection locations within the computational domain $[0, 1]$. To assess the spatial influence of therapeutic administration on anti-tumor efficacy, the injection sites are divided into two distinct regions: The boundary region ($\Omega_b = [0, 0.1]\bigcup[0.9, 1]$) for drug infusion through blood vessels and a region at the center ($\Omega_c = [0.45, 0.55]$) for direct intratumoral injection. To systematically compare the effects of drug injections at different sites, we modify Eqs (3.24) and (3.26) by adding injection terms as follow:

Here, $I_{\Omega_i} ( = I_{\Omega_i}(x), \; i = c, b)$ is an indicator function (giving 1 or 0 based on the subdomain $x \in \Omega_i$) and the injection rates are now expressed as a weighted distribution between the two regions: $u_K v_1 I_{\Omega_c} + u_K v_2 I_{\Omega_b}, \quad u_S v_1 I_{\Omega_c} + u_S v_2 I_{\Omega_b}$, where the total injection amount remains constant, i.e., $\int_{\Omega} v_1 I_{\Omega_c} + v_2 I_{\Omega_b} \; dx = v_0$, constant. Figure 19A shows the cancer cell population at final time in response to injection of the therapeutic drugs with various injection rates $v_1 = 0, 0.1, \dots, 0.9, 1$ on $\Omega_c$. The best anti-tumor efficacy is obtained when 70% of these agents is injected in the center and the remaining 30% is injected on the periphery of the tumor ($v_1 = 0.7$; Figure 19A; blue in Figure 19B)) while the worst case comes from the injection of all resources only at the periphery ($v_1 = 0$) without intratumoral injection. Figure 19B shows the time courses of cancer cell population for three representative cases ($v_1 = 0.0, 0.7, 1.0$). While the anti-tumor agents injected on $\Omega_b$ ($v_1 = 0$) (Figure 19C-1) can effectively kill tumor cells only on the periphery of the tumor mass (Figure 19C-2), the intratumoral injection (Figure 19E-1) is very effective in killing cancer cells in the central, local part of the tumor mass with limited anti-tumor efficacy on the periphery (Figure 19E-2). The systemic injection on both the boundary and center (Figure 19D-1) leads to balanced suppression of tumor growth on both sites, leading to better anti-tumor efficacy (Figure 19B).

4.   Discussion

It is well known that many immune cells such as macrophages, TANs, and NK cells in a TME play a significant role in controlling tumor growth by interacting with a tumor via various cytokines and chemokines. In this paper, we investigated the role of NK cells in regulating the well-known opposite effect of N1 and N2 TANs on tumor growth by mathematical models in the presence and absence of various therapeutic approaches, including ionizing radiation. While the role of NK cells in regulation of the TME can be either tumor-promoting or tumor-suppressive ^[74], tumor-secreted TGF-$\beta$ was shown to suppress NK cell activity, thus blocking N2$\to$N1 transition by down-regulation of the IFN-$\gamma$ level ^{[59,60,71,72]}. For example, TGF-$\beta$ may suppress the function of NK cells in lung cancer by the following mechanism ^[165]: (ⅰ) Changes the receptor spectrum of NK cells in lung cancer patients ^{[166,167,168]} and (ⅱ) induction of NK cell polarization toward angiogenesis ^[169]. Immune cells may also suppress NK cell function by secreting immunosuppressive agents ^{[170,171,172]}. For instance, studies found that the relative population of Tregs ^[173] in the lung TME is higher than that in normal tissues and peripheral blood, both of which can secrete TGF-$\beta$ ^{[53,168,174,175]}. On the other hand, radiation therapy can enhance critical infiltration of the NK cells in TME, leading to better clinical outcomes ^{[73,74,75,76]}. The release of cytokines such as IFN-$\gamma$ was shown to activate NK cells and induce the chemotactic movement toward the TME ^{[79,80,81,82,83]}. NK cells play a significant role in regulation of tumor growth and there is an urgent need to develop NK cell-based cancer immunotherapies through recent technologies ^[176].

While the mutual interaction between TANs and tumor cells is important in TME, the dynamical system introduced in ^[85] is valid only when TANs are the dominant major players in the tumor-immune communication such as in the in vitro system and this ignores the significant roles of NK cells in the immune microenvironment. In this work, we developed a new model by extending the framework developed in ^[85] to take into account the key biological mechanisms of NK cells in regulation of the TAN-tumor interactions as well as neutrophenia, leading to finding the new dynamical regimes that were not observed in the reduced model. This extension also enabled us to perform the mathematical analysis of the more comprehensive tumor-immune interactions and development of new combination treatment strategies involving radiotherapy. In this work, we developed a mathematical model by a system of ordinary differential equations (ODEs) to investigate the effect of NK cells on controlling the TME including N1 and N2 TANs and RT-induced neutropenia. Mathematical analysis showed the complex behaviors of phenotypic switches between N2 TANs and N1 TANs, thus nonlinear tumor growth, in response to fluctuating TGF-$\beta$ and NK cells (Figure 3). We developed several therapeutic strategies of slowing down tumor growth by using the network structure of TANs and STAT signaling. In particular, we developed various injection schedules of two drugs (IFN-$\gamma$ and TGF-$\beta$ inhibitor) and NK cells with a goal of determining the most effective treatment regimen with the lowest possible dose (Figures 6–11). Due to the competitive nature of three cell death programs (apoptosis via Bcl2-BAX, N1-mediated killing, NK-mediated killing) in tumor cell killing, synergistic effects were observed within a particular schedule bands and injection orders. For instance, injection of IFN-$\gamma$ and NK cells in the minimum amount required to induce the anti-tumorigenic phenotype (Figures 7 and 9), followed by induction of the apoptosis phenotype (Figure 11) can induce better strategic outcomes in reducing the tumor size.

Radiation therapy can induce upregulation of TGF-$\beta$, a tumor-promoting factor, despite its efficient killing of tumor cells ^{[15,81,82,83]}. The combination of radiation and the TGF-$\beta$ inhibitor is effective in mitigating the pro-cancer effect of TGF-$\beta$ (Figure 15) due to enhancement of NK cell activities by radiation. Under the same schedule of radiotherapy and TGF-beta inhibitor ^{[109,110,134,135,136]}, the phenotypic changes of tumor-associated neutrophils, thus anti-tumor efficacy, varied depending on the timing of TGF-beta inhibitor injection, leading to a better treatment strategy (Figure 15).

Our findings in the PDE-based model highlight the complex interplay between tumor diffusion dynamics and immune regulation in a tumor microenvironment where the different spatial distribution of N2-N2 ratio or different tumor dispersal speed and NK infiltration in a harsh environment affect the immunosuppressive influence of N2 TANs and NK-mediated cancer killing, and the overall tumor-immune dynamics (Figures 17 and 18). We also found that the strategic approach of injecting therapeutic agents at different locations (infusion through blood vessels or direct intratumoral injection) may lead to better clinical outcomes (Figure 19), as shown in oncolytic virus therapy ^[102,177].

This work should be a starting point of finding the fundamental mechanism of biochemical interactions of N1/N2 TANs-tumor-NK cells. In the TME, there are many other critical factors that we did not consider in this work, such as M1/M2 TAMs ^[9], signaling networks ^[178], angiogenesis ^{[179,180,181]}, tumor-associated fibroblasts ^[182], tumor ECM remodeling and restructuring ^[3,182,183], CSF-1 ^[8,184], and NET ^{[185,186,187]} which can also play a pivotal role in regulation of cancer progression. In addition, mutations in STAT1 and STAT3 molecules (loss of function (LOF) or gain of function (GOF)) were not considered ^[188,189]. While we focused on TANs due to their rapid phenotypic plasticity and influence on tumor-immune interactions ^[15], M1/M2 TAMs are another critical component of the tumor microenvironment ^[9,190]. Both TAMs and TANs contribute to immunosuppressive and pro-tumorigenic pathways, often governed by overlapping cytokines such as TGF-$\beta$ and IL-6 ^[6,191]. In future work, we will extend this model to incorporate TAM dynamics and their interactions with TANs, providing a more comprehensive framework for understanding tumor progression and immune modulation. These factors could impact tumor invasion and growth. We plan to investigate the role of these important factors in near future work.

Use of AI tools declaration

The authors declare they have not used artificial intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. RS-2024-00347728). This paper was written as part of Konkuk University's research support program for its faculty on sabbatical leave in 2024.

Conflict of interest

Yangjin Kim is an editorial board member for Mathematical Biosciences and Engineering and was not involved in the editorial review or the decision to publish this article. The authors declare no conflict of interest.

Appendix

A.   Nondimensionalization

We utilized dimensionless expression when deriving the results presented in the main text. Consequently, in our mathematical model, we proceed with a process of non-dimensionalization to remove units in all variables. The reference value for each variable is specified in Table A1.

Let variables $C(t), N_1(t), N_2(t), G(t), L(t), S(t)$, $K_{end}(t)$, $K_{ex}(t)$, $S_1(t), S_3(t), B(t)$, and $X(t)$ be densities of cancer cells, N1 and N2 TANs, TGF-$\beta$, TGF-$\beta$ inhibitor, IFN-$\gamma$, endogenous and exogenous NK cells, and activities of STAT1, STAT3, Bcl-2, and BAX, respectively, at time $t$. The governing equations for the rate change of all variables are

By defining following and using reference values of each variable in Table A1

we get the dimensionless equations as follows:

B.   Parameter estimation

The parameter estimation for the N2/N1 TAN modules and STAT signaling module is conducted based on biological observations, literature data, and mathematical assumptions to ensure an accurate representation of tumor-immune interactions. Key parameters were chosen to reflect the autocatalytic effects, inhibition mechanisms, and transition dynamics. The estimation process aimed to capture the essential behaviors of the system, including positive feedback loops, mutual inhibition, and phenotypic transitions, while maintaining biological realism. Below, we outline the key parameter values and their derivations.

$\mu_{N2}$, $\mu_{N1}$: The half-life of TANs within the circulation ($T_{1/2}$) was reported to be 20 hours ^[193,194]. By taking $T_{1/2} = 20\; h$, we get $\mu_{N2} = \mu_{N1} = \frac{ln(2)}{T_{1/2}}\sim0.035\; h^{-1}$.

$\lambda_6$: The reference value for N2 TANs is taken as $N_2^* = 1.0\times10^{-4}\; g/cm^3$, with a signal source of N2 TANs by IL6 ($3.5\times10^{-8}\; g/cm^3\cdot h$), we estimate the dimensionless value of signaling strength to the N2 TANs to be $\lambda_6 = \frac{3.5\times10^{-8}\; g/cm^3\cdot h}{1\; h\cdot1.0\times10^{-4}\; g/cm^3} = 3.5\times10^{-4}$.

$k_1$, $k_3$: We take the dimensionless parameter value $k_1 = 0.14$ by assuming that the autocatalytic rate of N2 TANs is at least 4-fold larger than its negative contribution from natural decay ($\mu_{N2}N2^*$) in the absence of the inhibitory by N1 TANs. This assumption is based on the ratio $\frac{k_1}{\mu_{N2}} = \frac{0.14}{0.035} = 4$, ensuring that the self-enhancement of N2 TANs significantly outweighs their decay dynamics. In a similar fashion, we get the autocatalytic rate of N1 TANs, $k_3 = 0.14$.

$k_2$, $k_4$: The Hill-type coefficients ($k_2, k_4$) are fixed to be unity.

$\mu_{S1}$, $\mu_{S3}$: The half-life of STAT family within the circulation ($T_{1/2}$) was reported to be 24 hours ^[203,204]. By taking $T_{1/2} = 24\; h$, we get $\mu_{S1} = \mu_{S3} = \frac{ln(2)}{T_{1/2}}\sim0.0289\; h^{-1}$.

$\mu_{B}$: The half-life of Bcl-2 within the circulation ($T_{1/2}$) was reported to be 20 hours ^[205,206]. By taking $T_{1/2} = 20\; h$, we get $\mu_{B} = \frac{ln(2)}{T_{1/2}}\sim0.0347\; h^{-1}$.

$\mu_{X}$: The half-life of BAX within the circulation ($T_{1/2}$) was reported to be 4.8 hours ^[207,208]. By taking $T_{1/2} = 4.8\; h$, we get $\mu_{X} = \frac{ln(2)}{T_{1/2}}\sim0.1445\; h^{-1}$.

$\lambda_2$, $\lambda_4$: The reference value for Bcl-2 is taken as $B^* = 0.01\; \mu M = 2.6\times10^{-7}\; g/cm^3$ (molecular weight of Bcl-2 is $26\; kDa$), with a signal source of Bcl-2 ($1.508\times10^{-9}\; g/cm^3\cdot h$), we estimate the dimensionless value of signaling strength to the Bcl-2 to be $\lambda_2 = \frac{1.508\times10^{-9}\; g/cm^3\cdot h}{1\; h\cdot2.6\times10^{-7}\; g/cm^3} = 0.0058$. Similarly, by taking the reference value of BAX ($X^* = 351\; \mu M = 7.371\times10^{-3}\; g/cm^3$), and taking the signal source of BAX ($4.3\times10^{-5}\; g/cm^3\cdot h$), we estimate the dimensionless value of signaling source to be $\lambda_4 = \frac{4.3\times10^{-5}\; g/cm^3\cdot h}{1\; h\cdot7.371\times10^{-3}\; g/cm^3} = 0.0058$.

$a_5$, $\lambda_3$: We take the dimensionless value $a_5 = 0.0289$ by assuming that the autocatalytic rate of Bcl-2 is same as its negative contribution ($\mu_{S1}B^*$) in the absence of the inhibitory by STAT1. This assumption is based on the ratio $\frac{a_5}{\mu_{S1}} = \frac{0.0289}{0.0289} = 1$, ensuring that the self-enhancement of Bcl-2 significantly same their decay dynamics. Additionally, we assume that the STAT3-induced signaling rate is at least 1.2-fold larger than its negative contribution ($\mu_{S1}B^*$). We estimate the dimensionless parameter value of STAT3-induced signaling rate $\lambda_3 = 1.2 \times \mu_{S1} B^* = 0.0347$.

$a_1$, $a_3$, $a_7$: The dimensionless parameter $a_1$ is estimated as $a_1 = 4.0\times\mu_{S1}S_1^* = 0.1156$ by assuming that the autocatalytic rate $a_1$ of STAT1 is at least 4-fold greater than its negative contribution, in the absence of inhibition by STAT3. Using the same assumption, we obtain: $a_3 = 4.0\times\mu_{S1}S_1^* = 0.1156$ and $a_7 = 4.0\times\mu_{S1}S_1^* = 0.1156$.

$a_2$, $a_4$, $a_6$, $a_8$: The Hill-type coefficients are fixed to be unity.

$\mu_{S}$: The half-life of IFN-$\gamma$ within the circulation ($T_{1/2}$) was reported to be 5 hours ^[197]. By taking $T_{1/2} = 5\; h$, we get $\mu_{S} = \frac{ln(2)}{T_{1/2}}\sim0.1386\; h^{-1}$.

$\mu_{G}$: The half-life of TGF-$\beta$ within the circulation ($T_{1/2}$) was reported to be 24 hours ^{[8,91,100,195]}. By taking $T_{1/2} = 24\; h$, we get $\mu_{G} = \frac{ln(2)}{T_{1/2}}\sim0.0289\; h^{-1}$.

$\mu_{L}$: The half-life of TGF-$\beta$ inhibitor within the circulation ($T_{1/2}$) was reported to be 3 hours ^[196]. By taking $T_{1/2} = 3\; h$, we get $\mu_{L} = \frac{ln(2)}{T_{1/2}}\sim0.231\; h^{-1}$.

$\mu_{K}$: The half-life of NK cells within the circulation ($T_{1/2}$) was reported to be 170 hours ^[198,199]. By taking $T_{1/2} = 170\; h$, we get $\mu_{K} = \frac{ln(2)}{T_{1/2}}\sim0.0041\; h^{-1}$.

$\mu_{N}$: The half-life of neupogen within the circulation ($T_{1/2}$) was reported to be 4 hours ^[200,201]. By taking $T_{1/2} = 4\; h$, we get $\mu_{N} = \frac{ln(2)}{T_{1/2}}\sim0.1773\; h^{-1}$.

C.   Stability analysis in key submodules

In this section, we focus on two key interacting subsystems that primarily govern tumor phenotype regulation in the tumor microenvironment (TME): (ⅰ) N1/N2 TANs dynamics in the presence NK cells in response to varying TGF-$\beta$, and (ⅱ) intracellular, apoptosis signaling pathway in the presence of NK cells. These two subsystems are linked through shared cytokines and NK cell interactions, but their stability properties with respect to cytokine and NK cell levels in the subsystem are a crucial part of comprehensive understanding of the whole system due to their distinct functional roles. This approach aligns with a quasi-steady-state assumption, where we evaluate stability by fixing certain variables at physiologically relevant levels to determine the possible stable phenotypic outcomes.

The dynamics of the tumor-associated neutrophils (TANs) is essentially governed by ordinary differential equations:

Then, the equilibrium state of Eqs (A27)–(A29) satisfies

for low ($G = 0$), intermediate ($G = 0.4$), and high ($G = 1$) values of $G$. The Jacobian matrix, $\mathcal{J}$, is given by:

Then, the characteristic polynomial is given by

where ${\mathbf I}$ is the identity matrix and defines $A = -\frac{2k_1k_2^2\alpha N_1^E}{(k_2^2+\alpha(N_1^E)^2)^2}$ and $B = -\frac{2k_3k_4^2\beta N_2^E}{(k_4^2+\beta(N_2^E)^2)^2}$. We get two real eigenvalues ($\Lambda_1 = -\mu_S$ and $\Lambda_2 = -\mu_K$) and remaining two real eigenvalue $(\Lambda_{3, 4} = \frac{1}{2}\left(-(\mu_{N2}+\mu_{N1})\pm\sqrt{(\mu_{N2}-\mu_{N1})^2+4AB}\right))$. Here, since $\mu_{N2} = \mu_{N1}$, we have $\Lambda_{3, 4} = -\mu_{N2}\pm\sqrt{AB}$ (or $-\mu_{N1}\pm\sqrt{AB}$). Since all of our parameters are positive, $AB > 0$. Therefore, $\Lambda_4\left( = -\mu_{N2}-\sqrt{AB}\right) < 0$, and we can determine the stability of the equilibrium point by examining the sign of $\Lambda_3\left( = -\mu_{N2}+\sqrt{AB}\right)$. If $AB < \mu_{N2}^2$, the $\Lambda_3 < 0$ and the equilibrium is stable. On contrary, if $AB > \mu_{N2}^2$, the $\Lambda_3 > 0$ and the equilibrium is unstable. By substituting parameters in the model, we get The Jacobian matrix, $\mathcal{J}$, is given by:

where $A = -\frac{21 N_1^E}{50((1+1.5(N_1^E)^2)^2)}$ and $B = -\frac{14 N_2^E}{50((1+(N_2^E)^2)^2)}$. When $G = 0$, we get $\Lambda_3 = -0.0245$, leading to a stable steady state. Also, the high level of TGF$\beta$ ($G = 1$) leads to a stable state ($\Lambda_3 = -0.0185 < 0$). For an intermediate level of TGF$\beta$ ($G = 0.4$ in the the bi-stability window), we get $\Lambda_3 = -0.0095$ when $(N_1^E, N_2^E, S^E, K_{end}^E)\approx(2.7635, 0.7311, 0.1569, 0.7852)$, $\Lambda_3 = 0.0053$ when $(N_1^E, N_2^E, S^E, K_{end}^E)\approx(1.1780, 1.7081, 0.1569, 0.7852)$, and $\Lambda_3 = -0.0090$ when $(N_1^E, N_2^E, S^E, K_{end}^E)\approx(0.4764, 3.3939, 0.1569, 0.7852)$, leading two stable steady state and one unstable steady state.

We first examine $\Lambda_3$, whose sign is determined by $A$ and $B$ while $\Lambda_1 < 0, \; \Lambda_2 < 0, \; \Lambda_4$ for various TGF-$\beta$ values. However, the sign of $\Lambda_3( = -\mu_{N2}+\sqrt{AB})$ can be either positive or negative, depending on the equilibrium point that is determined by $G$. Figure A1A shows the third eigenvalue ($\Lambda_3$) for various TGF-$\beta$ values. Figure A1B,C shows the steady state of the densities of N1 TANs and N2 TANs for various $G$s, respectively. Here, the red (or blue) circles represent unstable (or stable) steady states in Figure A1. For example, the steady states on the middle branch of the bifurcation curves for N2 TANs (Figure A1B) and N1 TANs (Figure A1C) are unstable due to the positive eigenvalues (Figure A1A).

We also investigated the bifurcation curves for NK cells in the main text (Figure 4 in main text). To investigate the bifurcation curves for NK cells, it is necessary to fix the level of TGF-$\beta$ ($G$) and obtain the eigenvalue using equation of tumor-associated neutrophils (N1 TANs and N2 TANs) and IFN-$\gamma$:

The Jacobian matrix, $\mathcal{J}_1$, is given by:

Then, the characteristic polynomial is given by

where $A = -\frac{2k_1k_2^2\alpha N_1^E}{(k_2^2+\alpha(N_1^E)^2)^2}$ and $B = -\frac{2k_3k_4^2\beta N_2^E}{(k_4^2+\beta(N_2^E)^2)^2}$. We get one real eigenvalue ($\Lambda_1 = -\mu_S$) and remaining two eigenvalue ($\Lambda_{2, 3} = -\mu_{N2}\pm\sqrt{AB}$ since $\mu_{N2} = \mu_{N1}$). Since all of our parameters are positive, $AB$ is also always positive as well. Therefore, $\Lambda_3\left( = -\mu_{N2}-\sqrt{AB}\right)$ is always negative, and we can determine the stability of the equilibrium point by examining the sign of $\Lambda_2\left( = -\mu_{N2}+\sqrt{AB}\right)$. In this case, $\Lambda_{1, 3}$ are always negative real number for various NK cells at the intermediate level of TGF-$\beta$ ($G = 0.5$). Figure A2A shows the second eigenvalue ($\Lambda_2$) for various NK cells at $G = 0.5$. Figure A2B,C shows the steady state of the densities of N1 TANs and N2 TANs for various NK cells at $G = 0.5$. Also, the second eigenvalue is positive at middle branch steady state of the bifurcation curves and the equilibrium is unstable.

Finally, we investigate stability analysis of intracellular signaling (STAT signaling) for various NK cells.

At the equilibrium ($S_1^E, S_3^E, B^E, X^E, S^E$) of Eqs (A42)–(A46), we have

Then, the characteristic polynomial is given by

where $A = -\frac{2a_1a_2^2\gamma S_3^E}{(a_2^2+\gamma(S_3^E)^2)^2}$ and $B = -\frac{2a_3a_4^2\delta S_1^E}{(a_4^2+\delta(S_1^E)^2)^2}$. We get three negative real eigenvalues ($\Lambda_1 = -\mu_S < 0$, $\Lambda_2 = -\mu_B < 0$, $\Lambda_3 = -\mu_X < 0$) and two remaining eigenvalues $\Lambda_{4, 5}$ ($\Lambda_{4, 5} = -\mu_{S1}\pm\sqrt{AB}$ since $\mu_{S1} = \mu_{S3}$). Since all of our parameters and steady states are positive, $AB > 0$. Therefore, $\Lambda_{1, 2, 3}$ and $\Lambda_5 ( = -\mu_{S1}-\sqrt{AB})$ are always negative, and we can determine the stability of the equilibrium point by examining the sign of $\Lambda_4( = -\mu_{S1}+\sqrt{AB})$. Figure A3A shows the fourth eigenvalue ($\Lambda_4$) for various $K_{end}$. The red (or blue) circles represent the unstable (or stable) steady states. Figure A3B,C shows steady states of apoptotic cells (Bcl-2 and BAX, respectively) for various $K_{end}$. The positive fourth eigenvalue along the middle branch of the bifurcation curves indicates the unstable status of the intracellular signaling network.

D.   Glossary of selected terms and abbreviations

Angiogenesis: the development of new blood vessels, especially in tissues where circulation has been impaired by trauma or disease, e.g., cancer.

Apoptosis: programmed cell death characterized by nuclear breakdown and removal of remains by phagocytes.

BAX: Bcl2 Associated X.

Bcl-2: B-cell lymphoma 2.

BED: Biologically effective dose.

c-MYC: c-Myelocytomatosis oncogene.

CSF-1: Colony stimulating factor 1.

Cytokine: extracellular signaling protein that acts as a local mediator in cellecell communication. Those involved in taxis are sometimes called chemokines.

ECM: Extracellular matrix.

GOF: gain of function.

Growth factor: an extracellular signaling molecule that stimu- lates a cell to grow or proliferate, e.g., vascular endothelial growth factor (VEGF) and fibroblast growth factor (FGF).

IFN-$\gamma$: Interferon gamma.

IGRT: Image guided radiation therapy.

IL-6: Interleukin 6.

IMRT: Intensity modulated radiation therapy.

IR: Ionizing radiation.

LOF: loss of function.

LUAD: Lung adenocarcinoma.

LUSC: Lung squamous cell carcinoma.

Metastasis: the process by which cancer spreads from the site of initiation of the primary tumor to distant locations in the body. This occurs via either the circulatory system or the lymphatic system.

MHC: Major histocompatibility complex.

mTORC1: mammalian target of rapamycin complex 1.

NETs: Neutrophil extracellular traps.

Neutropenia: Low neutrophil count.

NF$\kappa$B: Nuclear factor kappa-light-chain-enhancer of activated B cells.

NK cells: Natural killer cells.

NLR: Neutrophil to lymphocyte ratio.

NSCLC: Non-small cell lung cancer.

ODE: ordinary differential equation.

PDE: partial differential equation.

RT: Radiation therapy.

SREBP: Sterol Regulatory Element-Binding Protein.

STAT: Signal Transducer and Activator of Transcriptions.

TAMs: Tumor-associated macrophages.

TANs: Tumor-associated neutrophils.

TGF-$\beta$: Transforming growth factor-beta.

TILs: Tumor-infiltrating lymphocytes TME: Tumor microenvironment.

TNF-$\alpha$: Tumor necrosis factor-$\alpha$.

VMAT: Volumetric modulated arc therapy.