Modeling and Simulation of EGF-CSF-1 pathway to Investigate Glioma - Macrophage Interaction in Brain Tumors

Keywords

Glioblastoma; Signaling; EGF; CSF-1; Simulation; Paracrine.

Introduction

Glioblastoma multiform (GBM) is one of the most lethal forms of brain cancer in humans. Survival of patients could be extended up to 15 months with chemotherapy, radiation and surgery [1-4] yet, multiple challenge remains for better clinical outcomes [5, 6]. The biggest difficulties for diagnostic and treatment of GBM underlie in its complex macro and microenvironment. The macro environment of brain presents several complexities such as composing of several sensitive cell types to chemotherapeutic reagents, being surrounded by blood-brain barrier that limits delivery of drugs, and the skull that restricts growth of tumors. On the other hand, microenvironment of GBM is also highly complex, dynamic, hierarchical cell society due to presence of diverse cell types with distinct phenotypes and different proliferative potentials. GBM cells co-evolve with stromal and tumor-associated immune cells (microglia/macrophage) and form complex physical and chemical cell-cell communication network. TAMs are abundant and the predominant infiltrating immune cells in malignant GBMs, which are present at World Health Organization (WHO) grade II-IV gliomas [7-10]. During tumor progression, macrophages can stimulate tumor invasion, angiogenesis, and metastasis [11]. Although the role of infiltrated macrophages in tumor progression has been recognized, still the precise nature of the interaction mechanisms between tumor cells and macrophages has not been elucidated [12]. Mills et al. extended an in vivo model for the function of TAMs and suggested two states of TAMs as activated (M1) and alternatively activated (M2) macrophages [13]. M1 and M2 TAMs differ in activating signals, expression of receptors, cytokine production and biological behavior. This suggestion describes that TAMs with M1 polarization are foes and TAMs with M2 polarization are friends for tumors [7-13].

Tumor-derived molecules, such as colony-stimulating factor 1(CSF-1), can polarize glioma-infiltrating macrophages towards M2 polarization and result in production of anti- inflammatory molecules [14-16] and epidermal growth factor (EGF), which acts in return on EGF receptor (EGFR) on the carcinoma cells to promote invasion [17]. EGF-CSF-1 signaling affects the ratio of cell types in aggregates and enables glioma cells to infiltrate into the brain parenchyma [5]. Since GBMs are highly complex with unpredictable patterns, several mathematical models have been used to reveal its complexity and predict its progress [18, 19]. Particularly, compared to other scientific boards, neuro- oncology still requires more effort to propose predictive tools that could accurately simulate the behavior of malignant gliomas [20, 21]. Martirosyan and his co-workers summarized the mathematical models that describe different aspects of GBM growth and evaluastion such as spheroid models, metabolic and vascular models, morphological models, and treatment models [20].

Among them, spheroid models represent a powerful theoretical framework to study initial growth of GBM when proliferation and diffusion of glioblastoma cells are the major players in the tumor initiation. These types of models composed of reactiondiffusion models, simple discrete models and continuum models. Stein et al. used bright field image sequences to estimate number of cells in the tumor spheroids and described a continuum mathematical model to quantitatively interpret the data [22]. After fitting quantitative and experimental data, they observed that glioma cells with EGF receptor show less cell-cell adhesion and invade in a more biased manner and greater rate [22]. Banerjee et al. developed a mathematical model considering the interactive dynamics of glioma cells, macrophages, cytotoxic T-lymphocytes and T11 Target structure (T11TS),

which is a membrane glycoprotein and affects the functional state of immune cells. Their model concluded that T11T structure might be used as a drug target for effective treatment of brain tumors [6]. Aubert and Bandoual proposed a two dimensional model that quantified the strength of cell-cell adhesion using a probability threshold [23]. The agreement of mathematical modeling with experimental results approved that cell-cell adhesion is extremely important for the growth and behavior of glioma cells [23, 24]. Considering biased diffusion in glioblastoma, Fort and Sole’s improved standard reaction-diffusion-advection model pointed that glioma cells move in a bias towards the invasion front instead of moving equally in all directions. It provides a great agreement with experiments [25].

