A model of NSCLC microenvironment predicts optimal receptor targets

Abstract

Background: Tumor microenvironment plays an essential role in the growth of malignancy. Understanding how tumor cells co-evolve with tumor-associated immune cells and stromal cells is important for tumor treatment.

Methods: In this paper, we propose a logistic population dynamics model for quantifying the intercellular signaling network in non-small-cell lung cancer (NSCLC). The model describes the evolutionary dynamics of cells and signaling proteins and was used to predict effective receptor targets through combination strategy analysis. Then, we optimized a multi-target strategy analysis algorithm that was verified by applying it to virtual patients with heterogeneous conditions. Furthermore, to deal with acquired resistance which was commonly observed in patients with NSCLC, we proposed a novel targeting strategy — tracking targeted therapy, to optimize the treatment by improving the therapeutic strategy periodically.

Results: The synergistic effect when inhibiting multiple signaling pathways may help significantly retard carcinogenic processes associated with disease progression, compared with suppression of a single signaling pathway. While traditional treatment (surgery, radiotherapy and chemotherapy) tends to attack tumor cells directly, the multi-target therapy we suggested here is aimed to inhibit the development of tumor by emasculating the relative competitive advantages of tumor cells and promoting that of normal cells.

Conclusion: The combination of traditional and targeted therapy, as an interesting experiment, was significantly more effective in treatment of virtual patients due to a clear complementary relationship between the two therapeutic schemes.

Graphical abstract

Keywords

non-small-cell lung cancer / tumor microenvironment / intercellular signaling network / logistic population dynamics / drug resistance / multi-target therapy

INTRODUCTION

Lung cancer is a significant primary malignant tumor with high incidence rate and fatality rate worldwide. According to the US annual report to the nation, lung cancer leads to more than 1 million deaths per year globally [¹]. Non-small-cell lung cancer (NSCLC), the most universal type of lung carcinoma accounting for 85%, is classified into two main histologic types: adenocarcinoma and squamous carcinoma each accounting for 30% to 50% of NSCLC [²].

Over the past decades, many novel approaches for the treatment of carcinoma have been developed, such as minimally invasive, surgery, targeted therapy, radiation therapy, chemotherapy, and individualized therapy. The therapeutic effects are, however, not satisfactory, owing to the concept of “tumor cell-centered” [³]. In recent years, the rapid growth of tumors has come to be understood that tumors recruit numerous host cells, which are of mesenchymal origin, to constitute the tumor stroma, and these stromal cells greatly accelerate the dramatic increase in the number of malignant cells [⁴]. It also marks a shift in the concept of cancer treatment and a great deal of research has been devoted to developing therapies targeting the tumor microenvironment. Recent evidences suggest that targeted treatment of the tumor stroma contributes to anti-angiogenesis, anti-tumor metastasis and prognosis [⁵]. Therefore, the concept that the malignant microenvironment is a key factor in cancer growth and metastasis has been widely accepted [⁶,⁷].

The normal lung has a wide range of traits, including high degree of hyperoxia, vascularization and homing cell diversity, so it has a unique microenvironment. However, chronic inflammation in the microenvironment alters the normal physiology of the lung. For example, cigarette smoking, the primary cause of lung cancer, may cause long-term inflammatory responses in the pulmonary microenvironment, resulting in the imbalance of cytokine secretion and recruitment of inflammatory cells and further causing malignant transformation of epithelial cells [⁸]. Trials show that an adaptive response of the tumor to the host immune system can create a favorable microenvironment for tumor growth [⁹]. Saffiotti et al. have found in trails that intratracheal administration of silica provokes lung inflammation and causes lung adenocarcinoma in rats [¹⁰]. However, only a small group of animals with pulmonary tumors are created in inhalation assays with cigarette smoke [¹¹].

In recent years, the treatment of advanced NSCLC has long been not limited to traditional chemotherapy and radiotherapy, but also a wide variety of choices in cancer treatment. The emergency of inhibiting key molecular pathways of oncogenesis through mechanism-based targeted therapies has aroused considerable interest. The success of imatinib in the treatment of chronic myelocytic leukemia (CML) showed the effectiveness of targeted therapy [¹²]. Since then, for instance, agents targeting critical parts of cellular signaling pathways involved in tumor growth, angiogenesis, and metastasis have made improvement in the prognosis of advanced NSCLC [¹³]. Besides several targeted drugs for components of the tumor microenvironment have been in clinical trials, more research is devoted to the development of drugs targeting VEGF and its receptors to block angiogenesis in tumor microenvironment. In spite of clinical effectiveness achieved in some tumor types, including renal, colon and liver cancers, results obtained by VEGF-targeted drugs in NSCLC are not as promising [¹⁴]. Therefore, the treatment of NSCLC and other solid tumors requires extensive studies into new targets, such as PDGFR, IL-2, CD11b, CXCR4, c-Kit, GM-CSF and CCL2/CCL12, evaluating their antitumor efficacy [¹⁵].

Despite clinical efficacy from the targeted drugs, a number of studies indicate that solid tumors can acquire resistance and evade the suppression of a single signaling pathway [¹⁶]. For instance, it has been confirmed that previously addicted lung tumor cells are unleashed by an increased TGF-b-dependent IL-6 secretion from their EGFR dependency, which leads to primary and acquired resistance [¹⁷].

In other disease states, there is sufficient evidence to show that combination regimens will be more effective than monotherapy or dual therapy, such as the long-term management of HIV infection [¹⁸]. Compared with monotherapy targeting a single signaling pathway, no doubt combination therapy inhibiting multiple signaling pathways may achieve better outcomes [¹³]. Therefore, that makes it progressively more important to introduce combination therapy to tackle the emergence of resistance in cancer treatment [¹³].

A number of models have been established to study the mechanical interactions of cell-cell and cell-ECM in the solid tumor microenvironment. For instance, many multiphase mixture models have been used to simulate the interactions of a variety of tumor cell species, but these models in silico failed to describe a detailed cellular signaling network in tumor microenvironment [¹⁹–³³]. Furthermore, they are prone to study the impact of specific targets on tumor proliferation, progression and metastasis, and not meant to seek the optimal target through the target strategy analysis.

In this paper, we formulate a model that combines logistic population growth and the Hill functions, and use the model to describe the development of NSCLC microenvironment, as well as explore the curative effect of combination therapy for NSCLC.

RESULTS

Cellular and cytokine dynamics

We studied dynamics of intercellular signaling network of NSCLC microenvironment in silico and observed six cell types that co-evolved in a non-linear manner. In Figure 1, a typical simulation of in situ NSCLC growth is shown. The growth of lung cancer cells generally went through three unique stages, from initial slow growth, via rapid proliferation, and eventually to malignant phase. The starting cell populations are Epith ({\mathrm{2.8}\times \mathrm{10}}^{\mathrm{7}}/\mathrm{m}\mathrm{L}), M1 macrophages (\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}/\mathrm{m}\mathrm{L}) and M2 macrophages (\mathrm{1}\times {\mathrm{10}}^{\mathrm{4}}/\mathrm{m}\mathrm{L}). The lung cancer cells were initiated from neoplastic transformation of normal epithelial cells. During the growth of lung cancer cells within the space that lung epithelial cells occupied, the lung epithelial cells endeavored to maintain their abundance and functions until they were replaced by lung cancer cells in the late stage. For all experiments in this paper, the threshold for the clinical detection of lung cancer cells is assumed to be \mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}/\mathrm{m}\mathrm{L}, which is consistent with the data of clinical studies [³⁴]. The number of lung cancer cells had been steadily increasing from the initial slow growth to the malignant phase, with a small kink appearing at the beginning of rapid tumor expansion. In NSCLC microenvironment, macrophages are mainly activated into M1 and M2 phenotypes. While M1 macrophages are the dominating phenotype observed in the early stages of tumor development, M2 macrophages are increasingly differentiated from monocytes which are induced by tumor cells, and eventually overtake M1 as they accumulate comparable competitive advantages in the rapid proliferation stage. As shown in Figure2, in the tumor microenvironment, M2 macrophages occupy approximately 70% of the entire macrophage population, which is consistent with experimental observations [³⁵].

