Tumor cell density (= n(x, t))

The mass balance equation for the tumor cell density n(x, t) is

We assume that the tumor extracellular matrix is homogeneous and isotropic in tumor microenvironment, and that the flux due to the random motility is given by

In lung tissue, tumor cells are strongly attracted to chemotactic attractants [32] such as NE and neutrophils [22, 32] and migrate toward the up-gradient (∇E) of the chemo-attractant, NE, through the process called ‘chemotaxis’ [50]. The chemotactic flux is assumed to be of the form

Tumor cell invasiveness is enhanced by proteolytic degradation of the extracellular matrix via MMPs [4, 25, 53] and NEs [29, 54] that are produced by neutrophils. This results in local degradation of tumor ECM [49] and tumor cell movement in the direction of the up-gradient (∇ρ) of ECM via a cellular process called haptotaxis. This process is valid only in the ECM domain S, therefore, we include the characteristic function I_S, providing the on-off switch on the ECM membrane. We represent the haptotactic flux in a similar fashion:

The net production of tumor cells is due to active NE-stimulated growth [8, 22] and cell killing by N1 TANs [3, 21, 23], which we represent as follows:

Combining the several fluxes in Eqs (2)–(4) and growth term in Eq (5) leads to the governing equation for the tumor cell density

Densities of neutrophils: N1 (= N₁(x, t)) & N2 (= N₂(x, t)) types

We use a similar form of reaction-diffusion-advection equations for the evolution of the densities of neutrophils, based on mass balance as in the previous section. We assume that (i) Neutrophils are chemotactic to the CXCL secreted by tumor cells [55–57], and the chemotactic flux is of the nonlinear form (3), but with different chemotactic sensitivities (χ₁, χ₂). Since the N2 TANs produce NE and MMPs, the movement of activated neutrophils further enhances tumor invasiveness and growth via the NE-PI3K pathway described earlier. (ii) The anti-tumorigenic (N1) neutrophils transform into the active N2 type at the rate λ₁₂ in the presence of TGF-β, based on experimental evidences [3, 21, 43]. For instance, the N1→N2 transition of TANs with protumour properties was typically observed in a TGF-β-rich tumor microenvironment and the presence of IFN-β or TGF-β inhibitor can mediate the reverse transition (N2→N1) with anti-tumoral properties [3]. Therefore, TGF-β pathway inhibitors are under clinical trials since they were shown to promote the development of N1 TANs [58, 59]. (iii) N1 and N2 phenotypes proliferate at a rate, λ₁ and λ₂(G), respectively. Then, we have the following evolution equations:

Tumor ECM density (= ρ(x, t))

The tumor ECM provides structural foundation for efficient cell migration [60], but it also needs to be remodeled via proteolysis for tumor cell migration by microenvironmental proteases [55, 61–63]. In this work, we assume that the tumor ECM is degraded by the TAN-secreted NEs [64] and TAN-secreted MMPs [4, 25, 53, 55] as in the invasion experiments [22]. The rate of ECM change can be represented as

CXCL8 concentration (= C(x, t))

Tumor cells secrete CXCL in order to recruit the immune cells such as neutrophils [55–57]. CXCLs and corresponding receptors (CXCR) such as CXCL5 and CXCR6 are important prognostic factors, alone or in a combination with the TANs, for shorter overall survival and cumulative risk of recurrence [3, 65, 66]. Thus the governing equation for CXCL8 is

TGF-β concentration (= G(x, t))

TGF-β is a polypeptide that plays a major role in regulation of many human diseases including cancers [67] due to its capacity of maintaining tissue homeostasis and involving in most of the chronic inflammatory and wounding processes by activating its inactive form in ECM [68]. Tumor cells are the primary source of TGF-β in TME [3, 69]. TGF-β activates proinflammatory and antitumorigenic N1 neutrophils into the aggressive N2 type, which in turn stimulates tumor cell invasion [55, 70]. Thus the governing equation for TGF-β is as follows:

Concentrations of NET/NE (= E(x, t)) and its inhibitors (= D(x, t))

