Role of tumor-associated neutrophils in regulation of tumor growth in lung cancer development: A mathematical model

Neutrophils display rapid and potent innate immune responses in various diseases. Tumor-associated neutrophils (TANs) however either induce or overcome immunosuppressive functions of the tumor microenvironment through complex tumor-stroma crosstalk. We developed a mathematical model to address the question of how phenotypic alterations between tumor suppressive N1 TANS, and tumor promoting N2 TANs affect nonlinear tumor growth in a complex tumor microenvironment. The model provides a visual display of the complex behavior of populations of TANs and tumors in response to various TGF-β and IFN-β stimuli. In addition, the effect of anti-tumor drug administration is incorporated in the model in an effort to achieve optimal anti-tumor efficacy. The simulation results from the mathematical model were in good agreement with experimental data. We found that the N2-to-N1 ratio (N21R) index is positively correlated with aggressive tumor growth, suggesting that this may be a good prognostic factor. We also found that the antitumor efficacy increases when the relative ratio (Dap) of delayed apoptotic cell death of N1 and N2 TANs is either very small or relatively large, providing a basis for therapeutically targeting prometastatic N2 TANs.


Introduction
Lung cancer is the leading cause of cancer mortality worldwide, with an approximate 1.6 million deaths each year [1]. The most common (*85%) form of lung cancer in patients is nonsmall cell lung cancer (NSCLC), of which lung squamous cell carcinoma (LUSC) and lung adenocarcinoma (LUAD) are the most common subtypes [2]. Various groups of myeloid cells have been known to promote tumor development by direct inhibition of immune responses [3], as well as by secreting growth factors, angiogenic factors, or matrix-degrading enzymes [4,5]. For example, tumor-associated macrophages (TAMs), also known as M2 macrophages [3], have been shown to promote tumor growth [6,7]. There is growing evidence suggesting that PLOS ONE | https://doi.org/10.1371/journal.pone.0211041 January 28, 2019 1 / 40 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 neutrophils play a major role in tumor progression from establishment of tumor formation and throughout the progression to the malignant state [8][9][10][11][12]. For example, tumor associated neutrophils (TANs) have been associated with poor prognosis in many cancers including metastatic melanoma [13], bronchoalveolar carcinoma [14], and renal carcinoma [15]. Like TAMs, TANs infiltrate tumor tissue and can have two differential states in cancer progression [8,9,16]: (i) an antitumorigenic role (called N1) (ii) promoting tumor progression (called N2). How these two phenotypes are regulated is largely unknown but many experimental and clinical findings suggest the significant potential of therapeutic targeting of the prometastatic role of TANs [17]. TGF-β has been identified as a major cytokine within a tumor that skews neutrophil differentiation toward the N2 phenotype [16,18,19], while TGF-β blockade and type-1 IFN (α, β, ω) treatment are known to shift the balance toward the N1 phenotype [20,21]. IFN-β in tumor microenvironment can directly suppress tumor growth [22] by interacting with p53 [23-25]. IFN treated neutrophils were shown to upregulate PD-L1 and suppress T-cell proliferation [26]. After binding to interferon receptor type 1, IFNAR1 and IFNAR2, Type 1 IFN-β signals through TYK2 and JAK1, which in turn phosphorylate STAT family members (STAT1, STAT2, STAT3, and others) and activate its downstream functions to stimulate anti-tumor activities [27]. For example, vesicular stomatitis virus expressing IFN-β was shown to enhance anti-tumor immune responses in a murine model of NSCLC [28]. It is well established that cancer associated fibroblasts (CAFs) can promote tumor growth, aggressive invasion, and metastasis through mutual interaction in the tumor microenvironment [29,30]. Fibroblastsecreted IFN-β was also able to restrict expression of the p53 RNA stabilizer, WIG1, and bring down mutant p53 RNA levels, thus suggesting an alternative therapeutic agent for mutant p53 positive lung cancer patients [31]. There are multiple levels of crosstalk between neutrophils and many cells including other immune cells and Th17 cells [32]. Neutrophils may express many important factors such as IL-6, IL-17A, IL-17F and IFNγ [33]. How neutrophils are induced by a tumor is still poorly understood. It is well known that tumor cells interact with stromal cells such as fibroblasts, immune cells (neutrophils, macrophages, Th17, Tregs, T cells), and cytokines in the tumor microenvironment, and that these complex interactions play a critical role in tumor initiation, growth, angiogenesis, and metastasis. The mutual interactions between a tumor and immune system involving TANs are summarized in  Table 1.
In this work, we employed the framework of the differential equation model in order to take into account the role of N1 and N2 TANs in the regulation of tumor-neutrophil interactions and tumor growth. Based on a schematic diagram in Fig 2, the model consists of a system of ordinary and delay differential equations involving the following variables: CðtÞ ¼ density of the N2 complex at time t; IðtÞ ¼ density of the N1 complex at time t; TðtÞ ¼ density of tumor cells at time t: TGF-β and IFN-β play a central role in the critical transition between N1 and N2 TANs, resulting in either tumor-promoting or tumor-suppressing microenvironment. We summarize in Fig 3 the fundamental, phenotypic spectrum underlying how the oncogenic and tumor suppressive properties of TANs arise in relation to TGF-β, IFN-β, and delayed apoptotic cell death mechanisms. This model is then used (i) to investigate how TGF-β and IFN-β regulate the balance between N1 and N2 phenotypes by complex mutual interactions, (ii) to investigate how up-or down-regulation of these N1 and N2 phenotypes affect tumor growth and anti-tumor densities of the N2 module (N 2 ) and N1 module (N 1 ) leads to a system of ordinary differential equations (ODEs) as follows: The first term in Eq (1) represents the signaling from N2-promoting molecules such as IL-6 and TGF-β. Similarly, the first term in Eq (2) represents the signaling from N1-promoting molecules such as IFN-β. The third terms in Eqs (1) and (2) indicate natural decay/consumption of N2 and N1 complexes, respectively. The second term in Eq (1) represents autocatalytic Table 1. Summary of interaction networks (Fig 1) in invasion processes in lung cancer development.
The mass balance of the tumor volume (T) leads to an ordinary differential equation as follows: where r is a growth rate of the tumor and T 0 is the carrying capacity of the tumor population. All simulations were performed using Matlab (Mathworks) for ordinary differential equations (ode45) and delay differential equations (dde23).  [60], and 17.3 h in terminal cancer patients with glioblastoma and chronic lymphocytic leukemia [61]. For example, Andzinski et al. [62] recently reported that the life span of TANs can be remarkably prolonged in tumor microenvironment where IFN-β activity is suppressed [63]. Based on these biological observations in the tumor microenvironment, neutrophils in tumor tissue, we take the half-life of about T 1/2 = 20 h. In general, we can obtain the decay rate of a variable when we know the half-life. In more detail, by solving a typical differential equation of the form dy dt ¼ À ly with an initial condition y(0) = y 0 and a delay rate λ, and comparing the solution y(t) = y 0 e −λt at t = T 1/2 , we can obtain the decay rate l ¼ lnð2Þ
G � : The TGF-β level in the lung cancer microenvironment was measured to be in the range of 400 − 2, 000 pg/mL [38]. We take G � = 1, 000 pg/mL. S � : Deng et al. [68] investigated the effect of STING-dependent cytosolic DNA sensing on the IFN-β-dependent anti-tumor efficacy with the control IFN-β concentration of 10 ng/mL and IFNγ value of 20 ng/mL. We take S � = 10 ng/mL. Table 2 summarizes all the essential parameter values in the system (5)-(7).