The vascular and metabolic models are related to invasiveness and aggressiveness of the tumor that requires more nutrient supply consequently; these models oftentimes predict the onset of angiogenesis and creation of vasculature [20]. In this concept, some models analysed collective cell migration, tumor cell spatial distribution, morphology and viability using conservation laws [26-28]. Some models were compartmentalized via dividing tumor cell populations into normal, hypoxic, and necrotic cell groups to cover all dynamics of tumor microenvironment [19, 29, 31]. Some models investigated the phenotypic switch that occurs from proliferative state to invasive state in glioma cells as function of hypoxia [32]. The morphological models uses discrete models and reaction- diffusion models to investigate the microscopic and macroscopic morphological changes, glioma growth, invasion based on cell-cell and cell-extracellular matrix adhesion, hypoxia, chemotaxis, homotype attractions, substrate gradients (glucose, oxygen) and other microenvironmental parameters [33-36]. In the treatment models the ultimate goal is providing solutions for better treatment outcome, prolonging and improving patient life. As mentioned, treatment of glioblastoma is not very efficient compared to other cancer types and better strategies are urgently needed. The modeling strategies target better treatment regimen using radiotherapy, chemotherapy, patient MRI data and resection in conjunction with two-, three- and four-dimensional computer modeling systems [37]. Thus, chemotherapy and radiotherapy is able to simulated with different drugs [38] and dose schedules, partial resection of tumor and filling the ablated volume with different chemicals such as cerebrospinal fluid, chemoattractants [39].

Despite efforts to understand the dynamics of glioma cells and macrophage interactions, little data is available to suggest the partnership between glioma cells and M2 type TAMs. However, most of the findings are based on human end-stage tumor samples obtained from surgical secretions. In order to reveal the nature of interaction between glioma cells and macrophages, efforts may focus on studying the nature of interaction between glioma cells and macrophages arises at tumor onset [8]. Based on the need for better understanding of the macrophage-tumor cell interactions in tumor microenvironment, mathematical models, which reveal and simulate the nature of these interactions, are of high challenge and consideration. In this work, we present a computational model for further investigation of macrophage-glioma cell interactions focusing on concentration change of paracrine-acting agents (CSF-1 and EGF) in a defined microenvironment (domain) and on the cellular surfaces.

Model

Recent clinical experiments reported that TAMs facilitate invasiveness of GBM through EGF- CSF-1 paracrine signaling loop [17]. Macrophages secrete EGF and respond to CSF-1; similarly glioma tumor cells express CSF-1 and respond to EGF via chemotaxis, Figure 1a. This cooperation enables glioma cells to coordinate their aggregation and migration via macrophagefacilitated dissemination from primary tumor to surrounding healthy brain tissue [5, 16, 17]. To investigate the paracrine interaction between glioma cells and macrophages through the EGF and CSF-1 signaling [40]; we present a case study, as illustrated in Figure 1b. We propose that there is a macrophage accumulation at the vicinity of glioma cell population in the tumor microenvironment, that could appear where the blood brain barrier or the blood vessel is impaired and immune cells could colonized at the periphery of a solid tumor mass [41-44] in vivo. While Boyden Chamber assays could be considered as examples in vitro [44], where chemotaxis and cell migration studies have been frequently performed between two different cell types through 3-12 μm pore-sized filters [44]. In our model, the vicinity of two cell types in the tumor microenvironment is labeled as domain. The domain defines boundaries (5 μm x 25 μm) to determine the secreted and bounded signaling reagents and their concentration gradient during the different phases of glioma.

---

---

In our modeling, we numerically solve the equations that govern the movement and binding of CSF-1 and EGF. Fick’s second law explains diffusive transport, where D is the diffusion coefficient. C is the concentration of species, Δ is Laplacian and t is time.

Incorporating reaction-diffusion modeling, the nature of interaction in the paracrine-signaling loop is described as follows [45]:

CSF-1R represents the concentration of receptors on macrophage cell surface for binding of CSF-1 secreted by glioma cells and EGFR represents the concentration of receptors located on glioma cell surface for binding of EGF secreted by macrophages. Kon-C and Kon-E incorporate the rate constants for CSF-1R and EGFR binding, respectively. Moreover, CSF-1s and EGFs are surface species and represent the concentration of bound CSF-1 on the macrophage and bound EGF on the glioma cell surfaces. Kdeg-C and Kdeg-E are the rate constants for degradation of CSF-1s and EGFs. DPC and DPE are the degradation of products. The CSF-1d and EGFd are bulk species in the domain and introduced at the rate of ∂C and ∂E at specific locations and have effective diffusion coefficient of DC and DE, respectively. In our calculations, including surface reaction and bulk diffusion expressions (2) and (3) are described with the following equations [45]:

Equations (4) and (6) are surface-reaction expressions and include the concentrations of free species (CSF-1d and EGFd) and should be solved in combination with the mass balance of species in the domain. The coupling between bulk and surface expressions is obtained as boundary condition in the bulk’s mass expressions (equations (4) and (6)), which sets the flux of CSF-1d and EGFd at the active surfaces.

Our mathematical model is based on the following assumptions:

• The physical and chemical properties of domain and surface CSF-1 and EGF. eagents are uniform and continuous.
• The mass balance of the domain has been coupled to the mass balances of CSF-1 and EGF present on the active surfaces.
• The initial condition concentrations of CSF-1 and EGF are zero.

For the domain species, the boundary conditions at active surfaces couple the rate of the reactions at the surfaces with the concentration of free species in the domain:

In order to couple the reaction-diffusion expressions of CSF- 1 and EGF, we have assumed that Kon-E varies based on the concentration of macrophages and glioma cells in the domain and a linear correlation occurs between ∂C, ∂E, Kon-C and Kon-E as follows,

Therefore, the value of Kon-E for each glioma grade could be determined from equation (10). Definitions, default values and their references are provided in Table 1 [46-50].

---

---

The mathematical modeling deals with a diffusion occurring in a 2D domain, which is coupled, to a surface reaction phenomenon occur on a part of the domain’s boundary. The phenomenon in the domain refers to introduction of CSF-1d from the glioma and EGFd from the macrophages as sources and the surface phenomenon describes the binding of species from the domain to the active surfaces that reactions take place. Equations (4) and (6) are modeled using Transport of Diluted Species interface and equations (5 and 7) are described with General Form Boundary PDE interface in COMSOL Multiphysics 5. The equations of two interfaces are coupled considering expressions (8) and (9) as boundary conditions.

The concentration of species in the domain and on the surfaces of the cells has been affected from the interaction of macrophages and glioma cells that are present in the domain. In this study, WHO grade I-IV gliomas are considered based on the concentration of glioma and macrophages. Table 2 describes number of macrophages and glioma cells at each grade of glioma brain cancer [51]. It is important to point that, in our model given the nature of paracrine loop, the number of macrophages decreases with WHO grade I-IV gliomas, while the number of tumor-cells increases. One of the reasons for this assumption in our scenario is the tumor cells are continuously dividing and increasing their number, however once the monocytes differentiated and became active immune cells (macrophages) they cannot divide. Besides, it is assumed that all macrophages are infiltrating into tumor mass from the surrounding tissue, (vasculogenesis and angiogenesis are ignored) when the grade of the tumor increases, the number of macrophages that can diffuse inside the tumor decreases as well as the rate of replenishment of microglia from monocytes reduces. Moreover, microglia stimulates glioblastoma invasion [55] and proliferation; as a consequence glioma cells grow fast. To sum up, in our model the number of macrophages decreases with increased WHO grade of glioma [52-58]. We calculated the value of Kon-E using equation (10) and using the number of macrophages [52] and glioma cells from Table 2 [59].

Due to simplicity and being a general assumption for mathematical models a linear relation between the number of cells and Kon-E value has been chosen in Eq. 10. However, the model can be improved and more realistic studies could be performed considering amplification and overexpression of EGFR receptors, dynamics of receptor changes at the cellular surfaces due to dynamics of cell population such as epithelial to mesenchymal transition, receptor dysregulations, level of autophosphorylation and frequency of mutations [59-63].

---

---