NET and NE are highly associated with aggressive invasion, growth, EMT, and metastasis of cancer cells [4, 54, 71]. NE is produced by neutrophils [4, 64, 72] and used for degradation of extracellular matrix and tissue destruction [8, 54, 73]. It was also shown that NE inhibitors such as DNase I block this effect in growth models [8] and invasion assays [22]. In our framework, NET and the associated NEs are merged into one component. Thus, the governing equations for NET/NE and its inhibitors are

MMP concentration (= P(x, t))

Matrix metalloproteinases (MMPs) are highly associated with cancer cell invasion and metastasis [48, 74]. Neutrophils, not tumor cells [55, 75], were suggested to the primary source of MMPs [4, 25, 53, 55] including MMP-9 [53] in lung cancer development, showing strikingly predominant presence at the invasive fronts of metastatic cancers [53]. Thus the governing equation for MMPs is

TIMP concentration (= M(x, t))

Tissue inhibitors of metalloproteinases (TIMPs) play an important role in inhibiting tumor invasion and metastasis [77] by regulating major signaling pathways in pericellular proteolysis of various extracellular matrix and cell surface proteins [78]. In the model, TIMPs are injected for inhibition of the proteolytic activities of cancer cell invasion. Note, however, that this action can partially block cancer cell invasion since cancer cells can still execute the NE-mediated invasion. Thus, the governing equation of TIMP is

Boundary conditions and initial conditions

In the following simulations we prescribe Neumann boundary conditions on the exterior boundary Γ₁ (= ∂Ω; see Fig 2B) as follows:

Finally, we prescribe initial conditions,

Parameters are given in Tables 1 and 2. Nondimensionalization and parameter estimation of the system (6)–(19) are given in S2 and S3 Text, respectively. This non-dimensional form of governing equations was used for the simulations. Hereafter, the computational domain is restricted to one space dimension, and the computational domain is scaled to unit length.

Table 1

Parameters used in the tumor model.

	Description	Dimensional Value	Refs.
Diffusion coefficients (cm2 s−1)
Dn	Tumor cells	2.5 × 10−8	[33, 37, 129–131]
D1	N1 Neutrophil	1.1 × 10−8	[83]
D2	N2 Neutrophil	1.1 × 10−8	[83]
DC	CXCL8 (IL-8) chemokines	2.5 × 10−6	[57, 82, 83]
DG	TGF-β	1.0 × 10−6	[84–86]
DE	NET/NE	5.0 × 10−7	[132–137]
DP	MMP	5.0 × 10−10	[37, 138, 139]
DD	DNase I (NET/NE inhibitor)	7.374 × 10−6	[140], estimated
DM	TIMP (MMP inhibitor)	8.33 × 10−7	estimated
DA	TGF-β Anti-body	8.33 × 10−7	estimated
Production rates
r	Proliferation rate of tumor cells	3.3 × 10−4 s−1	[22, 33], estimated
rE	NE-mediated proliferation rate of tumor cells	7.0 × 10−2	estimated
kE	Hill type coefficient of tumor cell proliferation	2.15 × 10−9 gcm−3	[33], estimated
m	Hill coefficient of NE-mediated tumor proliferation	2	estimated
n0	Tumor cell carrying capacity	2.5 × 104 cells/cm3	[33], estimated
λ1	Proliferation rate of N1 neutrophil	4.38 × 10−6 s−1	[33, 141], estimated
λ12	Transformation rate from N1 to N2 neutrophils	4.08 × 103 cm3 g−1 s−1	[33], estimated
λ2	Proliferation rate of N2 neutrophils	2.65 × 10−5 s−1	[33, 141], estimated
λC	Production rate of CXC from tumor cells	4.44 × 10−11 s−1	[33, 142, 143], estimated
λG	Production rate of TGF-β from tumor cells	4.89 × 10−7 s−1	[33, 142, 143]
λE	Production rate of NET/NEs from Neutrophils	2.26 × 10−7 s−1	[144]
λP	Production rate of MMPs from Neutrophils	2.22 × 10−8 s−1	[144]
λD	Production rate of NET/NE inhibitor	9.0 × 10−13 gcm−3 s−1	estimated
λM	Production rate of TIMP	1.29 × 10−11 gcm−3 s−1	estimated
λA	Production rate of TGF-β anti-body	4.78 × 10−5 μMs−1	estimated
Degradation/Decay rates
μn	tumor cells degradation rate by N1	2.78 × 10−1 cm3 g−1 s−1	estimated
μρ1	ECM degradation rate by NEs	1.02 × 104 cm3 g−1 s−1	estimated
μρ2	ECM degradation rate by MMPs	3.19 × 105 cm3 g−1 s−1	estimated
μC	decay rate of CXCL8	6.42 × 10−5 s−1	[145–147]
μG	decay rate of TGF-β	8.02 × 10−6 s−1	[33], estimated
μE	decay rate of neutrophil elastase	8.02 × 10−6 s−1	[148]
μP	decay rate of MMP	5.0 × 10−5 s−1	[49]
μD	decay rate of NE inhibitor	9.627 × 10−5 s−1	[149–153]
μM	decay rate of TIMP	4.56 × 10−6 s−1	[154]
μED	degradation rate of NET/NE by its inhibitor	(2.8 × 10−4—2.8 × 10−2) s−1	estimated
μA	decay rate of TGF-β anti-body	6.42 × 10−5 s−1	[43, 155]
μAG	decay rate of TGF-β by anti-body	4.8 × 10−3 μM−1 s−1	[43], estimated
KD	Hill coefficient of NET/NE degradation by its inhibitor	3.2 × 10−9 g/cm3	[156], estimated
l	Hill coefficient of NET/NE degradation by its inhibitor	2	[156], estimated
μPM	degradation rate of MMPs by TIMP	(2.8 × 10−5—2.8 × 10−4) s−1	estimated