Cytokine dynamics are also non-linear and multi-stage. Figure3 shows the normalized dynamics curves of 16 cytokines, which can be classified into three types according to their time traces. Due to the non-monotonic change of M1 macrophages, the expression levels of TNF-a and MCSF reach the trough and crest respectively during the rapid proliferation stage; IFN-g experiences a rapid decay in the second stage because of the inhibition of M2 macrophages; All the other cytokines have the same curve characteristics as lung cancer cell dynamics.

Therapy for surgery, radiotherapy and chemotherapy

Current treatments for NSCLC are typically surgery, systemic cytotoxic chemotherapy and radiation therapy [³⁶]. As a result, following the traditional treatment regimen, we simulated treatment for three virtual patients (Supplementary Table S1) with the sequential combination of surgery, radiotherapy (five times a week for 6 weeks) and chemotherapy (once every three weeks for 15 weeks). When tumor cells reached the diagnostic threshold (\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}/\mathrm{m}\mathrm{L}), patients began to receive treatment, and the level of tumor cells would rapidly descend to a very low concentration, entering a steady growing period, the duration of which depends on molecular profiles of individual patients (green dashed line as in subsequent figure). However, among the three typical virtual patients, tumor cells rebounded and relapsed after a delay of at most about two years.

Therapeutic targeting of receptors in NSCLC microenvironment

In spite of long-term development and extensive application, traditional treatment still vary greatly in therapeutic effects among individuals mainly due to inter-patient heterogeneity, which may originate from tumor microenvironment and its intercellular signaling network. A patient specific treatment is desired to achieve more consistent and significant efficacy.

Here we traced to the root of heterogeneity and proposed an individualized therapy that targeting the heterogeneous intercellular signaling network. The key to this therapy is to find the most sensitive targets that can maximize the efficacy. Single target analyses show marked difference across three virtual patients (Figure 4). The optimal targets for virtual patient 1 and 2 are both “Ul” ({{\displaystyle \mathrm{u}}}_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{I}\mathrm{L}-\mathrm{10}}) targeting IL-10 receptors on lung cancer cells; the best target for virtual patient 3 is “u2” ({{\displaystyle \mathrm{u}}}_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}) targeting VEGF receptors on lung cancer cells. Unfortunately, none of the patients get their deterioration effectively controlled with single target therapy, suggesting a requirement of multi target therapy (Figure 5).

We further designed a four-target therapeutic strategy analysis algorithm. The multi-target therapy significantly prolongs the progression-free survival (PFS) to an extend way over the linear superposition of that with single target therapy, which implies that there are synergistic effects between targets (Figure 6). The trend of relative competitive advantage was reversed and the development of tumor was inhibited until the multi-target treatment was discontinued after a half-year treatment program for all three virtual patients (Figure 7). It is worth noting that, for virtual patient 3, a half-year multi-target treatment was even more effective than a sequential combination of surgery, radiotherapy and chemotherapy in terms of efficacy (Figure 7), suggesting that the multi-target therapy can serve as an effective strategy by reconstructing the tumor microenvironment instead of eliminating tumor cell directly.

Combination therapy: a combination of traditional treatment and targeted therapy

Since multi-target therapy can emasculate the competitive advantages of tumor cells by transforming the tumor microenvironment, it is possible to wipe out all the tumor cells at their source, while traditional treatment tends to attack tumor cells directly, which seems difficult to completely root out tumor cells. A combination of traditional treatment and targeted therapy is anticipated to significantly increase the efficacy. The results of simulation show that the extension of PFS of three virtual patients with a combination therapy could be over 50% more than the sum of that with two solo therapies (Figure 7), suggesting that there is a clear complementary relationship between the two therapeutic schemes.

Targeted therapy for dynamic patients

In the course of tumor development, the change of physical conditions and psychological reactions of patients may become more intense on average compared with healthy individuals. To evaluate the therapeutic effect of targeted therapy in response to this reality, we have designed a dynamic virtual patient with time varying biochemical parameters, reflecting the changes in physical conditions. A portion of parameters of ‘k’ (protein secretion rate) and ‘u’ (maximum up-regulation ratio of cellular activity) in the dynamic virtual patient are changed with the variation of lung cancer cell concentration, which is described by the Hill functions (Supplementary Table S3). To improve the efficacy of targeted therapy taking into account the intense fluctuation of biochemical parameters, we proposed a novel targeting strategy—tracking targeted therapy, which aimed to optimize the treatment by upgrading the treatment regimen periodically. The same dynamic virtual patients was treated with tracking and non-tracking targeted therapy, respectively. Both treatments last for one year, and the variations of their target combinations are noted in Supplementary Table S4. While getting well curative effects from both treatments, the dynamic virtual patient benefits more from the tracking targeted therapy (Figure 8).

Drug resistance for multi-target therapy

Regardless of drug resistance, multi-target therapy could be prolonged to one year or even longer so as to bring about more effective results (green dotted line, Figure 9). However, it is common in clinic that acquired resistance happens after about six months of continued targeted therapy, which leads to a huge reduction of therapeutic effects especially in the treatment of lung cancer (red dashed line, Figure 9). From this perspective, a targeted therapy is actually like a shooting sport with running targets instead of stationary targets. To capture the dynamic feature of moving targets, we put forward tracking targeted therapy again to continue eroding the competitive advantages of lung cancer cells by way of responding timely to fast-changing drug resistance. According to real-time changing biochemical parameters, virtual patients took on a strategy shift every two months (Table 1). The PFS can be well extended as long as the treatment was continued (Figure 10), suggesting the tracking targeted therapy proposed here can effectively avoid the failure caused by acquired resistance. Thus, it is quite promising that the PFS could be extended to a very surprising level with a carefully designed tracking multi-target therapy (data not shown).

DISCUSSION

This paper focuses on quantitatively describing the intercellular signaling network in NSCLC microenvironment and exploring key receptor targets for personalized therapy. We integrated a series of different kinds of cells and their intercellular signal transduction pathways into the intercellular signaling network and displayed the dynamics of cells and cytokines by means of dynamics modeling and simulation. All the in silico results show a good qualitative agreement with that from clinical research and animal experiment [³⁷]. To construct the intercellular signaling network of the tumor microenvironment enabled us to discover the mechanism of tumor-stroma competition and cooperation at the systems level and provided novel biological insights into the therapeutic strategies of NSCLC.

With only 24 ODEs of 6 kinds of cells and 18 cytokines, we depicted the malignant transformation, tumor proliferation and angiogenesis in NSCLC microenvironment from perspective of dynamic evolution and system interaction. Despite of significant inter-patient heterogeneity, dynamics of lung cancer cells in all virtual patients exhibited three typical phases: initial slow growth, rapid proliferation and malignant phase. Not only lung cancer cells, but also other tumor-associated cells in the tumor microenvironment underwent these three stages, indicating a co-evolution of lung cancer cells, tumor-associated immune cells and stromal cells. Remarkably, in the second phase, two critical events occurred in NSCLC microenvironment: First, lung cancer cells gained competitive advantages, eroded the space and nutrient resources of epithelial cells and arrived the stage of swift and violent growth; Second, due to paracrine stimulation from lung cancer cells, M2 macrophages seized the selection advantage in monocyte differentiation and dominated in macrophage population which was accompanied by a rapid decrease in immune response. These results did not deliberately cater to any specific clinical or experimental data, but rather reflected the mechanism of competition and cooperation underlying the growth of tumor.

Our model revisits tumor development from the perspective of tumor microenvironment and thus offer new cut-in points to anti-tumor therapy. Since the model quantitatively describes the intercellular signaling network with a variety of cells as nodes and intercellular signaling pathways as links, a targeted therapy that down-regulating/cutting off one or more of the links may change the topology of the entire network, thus interfere the development of lung cancer. We evaluated the efficacy of targeted therapy which can be realized in clinic by regulating specific signaling pathways with small-molecule drugs. Although the most sensitive receptor target in tumor microenvironment was identified, the efficacy of single-target therapy was trivial in all three virtual patients. Interestingly, we found that regulating multiple signaling pathways simultaneously could greatly increase the efficacy, most likely because of the synergistic effect, and thus proposed an algorithm to find the optimal targets combination. The multi-target treatment indeed effectively retarded the development of NSCLC in virtual patients.