Dynamics of the model
We first investigated the local dynamics of the system (5) and (6). By taking different initial values of N1 and N2 complexes and numerically solving Eqs (5) and (6), one can obtain solutions (C(t), I(t)) at time t and plot the relative locations of solutions in the C − I domain. In Fig  5 we illustrate three different flow patterns of the system (5) and (6) near the equilibrium point (steady state (SS); circles) in response to low (G = 0.1 in Fig 5A), intermediate (G = 0.4 in Fig  5B), and high (G = 1.0 in Fig 5C) TGF-β levels in the C − I phase diagram. We note that if the equilibrium point is stable, the solution converges to the stable SS (called attractor). On the other hand, when the equilibrium point is unstable, the solution does not converge to the unstable SS and the solution is repelled by the unstable SS (called repeller). By taking the thresholds, C th (= 1.5) for the N2 level and I th (=1.5) for the N1 level, we shall define the anti- tumorigenic (P a ) and tumorigenic (P t ) phases as follows: P a ¼ fðC; IÞ 2 R 2 : C < C th ; I > I th g P t ¼ fðC; IÞ 2 R 2 : C > C th ; I < I th g: Fig 5D illustrates the anti-tumorigenic (low activities of the N2 module, high activities of N1 TANs) and tumorigenic (high activities of the N2 complex, low activities of the N1 complex) in a C − I plane. A low TGF-β level (G = 0.1) induces only one stable SS ( Fig 5A) where suppression of N2 activities and promotion of N1 activities induce the anti-tumorigenic status (P a ). In this case, the environment eventually induces the anti-tumorigenic phenotype regardless of initial status of N1 and N2 complexes. In contrast, the increased activity of N2 TANs, and decreased activity of N1 TANs, leading to the tumorigenic phase (P t ), are induced under a high TGF-β condition (G = 1.0; Fig 5C; one stable SS). In this case, the environment eventually induces the tumorigenic phenotype regardless of initial status of N1 and N2 complexes. For an intermediate level of TGF-β signals (G = 0.4), there exist three equilibria: two stable SS (two filled circles; one in P t and one in P a ) and one unstable SS (empty circle) in the middle. This leads to a bi-stable system shown in Fig 5B. In this case, the environment will induce either anti-tumorigenic or tumorigenic phenotype depending on the initial status of N1 and N2 complexes. From these observations, we anticipate to see a bifurcation curve, especially hysteresis, i.e., the N2 level with TGF-β signal parameter G.

Effect of TGF-β on transitions between N1/N2 TANs and tumor growth
When the system (5) and (6) is in equilibrium, we can solve C as a function of the TGF-β signal (G). Fig 6 shows the graphs C = C(G) and I = I(G) as a hysteresis. While the upper and lower branches are stable, the middle branch is unstable (recall the G-dependent stability in Fig 5). If G is small, then the system is in the lower branch (C low, I high), leading to the anti-tumorigenic phase (P a ). This situation continues to hold as G is increased until it reaches the right knee point of the curve (when G = 0.53). At this point, the N2 level jumps to the high branch, with an elevated level of N2 TANs, leading to the tumorigenic phase (P t ). As G is decreased, the N2 level remains elevated, until G is decreased to the left knee point of the curve (when G = 0.33), at which time the N2 level jumps down to the lower branch, and the system returns to the cancer-free mode (P a ). We conclude that the effect of G is history dependent: when G is at an intermediate level (0.33 = G m < G < G M = 0.53; bi-stable mode in Fig 5B), the system induces the anti-tumor phase P a if G is in an increasing mode from the low G value, and leads to the tumor-promoting phase P t if G is in decreasing mode from the high G level. The bifurcation diagram in Fig 6 suggests that a state (G, C) with C > C th = 1.5 will be moved by the dynamical system (5) and (6) into the upper stable steady state branch, resulting in prevailing P t -status. On the other hand, if C < C th = 1.5, then (G, C) will stay in the lower stable state branch, leading to P a -status. The size of the bi-stability window (W G = [G m , G M ]) depends on other parameters and may even disappear for some parameter sets. Fig 7A shows how main variables (C, I) in the system respond to the slowly increasing TGF-β level (G). The solution C(t) (solid red line) initially follows the P a -status along the lower branch of the bifurcation curve until the right knee point of the bifurcation curve and transits to the upper branch solution, leading to the P t -mode. This trajectory of the solutions (C, I) along G is illustrated in Fig 7C. This illustrates that when G is increasing from a low value, the system will stay longer in the P a -phase. On the other hand, when G is slowly decreased from a high value, the solution C(t) follows the upper branch of the bifurcation curve in Fig 6 until the left knee point and drops down to the lower branch, staying longer in the P t -phase including the bistability region W G (Fig 6). The corresponding trajectory of the solutions (C, I) along G is shown in Fig 7D. In Fig 8, we test the effect of fluctuating TGF-β on tumor growth. A time-dependent TGF-β level was assigned for a periodic supply of TGF-β in tumor microenvironment as follows ( Fig 8A): In response to the initial increasing TGF-β level, the trajectory of the solution C(t) follows the lower branches of the hysteresis bifurcation loop (thick gray curve) and jumps to the upper branch ( Fig 8B). However, a decrease in G around t = 250 switches the moving direction of the solution and it follows the upper branch in the bi-stability region until it jumps down to the lower branch near the left knee of the hysteresis curve. This flow of the solution (G(t), C(t)) was marked in red arrows with initial position on the left lower corner ( Fig 8B). The blue dashed curve and blue arrows represent the solution (G(t), I(t)) and its flow, respectively, in the opposite side ( Fig 8B). The detailed time courses of solutions (C(t), I(t)) in response to the  Table 2. https://doi.org/10.1371/journal.pone.0211041.g007