Table 2

Continued from Table 1.: Dimensionless values were marked in *.

Parameters used in the tumor model.

	Description	Dimensional Value	Refs.
Degradation/Decay rates
KM	Hill coefficient of MMP degradation by TIMP	4.64 × 10−8 g/cm3	[156], estimated
m	Hill coefficient of MMP degradation by TIMP	2	estimated
Movement parameters (chemotaxis, haptotaxis)
χE	Chemotactic sensitivity to NE	1.117 × 10−9 cm.s−1	[22, 33, 130, 157]
δE	Scaling parameter of NE gradient (|E|)	6.46 × 10−8 g/cm4	estimated
σE	Scaling parameter of NE gradient (|E|)	1.0	estimated
χρ	Haptotactic sensitivity	3.5 × 10−10 cm.s−1	[22, 158–160], estimated
δρ	Scaling parameter of ECM gradient (|ρ|)	5.0 × 10−3 g/cm4	estimated
σρ	Scaling parameter of ECM gradient (|ρ|)	1.0	estimated
$$	Chemotactic sensitivity of N1 neutrophils to CXCL8	1.1 × 10−9 cm.s−1	[22, 33, 161], estimated
δ1	Scaling parameter of CXCL gradient (|C|)	1.0 × 10−11 g/cm4	estimated
$$	Chemotactic sensitivity of N2 neutrophils to CXCL8	1.1 × 10−9 cm.s−1	[22, 33, 161], estimated
δ2	Scaling parameter of CXCL gradient (|C|)	1.0 × 10−11 g/cm4	estimated
σC	Scaling parameter of CXCL8 gradient (|C|)	1.0*	estimated
Membrane parameters
γc	permeability of cells (Tumor cells, N1 neutrophils, N2 neutrophils)	0.01-785 (78.5)*	estimated
γ	permeability of chemicals (CXCL8, TGF-β, NET/NE, NE inhibitor, MMP, TIMP)	0.1-7850 (785)*	[33], estimated
μ	ECM width	0.6*	estimated