Surgery, radiotherapy and chemotherapy have been widely used in the treatment of NSCLC. While most of traditional therapies treat the patients by directly attacking tumor cells, the targeted therapy we proposed here inhibits the growth of tumor through regulation of tumor microenvironment, instead of directly killing tumor cells. Apparently, there are complementary effects between traditional and targeted therapies, a combination of the two is therefore anticipated to achieve better results. Thus we put forward a new therapy combined with surgery, radio-chemotherapy and multi-target therapy, and the expected satisfactory results have been achieved in all virtual patients.

In recent years, multi drug resistance sites are prone to occur around 6 months after the use of small-molecule drugs, suggesting we are dealing with running targets instead of stationary targets. To this end, we proposed an improved version of the multi-target therapy strategy- tracking targeted therapy, and successfully controlled the disease during treatment. Yet nowadays, targeted drug species are less and mostly focus on tumor cells themselves. This study provides a certain theoretical reference for targeting the tumor microenvironment and quantitatively assessing the effectiveness of multi-target therapy. In the future, with the development of targeted drug research and diagnosis techniques, the model can be treated as a powerful clinical tool for the guidance of personalized combination therapy.

METHODS

Constructing an intercellular signaling network of the NSCLC microenvironment

Although the recognitions and biochemical activities of signal molecules in NSCLC microenvironment have been experimentally investigated [⁴,⁷,³⁸–⁴²], it is not clear how these mediators coordinate and play a role in regulating tumor development at the system level. Here we first set up an intercellular signaling network by integrating all known autocrine/paracrine pathways in NSCLC, as shown in Figure 11. Six types of cells – lung cancer cells (LCC), cancer-associated fibroblasts (CAF), tumor-associated macrophages (M1, M2), mesenchymal stem cells (MSC) and epithelial cells (Epith) – and 18 kinds of cytokines/growth factors/chemokines were included in the signaling network.

Well-mixed model

We took 1ml control volume into consideration and derived a quantitative model using logistic population dynamics and the Hill functions. The whole model is composed of 24 ordinary differential equations (ODEs), which corresponds to each type of cell and cytokine, respectively. First, a class of basic population dynamics equation was used to calculate six types of cell growth rates as a function of their recruitment rate, proliferation rate, decay (apoptosis) rate, rate of transformation through direct mutation, rate of transformation through differentiation of their stem/progenitor cells, and rate of dedifferentiation. Second, the growth rate of each cell type is quantitatively regulated by the soluble signal transmission medium present in the tumor microenvironment via the Hill function. So, the cell growth rate can be described as

(1)$$\begin{gathered} \begin{array}{l}{\displaystyle \frac{\stackrel{\begin{array}{l}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\hspace{0.17em}\mathrm{o}\mathrm{f}\hspace{0.17em}\mathrm{c}\mathrm{h}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}\hspace{0.17em}\mathrm{o}\mathrm{f}\hspace{0.17em}\\ \mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\hspace{0.17em}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\end{array}}{\overbrace{{{\displaystyle \mathrm{d}x}}_i}}}{\mathrm{d}t}}=\stackrel{\begin{array}{l}\mathrm{i}\mathrm{m}\mathrm{m}\mathrm{i}\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n},\hspace{0.17em}\mathrm{e}\mathrm{m}\mathrm{i}\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n},\hspace{0.17em}\\ \mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{n}\hspace{0.17em}\mathrm{f}\mathrm{r}\mathrm{o}\mathrm{m}\hspace{0.17em}\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{o}\mathrm{r}\end{array}}{\overbrace{\underset{k=\mathrm{1}}{\overset{{{\displaystyle N}}_i}{{\displaystyle \Sigma }}}{{\displaystyle Y}}_k}}+\stackrel{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}{\overbrace{{{\displaystyle r}}_i{{\displaystyle x}}_i\left(\mathrm{1}-{\displaystyle \frac{\underset{k=\mathrm{1}}{\overset{N^{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}}{{\displaystyle \Sigma }}}{{\displaystyle x}}_k}{{{\displaystyle x}}_{\max }{{\displaystyle \alpha }}_{\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{s}}}}\right)\underset{l=\mathrm{1}}{\overset{N^{\mathrm{c}\mathrm{y}\mathrm{t}\mathrm{o}\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{e}}}{{\displaystyle \Pi }}}\underset{m=\mathrm{1}}{\overset{N^{\mathrm{c}\mathrm{y}\mathrm{t}\mathrm{o}\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{e}}}{{\displaystyle \Pi }}}\left(\mathrm{1}+{\displaystyle \frac{{{\displaystyle u}}_{il}{{\displaystyle y}}_l}{{{\displaystyle s}}_{il}+{{\displaystyle y}}_l}}\right)\left(\mathrm{1}-{{\displaystyle f}}_{im}+{\displaystyle \frac{{{\displaystyle f}}_{im}{{\displaystyle s}}_{im}}{{{\displaystyle s}}_{il}+{{\displaystyle y}}_m}}\right)}}\\ -\stackrel{\mathrm{d}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{y}}{\overbrace{{{\displaystyle d}}_i{{\displaystyle x}}_i\underset{n=\mathrm{1}}{\overset{N^{\mathrm{c}\mathrm{y}\mathrm{t}\mathrm{o}\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{e}}}{{\displaystyle \Pi }}}\left(\mathrm{1}-{{\displaystyle g}}_{in}+{\displaystyle \frac{{{\displaystyle g}}_{in}{{\displaystyle s}}_{in}}{{{\displaystyle s}}_{in}{{\displaystyle +y}}_n}}\right)}}+\stackrel{\mathrm{m}\mathrm{u}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n},\hspace{0.17em}\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n},\hspace{0.17em}\mathrm{d}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\hspace{0.17em}\mathrm{f}\mathrm{r}\mathrm{o}\mathrm{m}\hspace{0.17em}\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}{{\displaystyle x}}_k}{\overbrace{\underset{k=\mathrm{1}}{\overset{N^{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}}{{\displaystyle \Sigma }}}{{\displaystyle m}}_{ik}{{\displaystyle m}}_{ik}\underset{p=\mathrm{1}}{\overset{N^{\mathrm{c}\mathrm{y}\mathrm{t}\mathrm{o}\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{e}}}{{\displaystyle \Pi }}}\underset{q=\mathrm{1}}{\overset{N^{\mathrm{c}\mathrm{y}\mathrm{t}\mathrm{o}\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{e}}}{{\displaystyle \Pi }}}\left(\mathrm{1}+{\displaystyle \frac{{{\displaystyle u}}_{ikp}{{\displaystyle y}}_p}{{{\displaystyle s}}_{ikp}+{{\displaystyle y}}_p}}\right)\left(\mathrm{1}-{{\displaystyle h}}_{ikq}+{\displaystyle \frac{{{\displaystyle h}}_{ikq}{{\displaystyle s}}_{ikq}}{{{\displaystyle s}}_{ikq}+{{\displaystyle y}}_q}}\right)}}\end{array} \\ i=\mathrm{1},\mathrm{...},N^{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}} \end{gathered}$$

where, the first item of the right side of Equation (1) is the representation of discrete events which cover immigration, emigration, differentiation from progenitor. {{\displaystyle Y}}_k refers to the size of the k^{th} cell activity term. The second item denotes the proliferation component of cell {{\displaystyle x}}_i with the basal proliferation rate {{\displaystyle r}}_i, which can be up-regulated by cytokine {{\displaystyle y}}_l and inhibited by cytokine {{\displaystyle y}}_m. Parameter {{\displaystyle a}}_{angiogenesis} is the angiogenesis factor. {{\displaystyle x}}_{\max } is maximum cell population level due to limited nutrition and space. The third item describes the decay of cell {{\displaystyle x}}_i due to natural lifespan, which can be regulated by cytokine {{\displaystyle y}}_n. The last item is the mutation/differentiation/dedifferentiation from cell {{\displaystyle x}}_k that can be up-regulated by cytokine {{\displaystyle y}}_p and inhibited by cytokine {{\displaystyle y}}_q.

The cytokine growth rate can be described as