Results and Discussion

In this work, we have developed a continuum mathematical model that simulates the concentrations of EGF and CSF-1 paracrine reagents for the surface-bounded species on the active surfaces of macrophages, CSF-1s (mol/m2) and on the active surfaces of glioma tumor cells, EGFs (mol/m2), the secreted CSF-1d (mol/ m3) and EGFd (mol/m3) as bulk concentrations (free reagents) in the domain. The simulations were performed for 24 hours. Figure 2 shows the schematic view and the change of CSF-1d after 24 hours of interaction between macrophages and glioma tumor cells in the domain for all WHO grades of gliomas (I (a), II (b), III (c) and IV (d)). According to modeling results, all color tables of Figure 2 are assigned to have 1.76×10-26 as minimum (blue) and 3.45×10-16 (red) as maximum reference values for CSF-1d.

---

---

Likewise, Figure 3 shows the concentration gradient of EGFd after 24 hours interaction among macrophages and glioma tumor cells in the domain. The reference minimum value of EGFd is 4.25×10-30 (blue) and the reference maximum value is 2.67×10-17 (red) for the color table. The schematic images represent that the ratio for number of glioma cells to number of macrophages. Grade I glioma is close to 1; this ratio increases from grade I to grade IV, and reaches 260 at grade IV glioma.

One-dimensional plots in Figure 2 and Figure 3 illustrate the change in the concentration of CSF-1d and EGFd in x-direction of the domain for each stage. The x-axis represents the average displacement between the population of glioma cells and macrophages; it starts from the active surface of glioma cells (x = 0 μm) and ends at the active surface of macrophages (x = 5 μm). Consequently, the concentration of CSF-1d decreases through the x-direction of the domain and reaches the approximate value of zero at x = 5 μm.

---

---

When the concentrations of adsorbed species on the cellular surface increase with time, the concentrations of domain species (CSF-1d and EGFd) decrease due to adsorption by active surfaces and coupling of reaction and diffusion phenomena at all glioma grades. For the grade I glioma, the number of macrophages was 1.25 times higher than the number of glioma cells in the domain. At grade II gliomas, macrophages are 0.17 times of glioma cells in the domain and the concentration of CSF-1s and CSF-1d was approximately 5.5 times of the concentration of EGFs and EGFd. For grade III and IV gliomas the ratio of macrophages to glioma cells are 0.01 and 0.003 and the ratio of CSF-1 species to EGF species are approximately 32 and 260, respectively. From grade I to grade IV gliomas, the ratio of macrophages to glioma cells uniformly decreased, but based on the observations from the modeling this approximate uniformity could not be generalized to the ratio of domain species to surface species at each grade of glioma.

Figure 4 shows that for both species of CSF-1, grade I glioma has the minimum concentration and towards high-grade glioma, which means increasing in the ratio of glioma cells to macrophages, CSF-1 bulk and surface species experience higher concentrations. Despite the CSF-1 species, Figure 5 illustrates that EGF species have the minimum concentrations for grade IV gliomas where the ratio of macrophages to glioma cells is the minimum.

---

---

---

Conclusion