All the simulations were performed using a finite volume method (FVM; clawpack (http://www.amath.washington.edu/~claw/)) with fractional step method [87]. A nonlinear solver nksol was used for solving algebraic systems. The Eqs (6)–(19) were solved on a regular uniform grid with grid size 0.01 (h_x = 0.01). An initial time step of 0.0001 (or smaller) was used, but adaptive time stepping scheme based on the number of iterations can increase or decrease this step size.

Predictions of the mathematical model

Fig 3 shows the density profiles of all variables (n, N₁, N₂, ρ, C, G, E, P) at t = 0, 5, 14, 22h in the absence of DNase I and TIMP when neutrophils were added in the lower chamber. In each subframe, the right (or left) half of the computational domain represents the upper (or lower) chamber in the Boyden invasion assay (Fig 2A). By degradation of the tumor ECM on the membrane of the insert, tumor cells in the upper chamber were experimentally shown to have capacity of invading the lower chamber upon stimulus of N1/N2 neutrophils in the lower chamber [22, 31]. Tumor cells in the upper chamber secrete CXCL8 (Fig 3F), which then diffuses and attracts neutrophils in the lower chamber by chemotaxis. On the other hand, tumor cells produce TGF-β (Fig 3B), which diffuses and enhances the N1 → N2 transformation of neutrophils (Fig 3C) in the lower chamber. These activated N2 TANs in the lower chamber then secrete NET/NEs (Fig 3D) and MMPs (Fig 3E) to stimulate chemotactic and haptotactic movement of tumor cells in the upper chamber. Tumor cells break down the ECM component by proteolytic activities with the NE and MMP near the membrane and invade the left chamber (Fig 3A). As they invade, they can sense higher levels of NE, and proliferate at a higher rate (see Eq (6)).

Fig 3

Dynamics of the system.

The time evolution of the density of each variable. (A) tumor cells and ECM (B) TGF-β (C) N1/N2 neutrophils (D) neutrophil elastase (E) MMP (F) CXCL8. Here, ECM = [0.35, 0.65]⊂ Ω = [0, 1]. Note that the initial concentrations of CXCL8, TGF-β, neutrophil elastase and MMPs are uniformly zero, as in experiments. *x-axis = space (the dimensionless length across the invasion chamber), y-axis = the dimensionless density/concentration of the variables.

A comparison of computational results from the mathematical model with experimental data [22] is shown in Figs 4 and 5. Hereafter, in order to calculate the population of cells (tumor cells, N1 TANs, N2 TANs) and level of chemical variables (CXCL-8, TGF-β, NET/NE, DNase, MMPs) at various times in the mathematical model, we integrate the density and concentration over the space: density of tumor cells ($$ \widehat{n}\left(t\right)={\int }_{\Omega }n(x,t)dx$$), N1 TANs ($$ {\widehat{N}}_1\left(t\right)={\int }_{\Omega }{N}_1(x,t)dx$$), N2 TANs ($$ {\widehat{N}}_2\left(t\right)={\int }_{\Omega }{N}_2(x,t)dx$$), and concentrations of ECM ($$ \widehat{\rho }\left(t\right)={\int }_S\rho (x,t)dx$$), CXCL-8 ($$ \widehat{C}\left(t\right)={\int }_{\Omega }C(x,t)dx$$), TGF-β ($$ \widehat{G}\left(t\right)={\int }_{\Omega }G(x,t)dx$$), NET/NE ($$ \widehat{E}\left(t\right)={\int }_{\Omega }E(x,t)dx$$), DNase I ($$ \widehat{D}\left(t\right)={\int }_{\Omega }D(x,t)dx$$), and MMPs ($$ \widehat{P}\left(t\right)={\int }_{\Omega }P(x,t)dx$$). In the experiments, Park et al. [22] found that the presence of neutrophils in the lower chamber could enhance tumor cell invasion through NET formation and NE activities, and the DNase treatment abrogated the invasion-promoting effect of neutrophils in the lower chamber. Fig 4A–4C shows time courses of the tumor population, population of invasive tumor cells, and neutrophil population (N1 (red solid), N2 (blue dashed) TANs), respectively, in the absence (control) and presence (⁺TAN) of neutrophils. In the presence of neutrophils, both total (Fig 4A) and invasive (Fig 4B) tumor cell populations are increased relative to the control case due to neutrophil transition (Fig 4C) and NE activities (Fig 3D) in the system. After 22h the number of 4T1 tumor cells invading the lower chamber almost doubled (∼190%) in the co-culture with neutrophils (red bar (⁺TAN); left panel in Fig 4D) in the lower chamber as compared to the control (blue bar; left panel in Fig 4D) in experiments [22]. In the model simulations, the number of invading tumor cells increased (∼2-fold) in the presence of neutrophils in the lower chamber (red bar (⁺TAN); right panel in Fig 4D) relative to the control (absence of neutrophil (blue bar); right panel in Fig 4D). As Park et al. [22] note, several tumor cell lines (4T1, BT-549) invade the lower chamber even in the absence of neutrophils in the lower well, which indicates the intrinsic invasiveness of tumor cells.

Fig 4

TAN-promoted cancer cell invasion (Experiment [22] & simulation).

(A-B) Time courses of populations of total tumor cells (A) and invasive tumor cells (B). (C) Time courses of N1 (red solid) and N2 (blue dashed) neutrophils. (D) Experimental data from the invasion assay in [22] (left column; 4T1 cancer cells) and computational results from mathematical model (right column). The graph shows the (scaled) populations of invasive tumor cells at t = 22 h in the absence (control) or presence (⁺TAN) of neutrophils. Here and hereafter cell populations are derived from the continuum density.

Fig 5

DNase I treatment against NE can abrogate the invasion-boosting effects of neutrophils (Experimental data [22] and simulation results).

(A) NE levels in the system in the absence (control) and presence (⁺TAN) of neutrophils, and DNase treatment (⁺TAN⁺D) cases at t = 22 h. (B) Experimental data from the invasion assay in [22] (left column; 4T1 cancer cells) and computational results from mathematical model (right column). The graph shows the (scaled; %) population of invasive tumor cells at t = 22 h in the absence (control; blue) and presence (⁺TAN; red shaded) of neutrophils, and DNase treatment (⁺TAN⁺D; yellow dotted) cases. Addition of DNase I reduces the number of invading tumor cells by almost 50%.

In Fig 5, we investigate the effect of DNase I against NET/NE on tumor invasion. Our mathematical model predicts that injection of DNase I in the lower chamber of the transwell can inhibit NET/NE activities (Fig 5A) and reduce the TAN-induced invasiveness of tumor cells (right panel in Fig 5B). Park et al. [22] showed that NE inhibition or digestion of the DNA of the NETs by DNase I can effectively abrogated the invasion-promoting effect of TANs in the lower chamber, i.e., the number of invasive 4T1 tumor cells was reduced in the presence of the neutralizing DNase I (⁺TAN⁺D) when compared to the TAN case in the absence of the DNase I (⁺TAN) (left panel in Fig 5B). Thus, simulations are in good agreement with experimental data [22]. By definition, NETs are associated with neutrophil proteases with the extracellular histone-bound DNA [88]. In the experiments [22], pro-invasive effects of NETs were shown to be associated with protease activities of NET-associated protease, NE. Park et al. [22] found that the NE inhibitor reduced the extension of cancer cell-induced NETs and inhibit TAN’s ability to promote the invasion of 4T1 and BT-549 breast cancer cells. They also found that DNase I treatment can also prevent lung metastasis in mice. However, it is worth observing that this DNase I is not enough to completely inhibit the aggressive migration of tumor cells from the upper chamber to the lower chamber, since they are able to invade in the absence of neutrophils.

In Fig 6, we illustrate the TGF-β-mediated transition between N1 and N2 neutrophils, thus giving rise to two phenotypic states: (a) state I, dominated by non-invasive cancer cells, and (b) state II, dominated by invasive cancer cells. A conceptual schematic of cancer-immune interplay is shown in Fig 6A. Fig 6B–6D shows the scaled populations of invasive cancer cells and, N1 and N2 neutrophils for various levels of TGF-β (0, 0.0001, 0.001, 0.002, 0.01, 5, 10, 20, 50, 100, 200). When the TGF-β level is low, N1 neutrophils dominate the tumor microenvironment (Fig 6C and 6D), leading to non-invasive states of cancer cells (State I, Fig 6B). As the TGF-β level is increased, the N1-dominant system transits to the N2-dominant state (Fig 6C and 6D), resulting in the invasive state of cancer cells (State II, Fig 6B). These results illustrate that the positive feedback loop between N2 neutrophils and invasive phenotype via up-regulation of TGF-β, chemotaxis through CXCL8 and NET activities essentially determines the phenotypic transition between non-invasive and invasive states.

Fig 6

TGF-β-mediated cancer-TAN interplay can induce two types of phenotypic states: invasive and non-invasive types.

(A) Conceptual interaction network for the mathematical model. (B) Population of invasive tumor cells in the lower chamber as a function of TGF-β. (C,D) Populations of N1 (C) and N2 (D) neutrophils as a function of TGF-β.

The N1→N2 transition of TANs was shown to play a critical role in promoting tumor growth, angiogenesis, invasion [3, 24, 25, 89], and ultimately metastasis initiation [90, 91]. Fig 7A shows the spatial profiles of the tumor cells for various N1→N2 transition rates (λ₁₂ = 1.6 × 10⁻⁴ (red solid), 1.6 × 10⁻² (blue dashed), 1.6 × 10⁻¹ (pink with marks)) at the final time (t = 22 h). The N1- and N2-dominant spatial profiles of neutrophils in the lower chamber for the corresponding parameter set are shown in Fig 7B and 7C, respectively. If we increase the rate λ₁₂ (differentiation degree of anti-tumorigenic TANs to tumor-promoting TANs), the N2 population dominates the lower chamber (Fig 7B and 7C) with the higher population ratio N2:N1 of TANs (Fig 7F). Fig 7D shows time courses of NE levels for various values of the differentiation rate (λ₁₂ = 1.6 × 10⁻⁴, 1.6 × 10⁻³, 1.6 × 10⁻², 1.6 × 10⁻¹). The corresponding populations of invasive tumor cells and neutrophils (N1,N2) at final time (t = 22 h) are shown in Fig 7E and 7F, respectively. For a larger λ₁₂, more aggressive, tumor-promoting N2 TANs in the lower chamber can interact with tumor cells in the upper chamber (Fig 7A) by secreting more NE (Fig 7D). This leads to an increased tumor population and enhanced tumor cell invasion (Fig 7E). This increased invasiveness of the tumor cells is the result of the mutual interactions between tumor cells in the upper well and the neutrophils in the lower well. For instance, the NET/NE level increases as λ₁₂ increases (Fig 7D). For a large λ₁₂, the most of N1 TANs are converted into the N2 phenotype (4th column (λ₁₂ = 1.6 × 10⁻¹) in Fig 7F), leading to efficient tumor cell migration (Fig 7E). However, when this transition rate is small (λ₁₂ = 1.6 × 10⁻⁴), the less effective N1 TANs persist in the lower chamber (Fig 7B) with less population of the N2 phenotype (Fig 7C). This results in the slower (or close to zero) production of NE (Fig 7D) by TANs, and lower secretion of both TGF-β and MMP by tumor cells, which in turn reduces invasiveness of tumor cells by more than 48% (Fig 7E).

Fig 7

Effect of the N1→N2 transformation on tumor invasion and N1/N2 dynamics.

(A) Tumor density profiles on Ω = [0, 1] at the final time (t = 22 h) for various λ₁₂’s (λ₁₂ = 1.6 × 10⁻⁴, 1.6 × 10⁻², 1.6 × 10⁻¹). (B,C) Density profiles of the N1 and N2 TANs in the lower chamber ([0, 0.5]) for the corresponding λ₁₂’s in (A). (D) Time courses of NE levels for various values of the differentiation rate (λ₁₂ = 1.6 × 10⁻⁴, 1.6 × 10⁻³, 1.6 × 10⁻², 1.6 × 10⁻¹). (E,F) Scaled population of invasive tumor cells and neutrophils (N1 and N2) at the final time (t = 22 h) for various λ₁₂’s in (D).

Application of the model

Blocking TGF-β and its receptors was shown to inhibit tumor growth [92, 93] and critical cell invasion [34, 94, 95], decrease tumourigenic potential [92, 96], and reduce metastatic incidence [97] through many different pathways [67, 98]. In order to investigate the effect of TGF-β suppression on tumor cell invasion, we introduce a new variable A(x, t) for the TGF-β anti-body and derive the following equations including the modified equation of TGF-β from Eq (11)

Fig 8

The effect of TGF-β blocking (+Ab) and the combined therapy (+Ab+DNase I) on tumor cell invasion.

(A) The (relative) population of migrating tumor cells for various growth rates (r ∈ [1.2, 1.5]) and injection rates (λ_A ∈ [0, 10]) of the TGF-β antibody. When TGF-β activity is inhibited by antibody (r = 1.5), fewer cells (62% reduction) invade the lower chamber. (B) Population of migrating tumor cells when the TGF-β antibody was added in the absence (+Ab-D) and presence (+Ab+D) of DNase I relative to the control (-Ab-D). When proteolytic activity of tumor cells near the membrane is blocked by DNase I, fewer cells (66% reduction) invade the lower chamber in the presence of TGF-β inhibitor.

Next, we tested if MMP inhibition by TIMP can effectively reduce the TAN-mediated invasion through the transfilter. It has been shown that TIMP treatment can significantly reduce (> 50%) the number of migratory breast cancer cells through 8-μm pores in a Boyden invasion chamber assay [99]. In the mathematical model, inhibition of MMP is implemented by injecting TIMP (λ_M > 0) in Eq (15) which abrogates MMP production through a term of degradation of MMPs, $$ -{\mu }_{PM}\frac{P{M}^m}{K_M^m+M^m}$$, on the right hand side of Eq (14). We also tested if a combination therapy (TIMP+TGF-β inhibitor) can further reduce the invasion potential of tumor cells. Fig 9 shows the (relative) MMP levels, ECM levels, and number of migratory tumor cells at final time (t = 22 h) for control (-TIMP-Ab), TIMP alone (+TIMP-Ab), and combination treatment (+TIMP+Ab). One can see that TIMP treatment can inhibit the tumor cell motility by 26% (+TIMP-Ab in Fig 9C) through the significant reduction in MMP activities (82%; Fig 9A) and relatively intact ECM level (Fig 9B). An introduction of the TGF-β antibody to the system in addition to TIMP (+TIMP+Ab) significantly reduces the tumor cell migration through the filter (> 87%). Blocking tumor-TAN interactions by TGF-β inhibitor effectively lowered MMP levels (Fig 9A) and induced the intact ECM levels on the membrane (Fig 9B) [94, 97], contributing to this anti-invasion effect.

Fig 9

The effect of MMP blocking (+TIMP-Ab) and combined therapy (+TIMP+Ab).

(A,B) Levels of MMPs and ECM when MMP activity was blocked by TIMP in the absence (+TIMP-Ab) and presence (+TIMP+Ab) of the TGF-β antibody relative to the control (-TIMP-Ab). (C) Population of invading tumor cells corresponding to control (-TIMP-Ab), TIMP treatment (+TIMP-Ab), and combined therapy (+TIMP+Ab), respectively.

We now investigate the effect of CXCL-8 on tumor growth and invasion in Fig 10. In our model, CXCL-8 knockdown leads to critical reduction of chemotactic movement of the neutrophils, which in turn reduces tumor growth (blue curve; CXCL8 KO, Fig 10A) and invasion (red; CXCL8 KO, Fig 10B) compared to control due to decreased activities of NE (∼40%) and TGF-β (∼50%). These CXCL-KO-mediated reductions in invasive and proliferative potential of tumor cells were consistently observed in experiments. Kumar et al. [100] showed that inhibition of CXCL-8 leads to a drastic decrease (∼14-fold) in proliferation of LS174T cells (shCXCL8 in Fig 10A) and significant abrogation of tumor cell invasion (blue in Fig 10B) relative control.

Fig 10

Inhibition of CXCL8 reduces proliferation, viability and invasion (Experiments vs simulation results).

(A) Time course of cell proliferation shows a drastic decrease in the LS174T cell population with CXCL8 knockdown (shCXCL8; blue) compared to control (LS174T) in in vivo experiments [100]. Model simulation shows a consistent, significant reduction of the tumor cell proliferation in the CXCL knockdown case (CXCL8-KO; blue) compared to control, abrogating tumor growth. (B) CXCL8 knockdown significantly decreases the invasiveness of LS174T cells (blue) [100], as our model simulation illustrates (red).