(2)\begin{array}{l}{\displaystyle \frac{\stackrel{\begin{array}{l}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\hspace{0.17em}\mathrm{o}\mathrm{f}\hspace{0.17em}\mathrm{c}\mathrm{h}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}\hspace{0.17em}\mathrm{o}\mathrm{f}\hspace{0.17em}\\ \mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\hspace{0.17em}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\end{array}}{\overbrace{{{\displaystyle \mathrm{d}x}}_i}}}{\mathrm{d}t}}=\stackrel{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\hspace{0.17em}\mathrm{f}\mathrm{r}\mathrm{o}\mathrm{m}\hspace{0.17em}\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}{\overbrace{\underset{k=\mathrm{1}}{\overset{N^{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}}{{\displaystyle \Sigma }}}{{\displaystyle k}}_{jk}{{\displaystyle x}}_k\underset{l=\mathrm{1}}{\overset{N^{\mathrm{c}\mathrm{y}\mathrm{t}\mathrm{o}\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{e}}}{{\displaystyle \Pi }}}\underset{m=\mathrm{1}}{\overset{N^{\mathrm{c}\mathrm{y}\mathrm{t}\mathrm{o}\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{e}}}{{\displaystyle \Pi }}}\underset{n=\mathrm{1}}{\overset{N^{\mathrm{c}\mathrm{y}\mathrm{t}\mathrm{o}\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{e}}}{{\displaystyle \Pi }}}\left(\mathrm{1}+{\displaystyle \frac{{{\displaystyle u}}_{jkl}{{\displaystyle y}}_l}{{{\displaystyle s}}_{jkl}+{{\displaystyle y}}_l}}\right)\left(\mathrm{1}-{{\displaystyle v}}_{jkm}+{\displaystyle \frac{{{\displaystyle v}}_{jkm}{{\displaystyle y}}_m}{{{\displaystyle s}}_{jkm}+{{\displaystyle y}}_m}}\right)\left(\mathrm{1}-{{\displaystyle w}}_{jkn}+{\displaystyle \frac{{{\displaystyle w}}_{jkn}{{\displaystyle s}}_n}{{{\displaystyle s}}_{jkm}+{{\displaystyle y}}_n}}\right)}}-\stackrel{\mathrm{d}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{y}}{\overbrace{{{\displaystyle d}}_i{{\displaystyle y}}_j}}\\ j=\mathrm{1},\mathrm{...},N^{\mathrm{c}\mathrm{y}\mathrm{t}\mathrm{o}\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{e}}\end{array}

where, the first item on the right side of Equation (2) describes the secretion of cytokine {{\displaystyle y}}_j from cell {{\displaystyle x}}_k with basal secretion rate {{\displaystyle k}}_{jk}, which can be up-regulated by cytokine {{\displaystyle y}}_l, induced by cytokine {{\displaystyle y}}_m and inhibited by cytokine {{\displaystyle y}}_n. The last item is the decay component of cytokine {{\displaystyle y}}_j.

All differential equations and the initial values of all parameters are described in Supplementary Text S1 and Supplementary Table S5, respectively. The population dynamics equations of LCC and VEGF selected as the paradigms of 6 cell types and 18 cytokines are as follows.

LCC:

(3){{\displaystyle \dot{x}}}_{\mathrm{L}\mathrm{C}\mathrm{C}}={{\displaystyle R}}_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}}{{\displaystyle p}}_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}}+{{\displaystyle R}}_{\mathrm{L}\mathrm{C}\mathrm{C}}-{{\displaystyle d}}_{\mathrm{L}\mathrm{C}\mathrm{C}}\left(\mathrm{1}+{\displaystyle \frac{{{\displaystyle u}}_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{M}\mathrm{1}}{{\displaystyle x}}_{\mathrm{M}\mathrm{2}}}{{{\displaystyle s}}_{\mathrm{M}\mathrm{1}}+{{\displaystyle x}}_{\mathrm{M}\mathrm{1}}}}\right)\left({\displaystyle \frac{{{\displaystyle s}}_{\mathrm{M}\mathrm{2}}}{{{\displaystyle s}}_{\mathrm{M}\mathrm{2}}+{{\displaystyle x}}_{\mathrm{M}\mathrm{2}}}}\right){{\displaystyle x}}_{\mathrm{L}\mathrm{C}\mathrm{C}}

where

R_{\mathrm{L}\mathrm{C}\mathrm{C}}=r_{\mathrm{L}\mathrm{C}\mathrm{C}}{x}_{\mathrm{L}\mathrm{C}\mathrm{C}}(\mathrm{1}-{\displaystyle \frac{x_{\mathrm{L}\mathrm{C}\mathrm{C}}+x_{\mathrm{C}\mathrm{A}\mathrm{F}}+x_{\mathrm{M}\mathrm{1}}+x_{\mathrm{M}\mathrm{2}}+x_{\mathrm{M}\mathrm{S}\mathrm{C}}+x_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}}}{x_{\max }{A}_{\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{s}}}})

(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{E}\mathrm{G}\mathrm{F}}{y}_{EGF}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{E}\mathrm{G}\mathrm{F}}}}{s_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{E}\mathrm{G}\mathrm{F}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{E}\mathrm{G}\mathrm{F}}}+y_{\mathrm{E}\mathrm{G}\mathrm{F}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{E}\mathrm{G}\mathrm{F}}}}})(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{F}\mathrm{G}\mathrm{F}}{y}_{\mathrm{F}\mathrm{G}\mathrm{F}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{F}\mathrm{G}\mathrm{F}}}}{s_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{F}\mathrm{G}\mathrm{F}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{F}\mathrm{G}\mathrm{F}}}+y_{\mathrm{F}\mathrm{G}\mathrm{F}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{F}\mathrm{G}\mathrm{F}}}}})(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{H}\mathrm{G}\mathrm{F}}{y}_{\mathrm{H}\mathrm{G}\mathrm{F}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{H}\mathrm{G}\mathrm{F}}}}{s_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{H}\mathrm{G}\mathrm{F}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{H}\mathrm{G}\mathrm{F}}}+y_{\mathrm{H}\mathrm{G}\mathrm{F}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{H}\mathrm{G}\mathrm{F}}}}})

(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{I}\mathrm{G}\mathrm{F}-\mathrm{1}}{y}_{\mathrm{I}\mathrm{G}\mathrm{F}-\mathrm{1}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{I}\mathrm{G}\mathrm{F}-\mathrm{1}}}}{s_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{I}\mathrm{G}\mathrm{F}-\mathrm{1}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{I}\mathrm{G}\mathrm{F}-\mathrm{1}}}+y_{\mathrm{I}\mathrm{G}\mathrm{F}-\mathrm{1}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{I}\mathrm{G}\mathrm{F}-\mathrm{1}}}}})(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{I}\mathrm{L}-\mathrm{10}}{y}_{IL-\mathrm{10}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{I}\mathrm{L}-\mathrm{10}}}}{s_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{I}\mathrm{L}-\mathrm{10}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{I}\mathrm{L}-\mathrm{10}}}+y_{\mathrm{I}\mathrm{L}-\mathrm{10}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{I}\mathrm{L}-\mathrm{10}}}}})(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}{y}_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}}}{s_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}}+y_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}}}})

R_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}}=r_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}}{x}_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}}(\mathrm{1}-{\displaystyle \frac{x_{\mathrm{L}\mathrm{C}\mathrm{C}}+x_{\mathrm{C}\mathrm{A}\mathrm{F}}+x_{\mathrm{M}\mathrm{1}}+x_{\mathrm{M}\mathrm{2}}+x_{\mathrm{M}\mathrm{S}\mathrm{C}}+x_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}}}{x_{\max }{A}_{\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{s}}}})(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}\_\mathrm{I}\mathrm{L}-\mathrm{8}}{y}_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}\_\mathrm{I}\mathrm{L}-\mathrm{8}}^{n_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}\_\mathrm{I}\mathrm{L}-\mathrm{8}}}}{s_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}\_\mathrm{I}\mathrm{L}-\mathrm{8}}^{n_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}\_\mathrm{I}\mathrm{L}-\mathrm{8}}}+y_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}\_\mathrm{I}\mathrm{L}-\mathrm{8}}^{n_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}\_\mathrm{I}\mathrm{L}-\mathrm{8}}}}})(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{S}\mathrm{D}\mathrm{F}-\mathrm{1}}{y}_{\mathrm{S}\mathrm{D}\mathrm{F}-\mathrm{1}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{S}\mathrm{D}\mathrm{F}-\mathrm{1}}}}{s_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{S}\mathrm{D}\mathrm{F}-\mathrm{1}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{S}\mathrm{D}\mathrm{F}-\mathrm{1}}}+y_{\mathrm{S}\mathrm{D}\mathrm{F}-\mathrm{1}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{S}\mathrm{D}\mathrm{F}-\mathrm{1}}}}})(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{P}\mathrm{G}\mathrm{E}\mathrm{2}}{y}_{\mathrm{P}\mathrm{G}\mathrm{E}\mathrm{2}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{P}\mathrm{G}\mathrm{E}\mathrm{2}}}}{s_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{P}\mathrm{G}\mathrm{E}\mathrm{2}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{P}\mathrm{G}\mathrm{E}\mathrm{2}}}+y_{\mathrm{P}\mathrm{G}\mathrm{E}\mathrm{2}}^{n_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{P}\mathrm{G}\mathrm{E}\mathrm{2}}}}})