Role of TANs in cancer progression
periodic G input are shown in Fig 8C. The specific N1/N2 phenotypic transitions occur at time t = t 1 , t 2 , t 3 , t 4 . For example, the initial N1-dominant P a -status (white region in Fig 8C) switches to a N2-dominant P t -mode (pink region in Fig 8C) in the tumor microenvironment at time t = t 1 due to increasing TGF-β levels. The TAN phenotypic ratio, described by C Kþg 1 I in the second term of Eq (7), and N1 TAN immunity on tumor cells, described by δIT in the third term of Eq (7), determine either promotion or suppression of tumor growth in the model. Fig 8D shows time courses of tumor volumes in response to either fixed or varying TGF-β: (i) high TGF-β (G = 0.9, blue circle), (ii) low TGF-β (G = 0.1, dashed), and (iii) fluctuating TGF-β levels given in Fig 8A. The high TGF-β level induces N2-dominant P t -phase ( Fig  5C), therefore, leading to N2-mediated, faster tumor growth (blue circle). When G is low, the system induces the N1-dominant P a microenvironment (Fig 5A), resulting in N1-mediated, slower tumor growth (black dashed). On the other hand, fluctuating TGF-β induces alternating transitions between P a (white region) and P t (pink region) phases at t = t 1 , t 2 , t 3 , t 4 ( Fig  8C), which induces slow (white band in Fig 8D) or fast (pink band in Fig 8D) tumor growth in P a and P t -regions, respectively, and eventually leads to an intermediate tumor volume at final time ( Fig 8D). These results are in good agreement with experimental observations where N1 TANs significantly reduced the size of a tumor relative to N2 TANs [8]. The N2-to-N1 ratio (N21R) index and tumor volume corresponding to these three cases are shown in Fig 8E. An increase in the TGF-β level leads to an increase in the N2-to-N1 ratio, which in turn increase the tumor size. In Fig 8F, we investigate the correlation between N21R and the tumor size. The tumor volume and N21R were calculated from tumor and TAN populations in response to various TGF-β stimuli (G = 0 − 0.9) with two different initial conditions: C(0) = 1, I(0) = 2.5 in red asterisk; C(0) = 2.5, I(0) = 1 in blue circle. The case of the fluctuating G(t) was marked in black square. As N21R is increased, the tumour volume is increased, leading to a positive  correlation between N21R and the tumor size. The neutrophil-to-lymphocyte ratio (NLR) has been highly associated with cancer progression [69,70], becoming biomarkers for various cancers including non-small-cell lung cancer [71,72] and cervical cancer [73]. In our simulation, N21R is equivalent to NLR and may be suggested to be a prognostic factor. See S1 Appendix for analysis of the model. We investigated effects of N1 immunity (δ) on cancer cell killing. Depending on various values of δ, the killing rate of cancer cells by N1 TANs, time courses of tumor volume in response to fluctuating TGF-β may show different patterns: monotonic increase, monotonic decrease, or alternation between growth and shrinkage.

Dynamic transition to N1-dominant TANs by normalizing the complex immune system
In Fig 9A, we show the sensitivity of G-dependent N2 states (C s ) to changes in inhibition strength of the N2 module by the N1 module (α = 0.5, 1.5, 2.2, 6.0). The corresponding N1 states (I s ) are shown in Fig 9B. As α is increased, the bifurcation curve (C s , G) shifts to the right. This increase in α also moves the bi-stability window (W G ) to the right and leads to a decrease in the size of the window (|W G |). This indicates that the probability of switching to the anti-tumorigenic phase (P a ) is increased when α is increased. For example, the immune response is already in the P a -phase (C < C th , I > I th ) when G = 0.6 in the case of the higher α value (α = 2.2; green curve in Fig 9A) while it is still in the P t phase (C > C th , I < I th ) in the base case (α � = 1.5; red line in Fig 9A). On the other hand, the immune system would initiate the tumor-favoring progress even under low TGF-β conditions when α is decreased. For example, the system would be in the P t -phase for a low TGF-β signal (G = 0.1) when α is lowered (α = 0.5; blue curve in Fig 9A) while it should be in the P a -phase in the base case (α � = 1.5). We recall that, for the base case (α � = 1.5), the existence of the bi-stability regime W G illustrates that the TGF-β signal can regulate the forward transition from the P a -state to the P t -state and the recurrent process from the P t -state to the P a -state. However, when the N2 activity is highly enhanced with the decreased inhibition signal from the N1 complex (a decrease in α), the TGF-β signal may regulate only forward transition from the P a -state to the P t -state (one way switch) and a monotonic decrease in TGF-β signals does not push the immune system to the anti-tumorigenic phase. The phenotypic changes between P a -and P t -states for increased or decreased α are illustrated in Fig 9C. This mechanism of low α's therefore provides a way of keeping the immune system in a tumorigenic phase. Functional plasticity associated with a decrease in α induces a critical immunosuppressive switch to promote aggressive tumor growth and cancer progression [35,74].
In Fig 10 we show the effect of the inhibition strength of N1 activities by the N2 module (β = 0.1, 0.7, 1.0 � , 6.0) on the dynamics of the core control system. As β is increased, the bifurcation curves (G, C) shift to the left (Fig 10A). In contrast to the case of changes in α in Fig 9, this increase in β reverses its shifting direction i.e., it moves the bi-stability window (W G ) to the left and leads to an increase in the size of the window (W G ) despite disappearance of the Role of TANs in cancer progression bistable switch for larger β. Thus, an increase in β induces transient changes from a bistable-to a one-way-switch and to a mono-stability. This indicates that the probability of switching to the tumorigenic phase (P t ) is increased in response to fluctuating TGF-β levels (0 � G � 1.0) when β is increased. For example, the immune system is already in the tumorigenic phase (P t ; C > C th , I < I th ) in response to G = 0.12 in the case of higher β's (for example, β = 6.0; blue curve in Fig 10A) while it should still be in the anti-tumorigenic phase (P a ; C < C th , I > I th ) in the base case (β � = 1.0; red line in Fig 10A). On the other hand, the immune system would activate the anti-tumorigenic machinery even under high TGF-β signal conditions when β is decreased. For example, the system would increase probability of switching to the anti-tumorigenic phase by lowering N2 levels for a high TGF-β level (G = 1.0) when β is lowered (β = 0.1; cyan curve in Fig 10A) while it should be in the tumorigenic phase in the base case (β � = 1.0). While a decrease in β reduces the size of the bi-stability window (|W G |), an increase in β reduces the size of the one-way switch window. When N2 activities are highly enhanced with the increased inhibition signal of the N1 module by the N2 complex (an increase in β) due to mutual antagonism between N2 and N1 modules, TGF-β signals may also regulate only forward transition from P a -state to P t -state (one way switch) and a monotone decrease in TGF-β signals does not push the immune system to the anti-tumorigenic phase. Therefore, this βinduced mechanism has a similar effect on keeping high N2 activities as that of a decrease in α, providing a way of keeping the immune system in a tumorigenic phase.
Therefore, these results provide a scientific basis for targeting suppressive effect of N2 TANs by normalizing immune activities (increase in α or decrease in β) within the tumor microenvironment or keeping a favorable balance of the N1/N2 TAN's properties toward N1 TANs [35]. For example, blood normalization was suggested to be an alternative way of creating anti-tumor TAN microenvironment [75,76].