There is a growing recognition in the literature regards the role of microglia and macrophages in glioma maintenance and progression [66]; however, we still could not clearly reveal which factors trigger the interaction between glioma cells and macrophages. TAMs infiltrate into the glioma microenvironment, secrete various growth factors, cytokines and chemokine, perform both immune and tumor supportive behaviors. And, in this relation one of the most studied interaction ways for TAMs and glioma cells is CSF-1 and EGF paracrine loop. In glioma microenvironment, CSF-1 is a chemoattractant for TAMs and enrols in their differentiation, migration and survival [64]. In addition, TAMs release EGF, which facilitates glioblastoma proliferation, invasion and recurrence [5]. Although we could not entirely uncover the underlying interaction mechanisms between CSF-1 and EGF, due to urgently needed new targets for glioma treatment, researchers stared to develop TAM-targeted glioma therapies [66-73]. Pyonteck et al., found that CSF-1R inhibition blocks glioma unexpected growth, progression and invasion [64]. Therefore, CSF-1 has become one of the promising targets in the glioblastoma multiform and its reproduction and adsorption expected to be increased from low to high-grade gliomas. Our mathematical modeling approves the high reproduction of CSF-1 at high-grade gliomas as reported by Coniglio and his colleagues [5]. They observed that CSF-1 levels were elevated in higher-grade gliomas and approved that glioblastoma invasion completely depended on CSF-1R signalling [5]. Bender et al., proposed that in a genetic screen for oncogenes driving astrocytomas, CSF-1 was regulated in nearly 70% of spontaneous astrocytomas [64]. Our contribution is presenting a simple mathematical model that accurately determines paracrine loop between macrophages and glioma cells through CSF-1 and EGF interaction for grades I to IV human glioma disease. Our simulation results represent the difference in the surface-bounded (CSF-1s - EGFs) and bulk expressed (CSF-1d - EGFd) paracrine-signaling reagents of CSF-1 and EGF both at spatial and temporal resolutions. Most of the current models focus only on measuring the total concentration of signalling reagents in a constant domain. We improved these models and presented the behaviour of surface-bounded signalling molecules and free signalling reagents of the domain. The results clearly showed the concentration change of signalling molecules in time and their gradient distribution in the domain.

In our mathematical model, increasing the WHO grades of glioblastoma increases the concentration of signalling molecules both in the domain and on the surfaces of the cells. Glioma cells express more CSF-1 in the bulk and macrophages being affected from this high supply, adsorb more CSF-1 on their surfaces until they reach the steady state. Although in this model we used a linear relation between the malignancy grades of gliomas and signalling molecule-binding rates (Kon_E, Kon_C), the model can be improved to be more realistic considering the amplification and overexpression of the signalling receptors, receptor dysregulations, levels of autophosphorylation, and dynamics of cell population (frequency of mutations, receptor changes due to epithelial to mesenchymal transitions, etc.) [61, 62].

Moreover, in our model the increase in the concentration of CSF- 1 and decrease in the concentration of EGF do not obey the uniform change in the number of macrophages and glioma cells. In other words, we aimed to show that not only the increased number of glioma cells due to high-grade glioma contributes to increased CSF-1 concentration but also paracrine interactions [51] between macrophages and glioma cells modulate the CSF- 1 concentration, which is created via coupling the reactiondiffusion expressions of CSF-1 and EGF in our model. Likewise, in our scenario, the decreased number of macrophages compared to increased number of tumor cells with the grade of glioma is different compared to majority of the modeling studies. The reasons behind our assumption were once the monocytes differentiated and become macrophages they cannot proliferate like tumor cells. Besides, when the grade of tumor increases, concentrations of diffused reagents into the tumor mass decrease [74], the number of infiltrated macrophages decreases due to ignored vasculogenesis and angiogenesis in the tumor microenvironment.

Based on our scenario the model accurately generates consistent results with expectations. Our goal was mimicking a macrophage accumulation at the vicinity of glioma cell population in the tumor microenvironment, that could appear where the blood brain barrier or the blood vessel is impaired and immune cells could colonized at the periphery of a solid tumor mass [41-44] in vivo. While Boyden Chamber assays could be considered as examples in vitro [44]. However, this model allows introducing different types of cytokines (more signalling pathways) and cell types (glioma stem cells, astrocytes, and microglia) of glioma microenvironment while providing possibilities of upgrading cellular interaction domain from 2D to 3D microenvironment, where complexity, dynamic cellular distribution of different cell types and heterogeneity of glioma microenvironment will be more realistically mimicked and the obtained results will be extensively contributed to the development of personalized treatment and drug test models for human brain tumors.

Last but not least, in order to develop more predictive and realistic models to investigate tumor-immune cell interactions, we should experimentally generate more quantitative data with more sampling time. Unfortunately, as it has been reported in D. Hambardzumyan, D. H. Gutmann and H. Kettenmann’s review paper, majority of the obtained data in neuro-oncology are not consistent [66]. Therefore, more attention should be given when the models will be tested using the data from the literature.

Acknowledgments

I would like to acknowledge the Faculty of Engineering and Natural Sciences, Sabanci University fund and support.