(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}\_\mathrm{E}\mathrm{G}\mathrm{F}}{y}_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}\_\mathrm{E}\mathrm{G}\mathrm{F}}^{n_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}\_\mathrm{E}\mathrm{G}\mathrm{F}}}}{s_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}\_\mathrm{E}\mathrm{G}\mathrm{F}}^{n_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}\_\mathrm{E}\mathrm{G}\mathrm{F}}}+y_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}\_\mathrm{E}\mathrm{G}\mathrm{F}}^{n_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}\_\mathrm{E}\mathrm{G}\mathrm{F}}}}})

A_{\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{s}}=(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}{y}_{VEGF}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}}}{s_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}}+y_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}}}})(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{P}\mathrm{D}\mathrm{G}\mathrm{F}}{y}_{\mathrm{P}\mathrm{D}\mathrm{G}\mathrm{F}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{P}\mathrm{D}\mathrm{G}\mathrm{F}}}}{s_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{P}\mathrm{D}\mathrm{G}\mathrm{F}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{P}\mathrm{D}\mathrm{G}\mathrm{F}}}+y_{\mathrm{P}\mathrm{D}\mathrm{G}\mathrm{F}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{P}\mathrm{D}\mathrm{G}\mathrm{F}}}}})(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{F}\mathrm{G}\mathrm{F}}{y}_{\mathrm{F}\mathrm{G}\mathrm{F}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{F}\mathrm{G}\mathrm{F}}}}{s_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{F}\mathrm{G}\mathrm{F}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{F}\mathrm{G}\mathrm{F}}}+y_{\mathrm{F}\mathrm{G}\mathrm{F}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{F}\mathrm{G}\mathrm{F}}}}})

(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{I}\mathrm{L}\mathrm{6}}{y}_{\mathrm{I}\mathrm{L}\mathrm{6}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{I}\mathrm{L}\mathrm{6}}}}{s_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{I}\mathrm{L}\mathrm{6}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{I}\mathrm{L}6}}+y_{\mathrm{I}\mathrm{L}\mathrm{6}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{I}\mathrm{L}\mathrm{6}}}}})(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{I}\mathrm{L}\mathrm{8}}{y}_{\mathrm{I}\mathrm{L}\mathrm{8}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{I}\mathrm{L}\mathrm{8}}}}{s_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{I}\mathrm{L}\mathrm{8}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{I}\mathrm{L}8}}+y_{\mathrm{I}\mathrm{L}\mathrm{8}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{I}\mathrm{L}\mathrm{8}}}}})(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{S}\mathrm{D}\mathrm{F}\mathrm{1}}{y}_{\mathrm{S}\mathrm{D}\mathrm{F}\mathrm{1}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{S}\mathrm{D}\mathrm{F}\mathrm{1}}}}{s_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{S}\mathrm{D}\mathrm{F}\mathrm{1}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{S}\mathrm{D}\mathrm{F}\mathrm{1}}}+y_{\mathrm{S}\mathrm{D}\mathrm{F}\mathrm{1}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{S}\mathrm{D}\mathrm{F}\mathrm{1}}}}})

(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{M}\mathrm{C}\mathrm{P}\mathrm{1}}{y}_{\mathrm{M}\mathrm{C}\mathrm{P}\mathrm{1}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{M}\mathrm{C}\mathrm{P}\mathrm{1}}}}{s_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{M}\mathrm{C}\mathrm{P}\mathrm{1}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{M}\mathrm{C}\mathrm{P}\mathrm{1}}}+y_{\mathrm{M}\mathrm{C}\mathrm{P}\mathrm{1}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{M}\mathrm{C}\mathrm{P}\mathrm{1}}}}})(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{T}\mathrm{G}\mathrm{F}\mathrm{b}}{y}_{\mathrm{T}\mathrm{G}\mathrm{F}\mathrm{b}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{T}\mathrm{G}\mathrm{F}\mathrm{b}}}}{s_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{T}\mathrm{G}\mathrm{F}\mathrm{b}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{T}\mathrm{G}\mathrm{F}\mathrm{b}}}+y_{\mathrm{T}\mathrm{G}\mathrm{F}\mathrm{b}}^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{T}\mathrm{G}\mathrm{F}\mathrm{b}}}}})(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{T}\mathrm{N}\mathrm{F}-\alpha }{y}_{\mathrm{T}\mathrm{N}\mathrm{F}-\alpha }^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{T}\mathrm{N}\mathrm{F}-\alpha }}}{s_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{T}\mathrm{N}\mathrm{F}-\alpha }^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{T}\mathrm{N}\mathrm{F}-\alpha }}+y_{\mathrm{T}\mathrm{N}\mathrm{F}-\alpha }^{n_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\_\mathrm{T}\mathrm{N}\mathrm{F}-\alpha }}}})

where x is the concentration of cell, y is the concentration of cytokine, r is the basal proliferation rate, d is the basal decay rate, u is maximum up-regulation ratio of cellular activity, s is the concentration of cytokine producing half of the maximal regulation effect to cellular activity, n is the Hill coefficient, {{\displaystyle p}}_{\mathrm{L}\mathrm{C}\mathrm{C}\_\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}} is the probability of mutation from dividing Epith to LCC during each cell cycle. R is the logistic proliferation term. Parameter {{\displaystyle x}}_{\max } is the saturation concentration of cell due to limited resource, and {{\displaystyle A}}_{\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{s}} is the angiogenesis factor. So the product {{\displaystyle x}}_{\max }{{\displaystyle A}}_{\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{i}\mathrm{o}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{s}} represents the carrying capacity of the microenvironment.

Human lung cancer is derived from normal epithelial cells through a series of multi-step processes, including continuous genetic and epigenetic abnormalities [⁴³]. M1 macrophages engulf tumor cells and recruit other immune cells [⁴⁴], while M2 macrophages weaken anti-tumor immunity by suppressing cytotoxic T cell response [⁴⁵]. VEGF, EGF, FGF and HGF are highly overexpressed in NSCLC and play important roles in promoting lung cancer cell proliferation [⁴⁶–⁴⁹]. Abundant SDF-1, IGF-1 and IL10 in NSCLC microenvironment are reported to stimulate lung cancer cell growth [⁵⁰–⁵²]. In NSCLC cell lines, PGE2 and IL-8 promote cell proliferation [⁵³,⁵⁴].

VEGF:

{\dot{y}}_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}=k_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{L}\mathrm{C}\mathrm{C}}{x}_{\mathrm{L}\mathrm{C}\mathrm{C}}(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{T}\mathrm{G}\mathrm{F}-\beta \_\mathrm{L}\mathrm{C}\mathrm{C}}{y}_{\mathrm{T}\mathrm{G}\mathrm{F}-\beta }^{n_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{T}\mathrm{G}\mathrm{F}-\beta \_\mathrm{L}\mathrm{C}\mathrm{C}}}}{s_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{T}\mathrm{G}\mathrm{F}-\beta \_\mathrm{L}\mathrm{C}\mathrm{C}}^{n_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{T}\mathrm{G}\mathrm{F}-\beta \_\mathrm{L}\mathrm{C}\mathrm{C}}}+y_{\mathrm{T}\mathrm{G}\mathrm{F}-\beta }^{n_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{T}\mathrm{G}\mathrm{F}-\beta \_\mathrm{L}\mathrm{C}\mathrm{C}}}}})(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{L}\mathrm{C}\mathrm{C}}{x}_{\mathrm{L}\mathrm{C}\mathrm{C}}^{n_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{L}\mathrm{C}\mathrm{C}}}}{s_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{L}\mathrm{C}\mathrm{C}}^{n_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{L}\mathrm{C}\mathrm{C}}}+x_{\mathrm{L}\mathrm{C}\mathrm{C}}^{n_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{L}\mathrm{C}\mathrm{C}}}}})