Role of IFN-β signaling in regulation of the transition between N1-and N2-dominant system
In Fig 11, we investigate the role of the IFN-β signaling pathway in the phenotypic balance of N2 and N1 complexes for various TGF-β levels (G = 0.1, 0.5, 0.65). When G is low, the system induces a bi-stable system for low IFN-β levels and N1-dominant phenotypes in response to intermediate and high IFN-β levels ( Fig 11A). Therefore, the system is quite sensitive to relatively low IFN-β stimuli to induce the anti-tumor microenvironment. When G is increased to G = 0.5, the bistability window (W s ) in Fig 11A moves to the right and the system now can generate tumorigenic, bi-stable, and anti-tumorigenic status in response to low, intermediate, and high IFN-β levels, respectively. Therefore, the relative balance between TGF-β and IFN-β signaling determines the phenotypic preference (Fig 11B). On the other hand, in the presence of high levels of TGF-β (G = 0.65; Fig 11C), the bistability window disappears and the probability of inducing a N1-dominant microenvironment is decreased, shifting the balance toward the high N2-to-N1 ratio. Fig 11D summarizes the phenotypic representation of the system in response to various TGF-β and IFN-β stimuli. The majority of the bi-stable region resides in the lower left corner of the G − S plane. Therefore, the biochemical fluctuation of both TGF-β and IFN-β levels in the lower regime is likely to generate a selection process of either promoting or suppressing tumor growth based on the initial distribution and phenotypic status of TANs in tumor microenvironment.

Sensitivity analysis
In the model developed in this paper, there are parameters for which no experimental data are known and they may affect the simulation results. We take all parameters in the model (r, K, γ 1 , T 0 , λ, δ, G, S, k 1 , k 2 , k 3 , k 4 , α, β, and μ) for sensitivity analysis. We investigated the sensitivity of the activities of N2 and N1 modules at t = 10, 100, 200 to these parameters. We have chosen a range for each of these parameters and divided each range into 10,000 intervals of uniform length. For each of the fifteen parameters of interest, a partial rank correlation coefficient (PRCC) value is calculated. PRCC values range between -1 and 1 with the sign determining whether an increase in the parameter value will decrease (-) or increase (+) activities of the N2 and N1 complex, and tumor size at a given time. The sensitivity analysis to be described below was carried out using the method from [77] and Matlab files available from the website of Denise Kirschner's Lab: http://malthus.micro.med.umich.edu/lab/usadata/. Fig 12 shows the sensitivity analysis results for the N2 level and N1 activity, and tumor volume at t = 10, 100, 200 for Eqs (5)-(7) i.e., we compute PRCC values for C, I, T at time t = 10, 100, 200. We calculate PRCC values and associated p-value for fifteen perturbed parameters (r, K, γ 1 , T 0 , δ, λ, G, S, k 1 , k 2 , k 3 , k 4 , α, β, μ). We conclude that the activity of the N2 complex at t = 200 is positively correlated to the parameters λ, G, μ but is very little correlated (and thus not sensitive) to r, K, γ 1 , T 0 , α, β, k 1 , k 2 , k 3 , k 4 , S. Thus, in particular, the N2 level will increase significantly if the TGF-β signal is increased. Due to mutual antagonism between N1-and N2-dominant system, it is expected to see opposite correlations between N1 activity and those parameters. Indeed, the N1-dominant system is positively correlated to the parameters k 2 , S but is negatively correlated to the parameters G, β, k 4 , μ, i.e., there would be significant reduction in the N1 expression level if the TGF-β signaling strength is increased. The tumor volume For a constant IFN-β signaling, the system transits from P a -to P t -modes as the TGF-β level is increased. On the other hand, for a fixed TGF-β level, an increase in S induces a transition from P t -favoring mode to the P a -status. Other parameters are given in Table 2. https://doi.org/10.1371/journal.pone.0211041.g011

Role of TANs in cancer progression
is positively correlated to the growth rate (r) and carrying capacity (T 0 ) but is negatively correlated to the tumor killing rate by N1 TANs (δ).

Therapeutic approaches: Blocking the N2-dominant system by TGF-β inhibitors
In Fig 13 we investigate the anti-tumor effect of a TGF-β inhibitor on tumor growth. We introduce new variables: concentration of TGF-β (G), instead of a parameter, and TGF-β inhibitor concentration (L) whose dynamics are described as follows: Here, G s is the constant source of TGF-β, μ G is the decay rate of TGF-β, γ L is the degradation rate of TGF-β by its inhibitor, L s is the external periodic source of the TGF-β inhibitor Role of TANs in cancer progression inhibitor. Eqs (8) and (9) then are coupled with the original Eqs (5)-(7) to investigate the effect of TGF-β inhibitors on the transition behaviors between N1 and N2 phenotypes. In Fig 13(A) and 13(B), the dynamic system shows different mode changes as the TGF-β inhibitor was injected at the different periods (τ L = 1, 3). A frequent injection of the inhibitor (τ L = 1 in Fig  13A) induces a decrease in the TGF-β level (G < 0.4), resulting in the transition from the initial tumorigenic condition (black arrow in Fig 13C) to the anti-tumorigenic mode (blue arrowhead in Fig 13C). On the other hand, the TGF-β level still stays in the upper stable branch ( Fig  6) in response to a less frequent injection schedule (τ L = 3 in Fig 13B), keeping the system in the tumorigenic phase (red arrowhead in Fig 13C). These different cellular compositions in N1-or N2-dominant tumor microenvironment result in fast or slow tumor growth: slow and fast tumor growth in response to frequent (τ L = 1) and less frequent (τ L = 3) injection schedules, respectively (Fig 13D). Fig 13E shows Table 2. https://doi.org/10.1371/journal.pone.0211041.g013

Role of TANs in cancer progression
inhibitor needs to be injected into the tumor in order to maintain the same level of anti-tumor efficacy in the tumor microenvironment.

Therapeutic strategies by IFN-β injection
We investigate the effect of normalizing immune activities on tumor growth by injecting IFNβ. We consider the following conventional equations for time-dependent IFN-β administration:  Fig 14B), the N2-dominant system responds to the injected IFN-β in a periodic fashion and slowly transits to the N1-dominant system (pink box in Fig 14B) by decreasing the N2 level and increasing N1 activity (blue arrowhead). Fig 14C shows time courses of tumor volume (mm 2 ) in the absence (PBS; red solid curve, S s = 0) and presence (blue dashed curve, S s = 35) of intratumoral injection of IFN-β. IFN-β treatment significantly reduces the tumor size through repressing the N2-suppressing immune microenvironment (Fig 14B), which is in good agreement with experimental data. Note that even though the (C, I) returned back to the tumorigenic phase (blue arrowhead in Fig 14B) after the last IFN-β injection on day 9 (Fig 14A), the overall tumor growth is still very slow due to the earlier suppression of the N1-dominant environment, thus normalizing the immune activities against the tumor. The empty circles and triangles in Fig 14C represent Fig 14F. In experimental studies [28], the relative population of tumor-infiltrating CD8 T cells were significantly increased in tumor microenvironment in response to IFN-β treatment but the Treg population was reduced in non-small cell lung cancer. The model predicts a significant increase in the immune activities ( Fig 14E) and dramatic decrease in Treg populations, which is consistent with the experimental results [28].
In Fig 15A, we show anti-tumor efficacy for various dose schedules and injection doses of IFN-β. For a fixed IFN-β dose, the low (or high) dosing frequency induces N2-dominant (or N1-dominant) tumor microenvironment and thus faster (or slower) tumor growth. For a fixed injection frequency, the system transits from the N2-dominant phase to the N1-dominant phase, thus slower tumor growth, as the injection dose amount (S s ) is increased as shown in Fig 14. The mathematical model predicts that the larger injection period (τ s ) and smaller dose (S s ) would result in the N2-mediated, faster tumor growth while the frequent injection schedule and larger doses would lead to the phenotypic switch to the N1-dominant microenvironment and the greater anti-tumor efficacy. Unfortunately, IFN therapy has been shown to be associated with side effects and significant toxicity in cancer, including neuropsychiatric, constitutional, hepatic, and hematologic effects [78], as well as in multiple sclerosis [79]. Therefore, clinicians would want to optimize the overall administration schedule in order to minimize the dose amount while still avoiding high dose therapy and maintaining sustainable clinical results. In other words, the dose level needs to be minimized in order to avoid drug complication region in Fig 15B. On the other hand, too frequent injections are not desirable due to an  Table 2.  1 day). Here, Γ(t)(= I(t)/C(t)) and S(t) are the N1/N2 ratio and IFN-β concentration, respectively, at time t. An enlarged subpanel (black box on the left) in Fig 15C is shown in Fig 15D. As the injection strength is increased (10!35!100), the N1/N2 ratio is increased and the tumor microenvironment transits from N2-to the N1-dominant system (gray arrow), leading to better anti-tumor efficacy (Fig 15A). High doses of IFN-β were shown to suppress tumor cell proliferation and promote cell death in both in vitro and in vivo experiments [80][81][82]. However, a high IFN-β dose also increases IFNβ levels, leading to drug complications or extensive associated toxicity (gray zone on the top) As the injection strength is increased (10!35!100), the trajectory transits to the N1-dominant system (gray arrow). However, the high dose IFN-β may lead to drug complications or extensive associated toxicity (gray zone on the top) [83][84][85][86]. N1 and N2 phenotypes on the bottom of panels (C, D) are determined from the N1/N2 ratio. https://doi.org/10.1371/journal.pone.0211041.g015 Role of TANs in cancer progression [83][84][85][86]. Therefore, the injection schedule should be designed to avoid the associated toxicity but still maximize the N1/N2 ratio for optimal clinical outcome.