+k_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{C}\mathrm{A}\mathrm{F}}{x}_{\mathrm{C}\mathrm{A}\mathrm{F}}(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{T}\mathrm{G}\mathrm{F}-\beta \_\mathrm{C}\mathrm{A}\mathrm{F}}{y}_{\mathrm{T}\mathrm{G}\mathrm{F}-\beta }^{n_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{T}\mathrm{G}\mathrm{F}-\beta \_\mathrm{C}\mathrm{A}\mathrm{F}}}}{s_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{T}\mathrm{G}\mathrm{F}-\beta \_\mathrm{C}\mathrm{A}\mathrm{F}}^{n_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{T}\mathrm{G}\mathrm{F}-\beta \_\mathrm{C}\mathrm{A}\mathrm{F}}}+y_{\mathrm{T}\mathrm{G}\mathrm{F}-\beta }^{n_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{T}\mathrm{G}\mathrm{F}-\beta \_\mathrm{C}\mathrm{A}\mathrm{F}}}}})+k_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{M}\mathrm{2}}{x}_{\mathrm{M}\mathrm{2}}+k_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{M}\mathrm{S}\mathrm{C}}{x}_{\mathrm{M}\mathrm{S}\mathrm{C}}(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{H}\mathrm{G}\mathrm{F}\_\mathrm{L}\mathrm{C}\mathrm{C}}{y}_{\mathrm{H}\mathrm{G}\mathrm{F}}^{n_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{H}\mathrm{G}\mathrm{F}\_\mathrm{L}\mathrm{C}\mathrm{C}}}}{s_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{H}\mathrm{G}\mathrm{F}\_\mathrm{L}\mathrm{C}\mathrm{C}}^{n_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{H}\mathrm{G}\mathrm{F}\_\mathrm{L}\mathrm{C}\mathrm{C}}}+y_{\mathrm{H}\mathrm{G}\mathrm{F}}^{n_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{H}\mathrm{G}\mathrm{F}\_\mathrm{L}\mathrm{C}\mathrm{C}}}}})+k_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}}{x}_{\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}}(\mathrm{1}+{\displaystyle \frac{u_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{T}\mathrm{G}\mathrm{F}-\beta \_\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}}{y}_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}-\beta }^{n_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{T}\mathrm{G}\mathrm{F}-\beta \_\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}}}}{s_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{T}\mathrm{G}\mathrm{F}-\beta \_\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}}^{n_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{T}\mathrm{G}\mathrm{F}-\beta \_\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}}}+y_{\mathrm{T}\mathrm{G}\mathrm{F}-\beta }^{n_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}\_\mathrm{T}\mathrm{G}\mathrm{F}-\beta \_\mathrm{E}\mathrm{p}\mathrm{i}\mathrm{t}\mathrm{h}}}}})

(4)-d_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}{y}_{\mathrm{V}\mathrm{E}\mathrm{G}\mathrm{F}}

Cytokine VEGF is highly overexpressed in lung cancer cells [⁵⁵], which is induced by TGF-b and HGF [⁵⁶,⁵⁷]. TGF-b induces VEGF expression by CAFs and human lung epithelial cells [⁵⁶,⁵⁸]. MSCs and M2 macrophages express multiple angiogenic factors, such as VEGF [⁷,⁵⁹]. Also, hypoxia increases VEGF expression in NSCLC cell lines [⁵⁸].

Pharmacokinetics for drug treatment

We designed a novel treatment method by targeting cellular components of NSCLC microenvironment, and used pharmacokinetics to analyze the actual efficacy of the target therapy. First, concentration-time profile of a drug which follows one-compartment model can be described by the following equation.

(5)C_n={\displaystyle \frac{k_{a}F{X}_{\mathrm{0}}}{V(k_a-k)}}({\displaystyle \frac{\mathrm{1}-e^{-nk\tau }}{\mathrm{1}-e^{-k\tau }}}{e}^{-kt}-{\displaystyle \frac{\mathrm{1}-e^{-n{k}_a\tau }}{\mathrm{1}-e^{-k_a\tau }}}{e}^{-k_{a}t})

where {{\displaystyle C}}_n is the concentration at a given time t, n is the dosing frequency, \tau is the dosing interval, t is total drug treatment time, V is the volume of distribution of the central compartment, F is the fraction or percentage of the administered dose that is available to reach the general circulation, {{\displaystyle X}}_{\mathrm{0}} is the administered dose of drug, k is the elimination rate constant, and {{\displaystyle k}}_a is the absorption rate constant.

Lynch, TJ found that mutant EGFR strengthened EGF-dependent activation of mutant EGFR and enhanced self-sensitivity to gefitinib [⁶⁰]. According to their experimental results, we can get the relationship between the target inhibition effect and the blood concentration, which can be described by the following equation

(6)T={\displaystyle \frac{uC}{s+C}}\times \mathrm{100}\%

where T is the target inhibition rate, C=p\cdot {{\displaystyle C}}_n, p is the tissue/plasma drug concentration ratio, u is maximum down-regulation ratio of cellular activity, s is the concentration of drug producing half maximal regulation effect to cellular activity. We simulated the therapeutic effect of gefitinib (Figure 12). The parameter values about gefitinib, with references and meanings, are given in Table 2.

Single target and multi-target analysis

In recent years, advances in tumor-targeted therapy have entered a brand-new era with the development of molecular biology and the further acquaintance of the pathogenesis from the cellular and molecular levels.

For single target therapy, we independently carried out simulated treatment analysis for each target to obtain the optimal target for maximizing the following objective function:

(7)S_k=x_{\mathrm{L}\mathrm{C}\mathrm{C}}(u_k,t)-x_{\mathrm{L}\mathrm{C}\mathrm{C}}(u_{k-\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{p}\mathrm{y}},t),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{...},\mathrm{36}

where S is the difference of cancer cell concentrations, k is target parameter number, {{\displaystyle u}}_{k-\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{p}\mathrm{y}}={{\displaystyle u}}_k\times T, and t is the time of treatment. When t=0, the tumor has just been diagnosed.

To demonstrate the procedure of multi-target treatment, without loss of generality, we took 4-targets combination as an illustration. Each combination was composed of 4 out of 36 targets. The simulated treatment was carried out one by one to obtain the optimal combination for maximizing the following objective function:

(8)S_i=x_{\mathrm{L}\mathrm{C}\mathrm{C}}(u_i,t)-x_{\mathrm{L}\mathrm{C}\mathrm{C}}(u_{i-\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{p}\mathrm{y}},t),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{...},C_{\mathrm{36}}^{\mathrm{4}}

where u_i={[u_{\mathrm{1}},u_{\mathrm{2}},u_{\mathrm{3}},u_{\mathrm{4}}]}_i is the array of 4 targets, i is the serial number of target combination, {{\displaystyle u}}_{i-\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{p}\mathrm{y}}={{\displaystyle u}}_i\times T, all the remaining parameters have been previously defined.

Multi-target combination optimization

The efficacy of the multi-target combination therapy is not a linear superposition of each single target, due to the nonlinear correlations governed by the nonlinear coupled ordinary differential equations. A combination of several suboptimal targets may outplay a group of top-notch players due to synergistic effects (data not shown). Thus, the key rule to find the candidates of the optimal 4-target combination may not be looking for the top players who win the single target contest, but the ones that closely collaborate to demonstrate best efficacy in multi-target analysis. According to this rule, a straightforward procedure will be running through the targets in all possible combinations, however, it will cause large amount of calculation as the model covers 36 receptor targets. To increase the efficiency of the optimization algorithm, we proposed a grouping strategy to reduce the size of effective sampling library. The strategy is based on the assumption that, it is possible to have strong cross correlation only among the target parameters in the same ordinary differential equation (intra-equation) of the model, while the mutual influence of target parameters existing in different equations (inter-equation) are comparably weaker and the synergy can be ignored.

In accordance with the assumption, we organized an elimination tournament to determine the optimal 4-target combination (Figure 13). The 36 competitions were grouped according to intra- and inter-equation correlations (Table 3), and the winners in the knockout phase went on to play in the next round. The optimal 4-target combination was selected in the finals.

References

[1] Jemal, A., Ward, E. M., Johnson, C. J., Cronin, K. A., Ma, J., Ryerson, B., Mariotto, A., Lake, A. J., Wilson, R., Sherman, R. L., (2017) Annual Report to the Nation on the Status of Cancer, 1975-2014, Featuring Survival. J. Natl. Cancer Inst., 109, 22

[2] Langer, C. J., Besse, B., Gualberto, A., Brambilla, E. and Soria, J.-C. (2010) The evolving role of histology in the management of advanced non-small-cell lung cancer. J. Clin. Oncol., 28, 5311–5320

[3] Leong, D., Rai, R., Nguyen, B., Lee, A. and Yip, D. (2014) Advances in adjuvant systemic therapy for non-small-cell lung cancer. World J. Clin. Oncol., 5, 633–645

[4] Weinberg, R. A. (2008) Coevolution in the tumor microenvironment. Nat. Genet., 40, 494–495

[5] Villanueva, M. T. (2014) Microenvironment: the new midfielders in the tumour microenvironment. Nat. Rev. Cancer, 14, 765

[6] Pietras, K. and Ostman, A. (2010) Hallmarks of cancer: interactions with the tumor stroma. Exp. Cell Res., 316, 1324–1331

[7] El-Nikhely, N., Larzabal, L., Seeger, W., Calvo, A. and Savai, R. (2012) Tumor-stromal interactions in lung cancer: novel candidate targets for therapeutic intervention. Expert Opin. Investig. Drugs, 21, 1107–1122

[8] Lee, G., Walser, T. C. and Dubinett, S. M. (2009) Chronic inflammation, chronic obstructive pulmonary disease, and lung cancer. Curr. Opin. Pulm. Med., 15, 303–307

[9] Sautès-Fridman, C., Cherfils-Vicini, J., Damotte, D., Fisson, S., Fridman, W. H., Cremer, I. and Dieu-Nosjean, M.-C. (2011) Tumor microenvironment is multifaceted. Cancer Metastasis Rev., 30, 13–25

[10] Saffiotti, U., Daniel, L. N., Mao, Y., Shi, X., Williams, A. O. and Kaighn, M. E. (1994) Mechanisms of carcinogenesis by crystalline silica in relation to oxygen radicals. Environ. Health Perspect., 102, 159–163

[11] Coggins, C. R. E. (2010) A further review of inhalation studies with cigarette smoke and lung cancer in experimental animals, including transgenic mice. Inhal. Toxicol., 22, 974–983

[12] Druker, B. J., Talpaz, M., Resta, D. J., Peng, B., Buchdunger, E., Ford, J. M., Lydon, N. B., Kantarjian, H., Capdeville, R., Ohno-Jones, S., (2001) Efficacy and safety of a specific inhibitor of the BCR-ABL tyrosine kinase in chronic myeloid leukemia. N. Engl. J. Med., 344, 1031–1037

[13] Pennell, N. A. and Lynch, T. J. Jr. (2009) Combined inhibition of the VEGFR and EGFR signaling pathways in the treatment of NSCLC. Oncologist, 14, 399–411

[14] Sandler, A., Gray, R., Perry, M. C., Brahmer, J., Schiller, J. H., Dowlati, A., Lilenbaum, R. and Johnson, D. H. (2006) Paclitaxel-carboplatin alone or with bevacizumab for non-small-cell lung cancer. N. Engl. J. Med., 355, 2542–2550

[15] Zou, W. (2005) Immunosuppressive networks in the tumour environment and their therapeutic relevance. Nat. Rev. Cancer, 5, 263–274

[16] Engelman, J. A., Zejnullahu, K., Mitsudomi, T., Song, Y., Hyland, C., Park, J. O., Lindeman, N., Gale, C.-M., Zhao, X., Christensen, J., (2007) MET amplification leads to gefitinib resistance in lung cancer by activating ERBB3 signaling. Science, 316, 1039–1043

[17] Yao, Z., Fenoglio, S., Gao, D. C., Camiolo, M., Stiles, B., Lindsted, T., Schlederer, M., Johns, C., Altorki, N., Mittal, V., (2010) TGF-beta IL-6 axis mediates selective and adaptive mechanisms of resistance to molecular targeted therapy in lung cancer. Proc. Natl. Acad. Sci. USA, 107, 15535–15540

[18] Lange, J. M. A. (1995) Triple combinations: present and future. J. Acquir. Immune Defic. Syndr. Hum. Retrovirol., 10, S77–S82

[19] Ward, J. P. and King, J. R. (1997) Mathematical modelling of avascular-tumour growth. IMA J. Math. Appl. Med. Biol., 14, 39–69

[20] Ward, J. P. and King, J. R. (1999) Mathematical modelling of avascular-tumour growth. II: Modelling growth saturation. IMA J. Math. Appl. Med. Biol., 16, 171–211

[21] Breward, C. J. W., Byrne, H. M. and Lewis, C. E. (2002) The role of cell-cell interactions in a two-phase model for avascular tumour growth. J. Math. Biol., 45, 125–152

[22] Breward, C. J. W., Byrne, H. M. and Lewis, C. E. (2003) A multiphase model describing vascular tumour growth. Bull. Math. Biol., 65, 609–640

[23] Araujo, R. P. and McElwain, D. L. S. (2005) A mixture theory for the genesis of residual stresses in growing tissues I: A general formulation. SIAM J. Appl. Math., 65, 1261–1284

[24] Araujo, R. P. and McElwain, D. L. S. (2005) A mixture theory for the genesis of residual stresses in growing tissues II: solutions to the biphasic equations for a multicell spheroid. SIAM J. Appl. Math., 66, 447–467

[25] Macklin, P. and Lowengrub, J. (2007) Nonlinear simulation of the effect of microenvironment on tumor growth. J. Theor. Biol., 245, 677–704

[26] Wise, S. M., Lowengrub, J. S., Frieboes, H. B. and Cristini, V. (2008) Three-dimensional multispecies nonlinear tumor growth‒I Model and numerical method. J. Theor. Biol., 253, 524–543

[27] Macklin, P., McDougall, S., Anderson, A. R. A., Chaplain, M. A. J., Cristini, V. and Lowengrub, J. (2009) Multiscale modelling and nonlinear simulation of vascular tumour growth. J. Math. Biol., 58, 765–798

[28] Frieboes, H. B., Jin, F., Chuang, Y.-L., Wise, S. M., Lowengrub, J. S. and Cristini, V. (2010) Three-dimensional multispecies nonlinear tumor growth-II: Tumor invasion and angiogenesis. J. Theor. Biol., 264, 1254–1278

[29] Hawkins-Daarud, A., van der Zee, K. G. and Oden, J. T. (2012) Numerical simulation of a thermodynamically consistent four-species tumor growth model. Int. J. Numer. Methods Biomed. Eng., 28, 3–24

[30] Chen, Y. and Lowengrub, J. S. (2014) Tumor growth in complex, evolving microenvironmental geometries: a diffuse domain approach. J. Theor. Biol., 361, 14–30

[31] Albano, G., Giorno, V., Román-Román, P., Román-Román, S. and Torres-Ruiz, F. (2015) Estimating and determining the effect of a therapy on tumor dynamics by means of a modified Gompertz diffusion process. J. Theor. Biol., 364, 206–219

[32] Al-Mamun, M. A., Farid, D. M., Ravenhil, L., Hossain, M. A., Fall, C. and Bass, R. (2016) An in silico model to demonstrate the effects of Maspin on cancer cell dynamics. J. Theor. Biol., 388, 37–49

[33] Ng, C. F. and Frieboes, H. B. (2017) Model of vascular desmoplastic multispecies tumor growth. J. Theor. Biol., 430, 245–282

[34] Schreiber, H. and Rowley, D. A. (2010) Cancer. Awakening immunity. Science, 330, 761–762