Combination therapy
In this section, we investigate the effect of a combined therapy (TGF-β-inhibitor +IFN-β) therapy on tumor growth by solving Eqs (5)- (10). Fig 16A shows time courses of the activities of N2 (C) and N1 (I) modules, and concentrations of TGF-β (G), TGF-β inhibitor (L), and IFN-β (S) in response to a periodic injection of TGF-β inhibitor every day (τ L = 1 day) and IFN-β on 0, 3, 5, 7, 9 days (L s = 8, S s = 20). Injections of both two drugs induce a decrease in C and increase in I, generating the N1-favoring microenvironment, thus suppressing tumor growth. When these drugs are administered more often (τ L = 1 day, τ s = 1 day), the system quickly switches from the N2-favoring system to the N1-dominant system even with lower doses (L s = 5, S s = 15; See Fig 16B). Fig 16C shows the trajectories of solutions (C(t), I(t)) in response to control (PBS, red) and injections of IFN-β alone (blue), TGF-β inhibitor alone (TGF-β − , green), and combination therapy (IFN-β + TGF-β inhibitor, black): L s = 8, S s = 20. While solutions in control (PBS) and single drug (either IFN-β or TGF-β inhibitor) stay in the tumorigenic phase (arrowheads (red, blue, green) in the blue box), the solution in the combination therapy (TGF-β inhibitor+IFN-β) moves to the anti-tumorigenic region (black arrowhead in the pink box). Thus, these dynamics illustrate how the combination therapy may induce the transition from the P t -phase to the P a -phase even when either TGF-β inhibitor or IFN-β  Table 2. https://doi.org/10.1371/journal.pone.0211041.g016 Role of TANs in cancer progression therapy transits the system to the P t -mode. Time courses of the tumor volume corresponding to four cases in Fig 16C are shown in Fig 16D. Thus, the combination therapy leads to synergistic tumor killing effects via switching to the P a -phase (Fig 16C). Fig 16E shows anti-tumor efficacy of the combination therapy with various doses of IFN-β and TGF-β inhibitor at final time (t = 25 days). For a fixed dose level of either IFN-β or TGF-β inhibitor, the anti-tumor efficacy is increased in general. However, high doses of either IFN-β [78,[85][86][87] or TGF-β inhibitor [88,89] alone can induce toxicity, resistance and side effects. Therefore, careful dosing strategies need to be developed in order to efficiently deliver IFN-β to a tumor and maximize anti-tumor efficacy while minimizing critical, systemic side effects [90,91]. Therefore, double high doses of those two drugs may not be desirable due to economic costs as well as serious side effects. It has been shown that high doses of TGF-β inhibitors may lead to poor response to the combined 5-fluorouracil and IFNα-2b therapy in hepatocellular carcinoma patients [92]. The optimal drug dose schedule (red asterisk) can be obtained by avoiding drug complications (dotted region on the upper and right sides) and minimizing dose levels of both drugs while maintaining the overall high anti-tumor efficacy. Fig 16F shows the trajectories of solutions (L(t), S(t), Γ(t)) for the low (blue; (I) in Fig 16E), intermediate (red; (II) in Fig 16E), and high (green; (III) in Fig 16E) doses of both TGF-β inhibitor and IFN-β. Here, Γ(t) = I(t)/C (t) is the N1-to-N2 ratio. While the N1-to-N2 ratios in (II) and (III) are large, the N1-to-N2 ratio is small for low doses of IFN-β and TGF-β inhibitor (case (I)). Thus, the case (III) in Fig  16E is not desirable because of the high N1-to-N2 ratio.

Effect of time delays
We investigated the effect of time delays on regulation of phenotypic switches and tumor growth. We prescribe time delays in the suppression terms of the N2 and N1 complex. There are many cellular and molecular factors that may affect delayed processes of signaling pathways, apoptosis, and inhibition of neutrophils in the microenvironment [62,[93][94][95][96]. This may also be caused by partial blocking of signaling pathways of these TANs by drugs. Recall that the transition probability between P t and P a -phases can be modified by the mutual inhibition strength of both N1 and N2 complex (α, β; Figs 9 and 10). Therefore, we modify Eqs (5) and (6) by introducing time delays in the inhibition terms as follows where Δ 1 , Δ 2 are the time delays in the inhibitory pathway of the N2 module and suppression of N1 module, respectively (See Fig 4B). Fig 17A shows  Starting with the same initial condition (black arrow), the system adapts to changes in the microenvironment by inducing tumorigenic status (C high, I low; blue box) and anti-tumorigenic mode (C low, I high; pink box) in the absence and presence of time delays, respectively. This delay-induced P a mode results in a decrease in the tumor size (Fig 17C). This illustrates that the time delays in suppressing N1 or N2 TANs may stimulate a homeostatic imbalance, resulting in either promotion or suppression of tumor growth. Fig 17D shows the normalized tumor volume at day 25 for various strengths of the two time delays (Δ 1 , Δ 2 ). For a fixed Δ 1 , as Δ 2 is increased, the inhibition strength of the N1 module is decreased and N1 TANs become a dominant phenotype, decreasing the tumor volume. For a fixed Δ 2 , the tumor size is a nonlinear function of Δ 1 , inducing a small tumor for either high or low Δ 1 but assuming the maximum tumor size for an intermediate Δ 1 . These results illustrate the complex dynamics of the N1/N2 transitions and nonlinear tumor growth in the presence of time delays in mutual inhibitory properties. Now, we investigate the effect of time delays in apoptosis of TANs by modifying Eqs (5) and (6) by introducing time delays in the decay term as follows where Δ 3 , Δ 4 are the time delays in apoptosis of N2 and N1 modules, respectively. Role of TANs in cancer progression position was marked in black arrow. While the presence of time delays (Δ 3 , Δ 4 ) = (1, 1) induces initial fluctuations (Fig 18A), the solutions converge to the N2-dominant steady states of ODEs (red arrowhead in Fig 18C). However, a slightly different combination of time delays ((Δ 3 , Δ 4 ) = (0.9, 1)) leads to a phenotypic transition from the N2-dominant mode to the N1-dominant phase (blue arrowhead in Fig 18C) despite the initial fluctuation of solutions (blue curves in Fig 18B). Delayed apoptotic cell death leads to different tumor growth dynamics. Fig 18D shows time courses of the tumor volumes corresponding to three cases in Fig 18C. In the presence of time delays (Δ 3 , Δ 4 ) = (1, 1), the tumor grows slower than the control (case without time delays) initially but the delay-induced oscillations push the TAN system to the N2-dominant phenotype faster than the control case around t = 18 days, leading to the larger N2-to-N1 ratio and faster tumor growth (red curve in Fig 18D). On the other hand, in the presence of time delays (Δ 3 , Δ 4 ) = (0.9, 1), the tumor grows slower than the control due to strong suppression of the N2 module and enhancement of the N1 module. Fig 18E shows the normalized tumor size at t = 25 days for various strengths of the two time delays (Δ 3 , Δ 4 2 [0, 1]) with a uniform length of 0.1. A risk map of time delay-induced transition of TANs from N2 to N1 phenotypes in the Δ 3 − Δ 4 plane is shown in Fig 18F. For a fixed value of Δ 3 , an increase in Δ 4 results in the higher chance of transition to the P a -phase, thus slower tumor growth. On the other hand, for a fixed value of Δ 4 , an increase in Δ 3 leads to a higher probability of transition to the P t -phase and faster tumor growth. To characterize the complex dynamics of the apoptotic cell death mechanism, we introduce the following Δ 4 -to-Δ 3 ratio: time delay in apoptosis of N1 TANs time delay in apoptosis of N2 TANs : When D ap is large, the delayed apoptotic process of N1 phenotypes and minimal delay in the apoptotic process of N2 phenotypes allow accumulation of more N1 TANs within the Role of TANs in cancer progression microenvironment, pushing the system to the P a -phase and effective inhibition of tumor growth. On the other hand, a large delay in the apoptosis of N2 TANs with a relatively small delay in the apoptosis of N1 TANs, (i.e., D ap � 1), leads to significant accumulation of N2 TANs in the tumor microenvironment, leading to faster tumor growth. Fig 19A shows the tumor volume for various D ap at final time t = 25 day. As D ap is increased from the control value (D ap = 1), the tumor size is significantly decreased. On the other hand, a decrease in D ap also induces slower tumor growth, indicating the maximum tumor growth rate when D ap = 1. The absence of the time delays in apoptosis signaling also induces slower tumor growth. Fig 19B illustrates a dynamic antitumor killing model with the relative time delay function D ap as a key parameter. A relatively small delay ratio (D ap � 1) (basic step i) induces phenotypic perturbations in N1-or N2-responses, which in turn decreases the N2-assisted tumor growth rate. On the other hand, the depletion of time delays in the apoptosis pathways of TANs (step ii) improved antitumor efficacy and suggests that the major antitumor response was likely due to the persistent N2-dominant microenvironment with quick clearance of TANs (Fig 18A). Our study also showed that complementing TANs with adjuvant N1 immunotherapy improved overall antitumor efficacy (step iii), through the facilitation of a robust N1-mediated tumor cell killing response beyond the limits of endogenous N2 TANs. In spite of a typical short half life of neutrophils, delayed apoptosis and accumulation of neutrophils in a tumor environment have been observed in other biological studies [93,95,96]. In particular, TANs often have a longer lifespan in the tumor microenvironment due to delayed apoptosis of TANs by modifying its signaling pathways involving FAS, ROS, active Caspase 3, 9, and Bcl-2 family [62]. Our mathematical model predicts that (i) D ap can be used as a reliable diagnostic measure for aggressiveness of tumor progression; (ii) Anti-tumor efficacy can be improved by increasing the apoptosis decay of N1 TANs or by blocking the accumulation of N2 TANs in tumor microenvironment. In particular, drugs targeting apoptosis pathways of N1 TANs can effectively be used to enhance N1 neutrophil activities to reduce the tumor size.

Discussion
The immune system protects the host against various types of biological threats. However, it is well established that immune cells show functional plasticity in the context of cancer, and undergo phenotypic changes, promoting tumor growth and metastatic progression [74]. There have been substantial studies to pinpoint fundamental mechanisms through which TANs act on cancer progression [8, 17,34], and especially via the coordination of chemotaxisdriven recruitment and activation of distinct immune cells to the tumor microenvironment [97][98][99][100] as well as immune suppression [74,101,102]. A groundbreaking work by Fridlender and colleagues [16] had illustrated that TGF-β in tumor microenvironment can generate a functional transition from the proinflammatory neutrophil phenotype (termed "N1") to the anti-inflammatory phenotype (termed "N2"); [74]. IL-6 also contributes to the formation of the N2 microenvironment [103,104]. Neutrophil-to-lymphocyte ratio (NLR) was highly associated with progression in many cancer types [69] and has been suggested to be a prognostic factor (or biomarker) for various cancers including colorectal cancer [105,106], nasopharyngeal carcinoma [107], non-small-cell lung cancer [71], breast cancer [108,109], hepatocellular carcinoma [110], and melanoma [13]. Recently, it was shown that lung adenocarcinomas can promote bone stromal activity and increase bone mass even without bone metastasis, which in turn enhances tumor growth by remotely supplying tumor-infiltrating neutrophils [111]. However, how the systemic host environment can communicate with a tumor at a remote site is poorly understood.
Since TGF-β was identified as a master player that skews differentiation toward the N2 phenotype [16,18,19], TGF-β inhibitors are suggested to shift the TAN balance toward the N1 phenotype [20,21]. In cancer research, type I IFNs (α and β) are substantially studied and tried for cancer therapy due to their anti-tumor capabilities by acting directly on cancer cells and through immunoregulatory characteristic [112]. These Type I IFNs can induce cell cycle arrest and apoptosis in several models [112][113][114] by acting on several signaling pathways including upregulation of p53 transcription [24]. Type I IFNs are suggested to mediate initial recruitment, active proliferation, phenotypic differentiation, and activation of various immune cells [115][116][117][118]. Interestingly, recent evidence from animal models has implicated that IFN-β signaling pathways play a critical role in inducing antitumor response in chemotherapy [119] and radiotherapy [68,120]. Therefore, it is not surprising to see IFN-β gene transfer in extensive studies, due to the ability of IFN-β to improve (or normalize) immunological response in the tumor microenvironment [121][122][123][124][125]. IFN-mediated neutrophils was shown to upregulate PD-L1 and suppress T-cell proliferation [26].
We have developed a mathematical model by a system of ODEs and DDEs to study the effect of regulatory cytokines (TGF-β, IFN-β) and delayed apoptotic cell death mechanism on transitions between N1-and N2-dominant system, the N1/N2 mutual antagonism, and their impact on tumor growth (Fig 3). Mathematical analysis showed the complex, nonlinear behaviors of phenotypic switches between N1-and N2-dominant phases in response to fluctuating TGF-β and IFN-β, and tumor growth patterns (Figs 6-11). In particular, the bi-stability of the N1/N2 system generates a selection process of either promoting or suppressing tumor growth based on the initial distribution and phenotypic status of TANs in a tumor microenvironment (Figs 6 and 11). The microenvironmental pressure of inhibiting the phenotypes also determines the phenotypic switches (Figs 9-10). We also developed therapeutic strategies to control this complex network of TANs for slowing down tumor growth (Figs 13-16). Even though high-dose IFNs were suggested to the only available therapy for relapse-free benefits and overall survival (OS) [126], the administration of high-dose IFN in cancer patients in stage IIb/III needs to be considered carefully [112]. High-dose IFN is not always considered as the first line of treatment option for melanoma patients due to a lack of substantial evidence supporting a long-term OS merit and the extensive associated toxicity [83][84][85][86]. Therefore, it is necessary to establish optimal dose and treatment duration of IFNs despite extensive clinical experience [112]. In this study, we investigated the optimal dose schedule of IFN-β in the absence or presence of an alternative drug, the TGF-β inhibitor (Fig 15). In our mathematical framework, the high N2-to-N1 ratio (equivalent to NLR) was highly associated with faster tumor growth and worse clinical outcomes (Fig 15). Experimental studies [22-25, 36, 38] are in good agreement with the conclusions of the mathematical model. The optimal strategies of reducing the tumor size with minimal side effects [78,[85][86][87][88][89] and maximal antitumor efficacy were explored by using the TGF-β inhibitor only, IFN-β only, and the combination (TGF-β inhibitor + IFN-β) therapy (Fig 16).
The life span of neutrophils is strictly regulated in order to maintain tissue homeostasis due to their potential toxicity [127]. These cells are quickly removed from circulation after leaving the bone marrow with a very short half life. However, several proinflammatory cytokines were reported to influence their longevity [128,129]. There are also many cellular and molecular factors that may affect delays in signaling pathways, apoptosis, and inhibition of neutrophils in the given microenvironment [93][94][95]. For example, external (smoking and nicotine) and internal (ROS) factors were shown to delay neutrophil apoptosis by suppressing signaling pathways (InsP7, Akt, ROS) in lung diseases such as COPD and lung cancers [93,95,96]. Recently, it was also shown that ROS production is reduced in the absence of endogenous IFN-β, causing the delayed apoptosis of TANs [62]. Interestingly, a recent study showed that a radiation therapy (RT) can induce the polarization of N1 TANs and reactive oxygen species induced by RT damage tumor tissues, which overall improves anti-tumor efficacy [130]. These microenvironmental factors delay the inhibition and apoptosis of both types of TANs, inducing a complex imbalance between N1 and N2 phenotypes. The mathematical model predicted that various combinations of two time delays in inhibition of signaling pathways within those TANs will induce the relative imbalance between N1 and N2-dominant phases, leading to promotion or suppression of aggressive tumor growth (Fig 17).
We have shown that typical therapy can be exploited by manipulating the delay of apoptosis pathways of both types of TANs. By reducing the relative ratio of apoptosis delay of N1 and N2 TANs, or by increasing time delay of apoptotic pathways of N1 TANs, the efficacy of the antitumor treatment is increased. The mathematical model also predicts that the absence of time delays in apoptosis may allow a greater chance to induce N1 TANs in response to intermediate levels of TGF-β, thereby increasing antitumor efficacy (Fig 18). In particular, the model predicted that either strategy, i.e., induction of small (or none) delay of N2 TANs or large delay of N1 TANs, increased the survival rate (Fig 19), even though it is not clear which of the two strategies is clinically better. The answer should also depend on negative side effects associated with the delays of TANs (mixed N1/N2 immune response) and with the injection of apoptosistargeting drugs (toxicity and inflammation). Our results provide a scientific basis for targeting suppressive effect of N2 TANs by normalizing immune activities (increase in α or decrease in β) within the tumor microenvironment or keeping a favorable balance of the N1/N2 TAN's properties toward N1 TANs [35]. While blood normalization was suggested to be an alternative concept of creating an anti-tumor microenvironment including N1 TANs [75,76], more innovative platforms of normalizing the immune system need to be invented.
In this work, we did not specifically consider many other microenvironmental players such as signaling networks [131], human neutrophil elastase as a therapeutic target in the absence [132,133] and presence [134] of LPS, endogenous NK immune dynamics [40,135], angiogenesis from blood vessels [136][137][138], ECM remodeling [139], or growth factors such as epidermal growth factor (EGF) [140,141], and CSF-1 [7,142]. These factors may play critical roles in cancer progression. For example, neutrophil elastin (NE) was shown to indirectly promote the tumor growth by stimulating the PI3K signaling pathways in lung cancer cells [143]. Neutrophil extracellular traps (NETs) were abundant in tumour sections and induce the transformation of B cells [144], contributing to cancer progression [145]. However, another study [146] suggested that the NETs suppressed tumor growth in colonic adenocarcinoma. It is possible that unforeseen microenvironmental factors could have an effect on our approach (governing equations and the many estimated values for mathematical and biological variables). For example, in our model, it is assumed that the tumor microenvironment neutrophils can switch freely from N1 to N2 TANs or from N2 to N1 TANs, given a specific biochemical stimulus. However in vivo there may be unknown biochemical limits to these transitions, such as factors that stabilize one state or another after it has been established. This could be the presence of NETs and associated pro-inflammatory factors. The exact role of this biochemical structural device in regulation of promotion [147][148][149][150] or suppression of tumor progression [146] still remains controversial. Our model, however, provides a starting point of the general framework of those key transitions between N1 and N2 TANs in response to known key players. We also plan to investigate the role of the possibly continuous spectrum of the N1 ! N2 transition. Further investigation and experimental validation need to be done as more experimental data are available. A multi-scale mathematical model [64,[151][152][153][154][155][156] could be used to take into account inter and intra-cellular signaling regulation at the microscale level and integrating those to stromal cells and tumor cells at the cellular level. However, the mathematical model in this paper is a first step toward further understanding of lung cancer development in a tumor microenvironment and further experimental investigation. We hope to address these issues in future work.