[35] Ma, J., Liu, L., Che, G., Yu, N., Dai, F. and You, Z. (2010) The M1 form of tumor-associated macrophages in non-small cell lung cancer is positively associated with survival time. BMC Cancer, 10, 112

[36] Spira, A. and Ettinger, D. S. (2004) Multidisciplinary management of lung cancer. N. Engl. J. Med., 350, 379–392

[37] Abulaiti, A., Shintani, Y., Funaki, S., Nakagiri, T., Inoue, M., Sawabata, N., Minami, M. and Okumura, M. (2013) Interaction between non-small-cell lung cancer cells and fibroblasts via enhancement of TGF-b signaling by IL-6. Lung Cancer, 82, 204–213

[38] Anderson, A. R. A., Weaver, A. M., Cummings, P. T. and Quaranta, V. (2006) Tumor morphology and phenotypic evolution driven by selective pressure from the microenvironment. Cell, 127, 905–915

[39] Bremnes, R. M., Camps, C. and Sirera, R. (2006) Angiogenesis in non-small cell lung cancer: the prognostic impact of neoangiogenesis and the cytokines VEGF and bFGF in tumours and blood. Lung Cancer, 51, 143–158

[40] Tlsty, T. D. and Coussens, L. M. (2006) Tumor stroma and regulation of cancer development. Annu. Rev. Pathol., 1, 119–150

[41] Polyak, K., Haviv, I. and Campbell, I. G. (2009) Co-evolution of tumor cells and their microenvironment. Trends Genet., 25, 30–38

[42] Hanahan, D. and Weinberg, R. A. (2011) Hallmarks of cancer: the next generation. Cell, 144, 646–674

[43] Wistuba, I. I., Behrens, C., Virmani, A. K., Mele, G., Milchgrub, S., Girard, L., Fondon, J. W. 3rd, Garner, H. R., McKay, B., Latif, F., (2000) High resolution chromosome 3p allelotyping of human lung cancer and preneoplastic/preinvasive bronchial epithelium reveals multiple, discontinuous sites of 3p allele loss and three regions of frequent breakpoints. Cancer Res., 60, 1949–1960

[44] Sánchez-Martín, L., Estecha, A., Samaniego, R., Sánchez-Ramón, S., Vega, M. Á. and Sánchez-Mateos, P. (2011) The chemokine CXCL12 regulates monocyte-macrophage differentiation and RUNX3 expression. Blood, 117, 88–97

[45] Noy, R. and Pollard, J. W. (2014) Tumor-associated macrophages: from mechanisms to therapy. Immunity, 41, 49–61

[46] Tsakiridis, T., Cutz, J.-C., Singh, G., Hirte, H., Okawara, G., Corbett, T., Sur, R., Cai, W., Whelan, T. and Wright, J. R. (2008) Association of phosphorylated epidermal growth factor receptor with survival in patients with locally advanced non-small cell lung cancer treated with radiotherapy. J. Thorac. Oncol., 3, 716–722

[47] Puri, N. and Salgia, R. (2008) Synergism of EGFR and c-Met pathways, cross-talk and inhibition, in non-small cell lung cancer. J. Carcinog., 7, 9

[48] Marek, L., Ware, K. E., Fritzsche, A., Hercule, P., Helton, W. R., Smith, J. E., McDermott, L. A., Coldren, C. D., Nemenoff, R. A., Merrick, D. T., (2009) Fibroblast growth factor (FGF) and FGF receptor-mediated autocrine signaling in non-small-cell lung cancer cells. Mol. Pharmacol., 75, 196–207

[49] Huang, H. P., Feng, H., Qiao, H. B., Ren, Z. X. and Zhu, G. D. (2015) The prognostic significance of fibroblast growth factor receptor 4 in non-small-cell lung cancer. OncoTargets Ther., 8, 1157–1164

[50] Li, J., Guan, J., Long, X., Wang, Y. and Xiang, X. (2016) mir-1-mediated paracrine effect of cancer-associated fibroblasts on lung cancer cell proliferation and chemoresistance. Oncol. Rep., 35, 3523–3531

[51] Fritz, J. M., Dwyer-Nield, L. D. and Malkinson, A. M. (2011) Stimulation of neoplastic mouse lung cell proliferation by alveolar macrophage-derived, insulin-like growth factor-1 can be blocked by inhibiting MEK and PI3K activation. Mol. Cancer, 10, 76

[52] Hsu, T. I., Wang, Y. C., Hung, C. Y., Yu, C. H., Su, W. C., Chang, W. C. and Hung, J. J. (2016) Positive feedback regulation between IL10 and EGFR promotes lung cancer formation. Oncotarget, 7, 20840–20854

[53] Zhong, X., Fan, Y., Ritzenthaler, J. D., Zhang, W., Wang, K., Zhou, Q. and Roman, J. (2015) Novel link between prostaglandin E2 (PGE2) and cholinergic signaling in lung cancer: The role of c-Jun in PGE2-induced a7 nicotinic acetylcholine receptor expression and tumor cell proliferation. Thorac. Cancer, 6, 488–500

[54] Luppi, F., Longo, A. M., de Boer, W. I., Rabe, K. F. and Hiemstra, P. S. (2007) Interleukin-8 stimulates cell proliferation in non-small cell lung cancer through epidermal growth factor receptor transactivation. Lung Cancer, 56, 25–33

[55] O’Byrne, K. J., Koukourakis, M. I., Giatromanolaki, A., Cox, G., Turley, H., Steward, W. P., Gatter, K. and Harris, A. L. (2000) Vascular endothelial growth factor, platelet-derived endothelial cell growth factor and angiogenesis in non-small-cell lung cancer. Br. J. Cancer, 82, 1427–1432

[56] Pertovaara, L., Kaipainen, A., Mustonen, T., Orpana, A., Ferrara, N., Saksela, O. and Alitalo, K. (1994) Vascular endothelial growth factor is induced in response to transforming growth factor-beta in fibroblastic and epithelial cells. J. Biol. Chem., 269, 6271–6274

[57] Saijo, Y., Tanaka, M., Miki, M., Usui, K., Suzuki, T., Maemondo, M., Hong, X., Tazawa, R., Kikuchi, T., Matsushima, K., (2002) Proinflammatory cytokine IL-1 beta promotes tumor growth of Lewis lung carcinoma by induction of angiogenic factors: in vivo analysis of tumor-stromal interaction. J. Immunol., 169, 469–475

[58] Boussat, S., Eddahibi, S., Coste, A., Fataccioli, V., Gouge, M., Housset, B., Adnot, S. and Maitre, B. (2000) Expression and regulation of vascular endothelial growth factor in human pulmonary epithelial cells. Am. J. Physiol. Lung Cell. Mol. Physiol., 279, L371–L378

[59] Schmidt, T., Ben-Batalla, I., Schultze, A. and Loges, S. (2012) Macrophage-tumor crosstalk: role of TAMR tyrosine kinase receptors and of their ligands. Cell. Mol. Life Sci., 69, 1391–1414

[60] Lynch, T. J., Bell, D. W., Sordella, R., Gurubhagavatula, S., Okimoto, R. A., Brannigan, B. W., Harris, P. L., Haserlat, S. M., Supko, J. G., Haluska, F. G., (2004) Activating mutations in the epidermal growth factor receptor underlying responsiveness of non-small-cell lung cancer to gefitinib. N. Engl. J. Med., 350, 2129–2139

[61] Li, J., Brahmer, J., Messersmith, W., Hidalgo, M. and Baker, S. D. (2006) Binding of gefitinib, an inhibitor of epidermal growth factor receptor-tyrosine kinase, to plasma proteins and blood cells: in vitro and in cancer patients. Invest. New Drugs, 24, 291–297

[62] Peters, S., Zimmermann, S. and Adjei, A. A. (2014) Oral epidermal growth factor receptor tyrosine kinase inhibitors for the treatment of non-small cell lung cancer: comparative pharmacokinetics and drug-drug interactions. Cancer Treat. Rev., 40, 917–926

[63] Swaisland, H. C., Smith, R. P., Laight, A., Kerr, D. J., Ranson, M., Wilder-Smith, C. H. and Duvauchelle, T. (2005) Single-dose clinical pharmacokinetic studies of gefitinib. Clin. Pharmacokinet., 44, 1165–1177

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature