Modeling invasion patterns in the glioblastoma battlefield

Glioblastoma is the most aggressive tumor of the central nervous system, due to its great infiltration capacity. Understanding the mechanisms that regulate the Glioblastoma invasion front is a major challenge with preeminent potential clinical relevances. In the infiltration front, the key features of tumor dynamics relate to biochemical and biomechanical aspects, which result in the extension of cellular protrusions known as tumor microtubes. The coordination of metalloproteases expression, extracellular matrix degradation, and integrin activity emerges as a leading mechanism that facilitates Glioblastoma expansion and infiltration in uncontaminated brain regions. We propose a novel multidisciplinary approach, based on in vivo experiments in Drosophila and mathematical models, that describes the dynamics of active and inactive integrins in relation to matrix metalloprotease concentration and tumor density at the Glioblastoma invasion front. The mathematical model is based on a non-linear system of evolution equations in which the mechanisms leading chemotaxis, haptotaxis, and front dynamics compete with the movement induced by the saturated flux in porous media. This approach is able to capture the relative influences of the involved agents and reproduce the formation of patterns, which drive tumor front evolution. These patterns have the value of providing biomarker information that is related to the direction of the dynamical evolution of the front and based on static measures of proteins in several tumor samples. Furthermore, we consider in our model biomechanical elements, like the tissue porosity, as indicators of the healthy tissue resistance to tumor progression.

Introduction Glioblastoma is the most common, aggressive, and lethal tumor of the central nervous system. It has a glial origin and it is characterized by rapid cell proliferation, great infiltration capacity, and neurological impairment [1]. GB is composed of a heterogeneous genetic landscape of tumor cells [2] that reduces the efficiency of clinical treatments. Current treatments for GB include surgical resection of the solid tumor, radiation therapy, and chemotherapy with Temozolomide. However, GB is resistant to treatment and recurrence is the main cause of mortality [3,4]; the median survival is below 16 months and the incidence is 3/100.000 per year.
GB infiltration in the human brain is a complex phenomenon, influenced by different tumor properties and the tumor microenvironment, including brain-resident cells, bloodbrain barrier, and the immune system [5]. GB develops fronts of invasion towards uncontaminated areas of the brain, and GB cells modify their cytoskeleton components [6] to extend protrusions, known as Tumor Microtubes (TMs) [7]. TMs mediate GB progression, and they interact with the synapses of the neighboring neurons [8,9]. But, how do these fronts occur? Which are the mechanical and chemical processes that characterize these front areas of GB progression? Which are the main agents involved? What is the role of TMs in the tumor front development? Is there any heterogeneity in the spatial distribution of the tumor activity across the tumor domain?
These are some of the questions that we address throughout the paper, providing insights into the role of integrins, metalloproteases, and TMs in GB growth and spread. Our results are supported by experimental measurements in a Drosophila model of human GB in continuous feedback with the evolutionary mathematical model, which predicts behaviors and interactions of the biological agents.

GB in Drosophila
The most frequent genetic lesions in GB include the constitutive activation of phosphatidylinositol 3-kinase (PI3K) and Epidermal Growth Factor Receptor (EGFR) pathways, which drive cellular proliferation and tumor malignancy [10]. Drosophila melanogaster has emerged as one of the most reliable animal models for GB. It is based on genetic mutations in EGFR and PI3K pathways equivalent to the ones found in patients [11]. Glial cells respond to these oncogenic transformations and reproduce all main features of the disease, including glial expansion and invasion, but in a shorter time. This Drosophila model has been used for drug and genetic screenings and the results have been validated in human GB cells [12][13][14].

Agents involved in the battlefield
GB growth and migration are driven by specific signaling pathways as well as interactions between the tumor and its extracellular microenvironment. In our study, we include cell membrane protrusions as driving factors of tumor progression, as well as the cell response to signaling gradients and the interplay with the Extra Cellular Matrix (ECM). A schematic diagram of the described dynamics is given in Fig 1. The dynamics of the tumor cell membrane, including cell protrusions, are fundamental in several processes, among others cell movements, cell transport, and cell exposure to molecular interactions with substrates during GB progression. Generally, cell protrusions are highly dynamic extensions of the plasma membrane involved in cell migration and invasion through the ECM. Different types of protrusions have been identified to contribute to cell spreading, depending on specific contexts, cell types, and microenvironment. Precisely, filopodia are thin, finger-like, and highly dynamic membrane protrusions that have a significant role in mediating intercellular communication and in modulating cell adhesion. These protrusions appear to be required for haptotaxis and chemotaxis [15], i.e., for the cell response to the gradient of insoluble (haptotaxis) and soluble (chemotaxis) components of the tumor microenvironment. Cytonemes are one type of filopodia, first observed in the Drosophila wing imaginal disc [16] (see also [17][18][19][20][21][22][23]). In addition to the referred intensive studies in Drosophila, the importance of cell communication and motility mediated by cytonemes is known in several living systems and it has been largely studied in relation to various signaling pathways in vertebrates (see the seminal paper by Sanders et al [21], the recent results [24][25][26][27] and the references therein).
Also known as TMs in the context of GB, cytonemes have an important role in tumor development [12]. Specifically, we focus on TM involvement in the GB front progression and on the relationship of TMs with integrins and proteases.
Tumor propagation is characterized by a sharp invasion front located in the front area of GB progression, where TMs are mostly concentrated. The tumor front is a part of the tumor region in contact with the healthy tissue. It is characterized by the presence of TMs/cytonemes and collects a wide activity related to cellular communication signals and cell-ECM interactions. The agents involved in GB invasion, such as integrins, proteases, or the tumor cells themselves, are neither scattered nor randomly moving, but rather there is self-organization that determines invasion patterns around the tumor front. This crucial aspect of the GB progression can be observed in Fig 2. In this figure, we notice how proteases, identified by the anti-MMP1 (Matrix Metalloprotease 1) antibody, and integrins, identified by anti-FAK (Focal Adhesion Kinase) and anti-Talin antibodies, co-localize with the TM region at the GB front. The prediction of these expansion patterns is a need to envisage possible tumor outcomes in patients and for the development of potential therapies. Therefore, any mathematical model that attempts to predict GB dynamics and reproduce the formation of these evolutionary patterns must face these challenges.
Integrins are transmembrane glycoproteins known to mediate dynamic interactions between the ECM and the actin cytoskeleton involved in cell motility. The binding specificity is determined by the integrin extracellular domain, which allows the recognition of different matrix ligands (fibronectin, collagen, laminin) and other cell surface receptors [28]. In Drosophila, the gene rhea encodes Talin, a large adaptor protein required for integrin function. Talin links clusters of integrins, bound to the ECM, and the cytoskeleton ant. Therefore, it is essential for the formation of focal adhesion-like clusters of integrins [29]. Integrin activity and cytoskeleton dynamics are mediated by focal adhesion kinase, a cytoplasmic tyrosine kinase driving the cellular adhesion and spreading processes. Integrin interactions with the ECM components support cell migration [30], but affect also cell growth and division through additional interactions with some extracellular proteins and enzymes that control the cell Glial membrane is shown in red (A 1 ), and MMP1 protein in green (A 2 ). MMP1 accumulation in the GB front is observed in the merged images (A 3 ). B) Glial membrane (B 1 , red) and FAK distribution (B 2 , green) signals are shown through the brain. FAK accumulation in the GB front is detailed in the merged images (B 3 ). C) Co-staining of FAK (C 1 , green) and MMP1 (C 2 , magenta), and merged images (C 3 ). FAK and MMP1 signals accumulate in the same region of GB front (inset in C 1 , red). White arrows in the glial membrane images indicate the tumor front, where protein accumulation is observed. https://doi.org/10.1371/journal.pcbi.1008632.g002

PLOS COMPUTATIONAL BIOLOGY
cycle [31,32]. Integrins are generally present in the cellular membrane, but they are mostly concentrated at the tumor front, where these receptors are directly involved in the signaling leading the migration process.
Proteases are enzymes that catalyze the breakdown of proteins into smaller polypeptides or single amino acids. Some types of proteases are localized around the cell membrane and play a key role in promoting tumor invasion and tissue remodeling. They induce proteolysis of the ECM components [33] and maintain a microenvironment that facilitates tumor cell survival. Protease profiling studies have indicated that the expression of serine proteases, cysteine proteases, and matrix metalloproteases, the most prominent family of proteinases associated with tumorigenesis, increases in high-grade astrocytoma compared with the normal brain. In particular, increased MMP levels and tumor invasiveness in human GBs show a strong correlation [12,34,35]. Several studies reported the protease localization in TM membrane [36,37], especially for MMP-2 and MT1-MMP.
Our interest arises from the analysis of the role of TMs in tumor progression and the effect of different drugs on TM structures and dynamics [38,39], although we are not modeling drug effects in the present work. The inhibition of TMs leads to a decrease in cell migration and proliferation, with obvious consequences on GB treatment [38][39][40]. The process of growth and retraction of TMs plays a fundamental role in cell-to-cell communication routing [17,18,[20][21][22][23], as well as it constitutes a transient binding platform for essential proteins that regulate tumor dynamics [41]. Among these proteins, the presence of the microtubule plus end-binding protein EB1 correlates with GB progression and poor survival [40,42]. This discover has led to a great interest in the development of drugs aimed at affecting end-binding protein expression and, consequently, TM dynamics [38]. Moreover, several experiments have reported the influence of proteases on GB development and the role of several integrin receptors in tumor progression (see [43,44], S1 and S3 Text and S1 Fig for further details).

Integration of models and data
The ability of mathematical models to predict and guide biological experiments has gathered a great attention. The number of possible factors involved in GB progression is large and the relevance of each of them is still unknown. Thus, exploiting the advances in microscopy and protein signaling, there is a need to integrate the large amount of provided experimental data into mathematical models.
Different mathematical models describe the dynamics of tumors by linear (Fick's law) or nonlinear (Darcy's law) diffusion. Some approaches couple the dynamics for tumor cells and environmental agents involved in its evolution, some others consider them independently (e.g. see [45,46] and references therein). However, the arising of sharp profiles is not compatible with a movement induced by linear diffusive dynamics, although some mathematical models explored this possibility as a first approach.
We focus on the dynamics of the tumor propagation front and the emergence of the coordination between self-organized collective processes. This coordination is strongly non-linear and the different agent dynamics adapt to each other. The propagation front arises and locates in an area that occupies about 5-7 cell diameters ahead of the tumor main body, where we report an enhancement of the tumor activity and the integrin binding (see Fig 2).
We present here a novel mathematical model that covers specific dynamic aspects of tumor progression. From the modeling side, the novelties of our approach include the description of the dynamics of the tumor front and the link of cell membrane movements with the dynamics of some proteins, their concentrations, and their locations in the TM region. In this model, tumor density is governed by haptotactic and chemotactic processes, induced by MMP and active integrins, as well as by a flux saturated mechanism (see [17,47,48] and references therein). The latter allows the incorporation of biological features related to tumor dynamics (i.e., the viscosity of the medium and the speed of propagation) and the definition of sharp invasion profiles [47]. Therefore, this model acquires an excellent predictive capability, both quantitatively and qualitatively, since advanced mathematical concepts are confluent with precise biological experiments. Moreover, the outcomes of the model have partly guided the experiments and this proves the model consistency.

Localizing the front of GB
We focus on the dynamics of the active front in GB. Precisely, we analyze the distributions of MMP1, integrins, and FAK and their localization in relation to the tumor membrane density. We hypothesize a spatial heterogeneous activity of tumor cells, in terms of proteolysis and binding receptor activation. This spatial heterogeneity might determine the accumulation of MMP1 and the activation of FAK at the tumor front.
First, we show the presence of MMP1 protein at the active front of GB tissue. We dissected Drosophila brain samples with a genetically induced GB (repo>PI3K; EGFR). To delimitate the tumor front, we induced the co-expression of a membrane bound (myristoylated form) version of the red fluorescent protein (UAS-myrRFP) and, accordingly, all GB cell membranes were marked in red (Fig 2A 1 , red). MMP1 protein was detected by a specific monoclonal anti-MMP1 antibody (Fig 2A 2 , green). The confocal microscopy images show that the MMP1 signal is heterogeneous through the brain, and stronger in specific regions of the GB front, indicated with the white arrows in Fig 2A 1 . Then, we visualized focal adhesions with a specific monoclonal anti-FAK antibody (Fig 2B 2 , green). The confocal images show that FAK staining concentrates at the GB front, indicated with the white arrows in Fig 2B 1 . Finally, we co-stained for FAK (Fig 2C 1 , green) and MMP1 (Fig 2C 2 , magenta) proteins, and visualized the active front of GB cells (Fig 2C 1 , red inset). The confocal fluorescent images show that MMP1 and FAK signals are visible at the tumor front indicated with the white arrows in inset of Fig 2C 1 . These results are compatible with our hypothesis that MMP1 and FAK are characteristics features of the leading edge of migrating GB cells, with a central role in the process of GB expansion. We remark that in the figures visualizing the GB membrane, the healthy tissue is the only component not marked with the immunostaining. Thus it is represented by all the black regions in the images Fig 2A 1 and 2B 1 , and the inset of Fig 2C 1 .
Mathematically, we model the dynamics of the tumor front in a 1D scenario, i.e., we con- where h p is the maximum length of a microtube, while the functional for the tumor activity is given by Here, � indicates the convolution operator, I ½À h p ;h p � represents the identity function on the interval of semi-amplitude h p , and � and a F are parameters used for modulating the tumor activity.
These parameters are incorporated into the model on the base of the experimental results and this novel aspect enables for a better fitting of the biological data. With these basis, the dynamics of the tumor cells density N(t, x) is modeled as: The total flux of tumor cells (denoted by J N ðN; P; AÞ) is described by a combination of three main factors. First, the dynamics that J N exerts on cell movements include a saturated flow, which allows to define the movement of a sharp (non-diffusive) profile [47] and to incorporate the experimental data about the propagation front velocity and the porosity of the medium.
Then, J N collects information about the cell response to the gradient of soluble and insoluble components of the tumor microenvironment, precisely gradients of MMP1 P(t, x) and active integrins A(t, x). The former degrades the ECM, creating space for the tumor to migrate, while the latter is used to describe the migration toward a gradient of recognized adhesion sites. The proliferation term P N ðNÞ in (3) simply describes a logistic growth of glioma cells. However, several improvements of this growth term could be taken into account. Major details about the different terms and parameters involved in the glioma cell equation and about the possible improvements of the proliferation term are provided in the S1 Text. The combination of the different fluxes and the proliferation term leads to the following equation governing glioma cell dynamics: N m ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi À @ @x N a 1 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 þ @P @x � � 2 s @P @x þ N a 2 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Using this modeling approach, we demonstrate not only the role of proteases (in the specific, of the MMP1 family) and integrins in GB motility, but also their co-localization and spatial distribution with respect to the location of the tumor front. In the specific, this front region is characterized by a higher tumor membrane density compared to the bulk tumor (a region of higher tumor density).

MMP1s distribution
MMP1s facilitate the tumor invasion process by degrading the extracellular matrix. These proteins are produced by GB cells and released in the extracellular space, mostly in the area around the tumor invasion front. Here, TMs are present in a large number, thus, enhancing the proteolytic activity. Therefore, we propose the localization of MMP1s on TMs and the mathematical modeling of MMP1 dynamics in relation to the tumor front.
Experimentally, to analyze the protein distribution in the tumor we quantified the MMP1 signal in the inner GB mass and at the GB front. We visualized Drosophila brains with GB ( Fig  3A 1 and 3B 1 ) and immunostained with anti-MMP1 ( Fig 3A 2 and 3B 2 ). The results of our analysis indicate that, at the front region of GB (Fig 3A 1 ) where the membrane density (red curve in Fig 3A 3 ) is higher, MMP1s (green curve in Fig 3A 3 ) accumulate, showing a peak of concentration in the region corresponding to the TMs. Interestingly, this MMP1 maximum appears slightly shifted in the direction of GB migration with respect to the peak of GB membrane density. We analyzed the inner GB mass by taking a wider region of GB brains to measure MMP1 (Fig 3B 2 ) and GB density (Fig 3B 1 ). The confocal images and their analysis show homogeneous lower levels (basal levels) of GB membrane and MMP1 protein in the inner region (Fig 3B 3 ) compared to the front (Fig 3A 3 ). Therefore, these results suggest that MMP1 protein accumulation occurs specifically in the front region of GB, correlating with the peaks of GB cell membrane density in the TM region.
We model the dynamics and the distribution of MMP1s P(t, x) by taking into account three main phenomena: MMP1 production, diffusion in the extracellular space, and degradation. MMP1s are produced by tumor cells and this production is localized along the TM region, where the neoplastic tissue is in contact with the healthy tissue. The production of MMP1 is directly dependent on the heterogeneous tumor proteolytic activity. After MMP1 release into the extracellular space, their flux is limited in the areas surrounding the tumor mass, generating a sharp front ahead of the tumor. In particular, the protease flux is described using the same flux-saturated mechanism used for the tumor cells. In fact, saturated diffusion is not a mechanism exclusive of cells and it can affect species at different scales as long as the differences in size are taken into account in the estimation of velocities and viscosities (further details about tumor cells and proteases related parameters are provided in the S2 Text). Moreover, as the tumor front advances, the remaining MMP1s are degraded and maintain basal levels in the inner tumor region. The overall governing equation for MMP1s is given by: @P @t ¼ n P @ @x P m ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi MMP1 dynamics determines a dynamical degradation of the ECM E(t, x), described as We also include a basal level of ECM inside the main tumor mass, as some residual ECM material partially remains in the main tumor region after degradation. Further details about the derivation of the protease and ECM equations and the description of the involved parameters are given in the S1 Text.

Focal adhesions and integrin dynamics
Integrins are transmembrane receptors that allow cells to bind ECM ligand facilitating tumor cell movement. Upon activation, integrins mediate the organization of the cytoskeleton, cell cycle, and cell migration. In GB cells, the activation process occurs predominantly in the TM region, and active integrins result homogeneously distributed throughout TMs. Moreover, the conversion into the inactive (no bound) state occurs in the proximal region of the TMs with respect to the main GB mass, as shown by our experimental and mathematical results.
To determine the molecular changes at the GB front in relation to the activity of integrins, we immunostained GB brain samples with Talin, a mediator of integrin adhesivity [49], and with focal adhesion kinase, which is involved in signaling and cytoskeleton dynamics associated with integrin activity [50,51]. We analyzed confocal microscopy images of GB front regions, and compared fronts with low or high GB membrane signal. The quantification of the signals for anti-Talin and anti-FAK (Fig 4) shows that Talin (black curve in Fig 4A 4 and 4B 4 ) is reduced in the TM region, i.e., where the GB membrane signal is high (red curve in Fig 4A 4 and 4B 4 ), while FAK (magenta curve in Fig 4A 4 and 4B 4 ) is increased. Thus, the pattern of Talin and FAK expression is inverted at the front. Additionally, both the GB membrane and the FAK signals correlate. The GB/FAK signal correlation is maintained in different fronts, irrespective of their low (Fig 4A) or high (Fig 4B) levels of GB membrane signal. These results suggest an equivalent correlation for the reduction of Talin and the increase of FAK signals at the GB front, consistent with the relation of integrin dynamics and cell motility.
To confirm the functional contribution of integrins to GB progression and to validate our suggestions, we used specific RNAi constructs to knockdown myospheroid (mys), the Drosophila Integrin B subunit, or rhea, the Drosophila Talin, two key players for integrin function. The data show that GB cells require integrins to progress and expand (S1A-S1C Fig). Moreover, mys or rhea knockdown rescues the lethality caused by GB (S1D Fig). Finally, to confirm Talin and FAK inverse correlation, we compared inner GB mass and GB front areas. The data (S2 Fig) confirm that Talin and FAK maintain an inverse correlation, and they can be consider as indicators of the migratory status of the GB cells (see S3 Text for further details).
Mathematically, we split the integrin population into two subpopulations, referring to the active A(t, x) and the inactive state I(t, x) of the integrin receptors. Their dynamics take into account these processes of activation, and consequent inactivation, as well as a flux term describing the transport process of integrins due to the internal movement of GB cells.
Activation is mediated by the presence of the ECM, to which cells bind, and depends on the heterogeneous tumor activity along the TM region. The inactivation of integrins happens once the tumor has crawled on the ECM and moved forward. Since integrin receptors are located on the cell membrane and due to the movement of tumor cells, during the process of cell migration the receptors themselves are also transported, determining a flux of active and inactive integrins with estimated velocity (v Int ) depending on the tumor one. Precisely, since we model the dynamics of the bulk tumor, these transport terms allow us to translate the tumor dynamics to the front, where the interaction between integrins and the ECM occurs. Finally, assuming that, initially, the inactive integrins along the cell membrane have not reached their saturation value, a process of integrin production by exocytosis is also included in the model. The equations governing active and inactive integrins dynamics read Considering the experimental results shown in the S2 Fig, a basal level of active integrins in the tumor main body is included, as well as a basal level of inactive integrins along the TMs. A detailed description of these equations is provided in the S1 Text.

Motility features of GB cells at the front
The results shown in Fig 2, regarding the features that characterize the tumor front, together with the data on MMP1 and integrin dynamics, support our hypothesis about the processes that lead GB motility. Experimentally, we analyzed MMP1 and FAK concentration to prove their role as cues for cell motility, invasiveness, and migration [52,53]. We co-stained GB brain samples with anti-MMP1 and anti-FAK (Fig 5) and quantified the intensity of the fluorescent signals. We compared the signal in low and high tumor density fronts (Fig 5A and 5B, respectively). The results of the analysis are shown in Fig 5A 4 and 5B 4 and illustrate that, at the GB front ( Fig  5A 1 and 5B 1 ), the maximum signal for FAK ( Fig 5A 2 and 5B 2 ) occurs before the peak of MMP1 (Fig 5A 3 and 5B 3 ). This correlation occurs for both types of GB fronts, suggesting that this phenomenon is not a consequence of higher or lower GB membrane density. These results point towards a coordinated function between MMP1 activity and FAK dynamics and suggest that, at the front, FAK acts closer than MMP1 to the GB mass. By contrast, MMP1 plays its role further away from the front. Mechanistically, these results indicate that MMP1 activity in the proteolysis of the ECM is prior to the increase of integrin dynamics, it is associated with the presence of FAK, and, thus, together they contribute to the motility of GB cells.

Modeling results
From a mathematical viewpoint, we numerically solve the system given in Fig 6 with no flux boundary conditions. The results, represented in Fig 7, remarkably show how the model predicts that the region of greatest interest for the whole GB development process is the front, where TMs are located. The observed biological features correlate with the generation and progression of this critical tumor region, i.e., the front. Fig 7A shows the initial condition of the system, and it is accompanied by the quantifications of the basal levels of MMP1s with respect to the tumor density (Fig 7AE 1 ), and those of the active versus the inactive integrins (Fig 7AE 2 ). In both cases, the measurements are taken from areas inside the bulk tumor and the resulting data are incorporated into the model. Fig 7B and 7C show the time frames at 5 and 10 hours, respectively. The numerical simulations show that the model is capable of capture and predict the dynamics of the fronts and specific invasion profiles for each of the agents involved in the migration process. Precisely, Fig  7B and 7C illustrate how the model collects the exchange between active and inactive integrins at the beginning of the tumor front area, and how this corresponds to the experimental results (Fig 7BE). A plateau-like profile arises for A(t, x) in the TM region and, comparing Fig 7B 1 with Fig 7BE, these numerical results show a good agreement with the experimental data. Moreover, the model outputs regarding the MMP1 front profile correlate with the results obtained experimentally, as shown in the comparison between Fig 7C 1 and 7CE. MMP1 dynamics are characterized by an increased protein concentration along the TM region, indicating an enhanced tumor proteolytic activity, and by a steep profile ahead in the direction of tumor migration. Besides, we observe that the MMP1 external front is slightly shifted ahead with respect to the tumor front region. In fact, although MMP1s are produced along the protrusions, they are released in the extracellular space and spread in the areas around the tumor front. This behavior perfectly fits with our results.  The results also show a displacement of the MMP1 distribution with respect to active integrin and tumor membrane distributions. The model outcomes coincide with the experimental data in Fig 7DE. The previous considerations about the relative positions of these three elements (GB membrane density, FAK, and MMP1 distribution) might be seen as an indicator of the direction of tumor migration. In particular, peaks in the RFP distribution indicate areas with high tumor membrane density, i.e., the TM regions, while analogous peaks in the numerical simulations refer to high tumor density areas, i.e., the main tumor mass. The images B 1 ) and C 1 ) stand for a zoom of the indicated areas and they are correlated with the experimental data BE and CE, respectively. The parameters used for these simulations are listed in the S1 Table. https://doi.org/10.1371/journal.pcbi.1008632.g007

PLOS COMPUTATIONAL BIOLOGY
As we developed in the S2 Text, the MMP1 activity modifies the porosity of the medium, facilitating the spread and, therefore, the transport and progression of the tumor mass. In particular, tumor cell advance velocity might not be constant, but dependent on the porosity of the medium. This dependency could produce modifications of the tumor front, as shown in the S3 Fig. Furthermore, depending on the concentration of MMP1s, even in the case of constant tumor cell velocity, the front of the tumor may lose regularity and sometimes break into two different fronts. Then, if the tumor growth is not homogeneous and there is a certain heterogeneity in the growth of the invasion fronts, modifications and splittings of the frontal structure might arise, in agreement with experimental evidence (see S4 Fig. for further details).

Discussion
The location of the GB advancement front is essential to predict its evolutionary dynamics and it should help to design personalized therapies. In this work, we model and experimentally demonstrate that the active front of GB is a well-defined and consolidated structure. It harbors a considerable biochemical activity, defined by the protein ratio of integrins, proteases, and focal adhesions. Mathematically, the sharpness of the front results from a combination of the specific form of the flux term and the complex dynamical relationship between the agents involved in our system. This relationship also depends on the strength of the tactic mechanisms, due to integrin and protease activity and that drive GB migration. To ensure the robustness of the results, we analyze the possible influence of these tactic mechanisms on the tumor front profile as well as the effects of biomechanical changes in the porosity and stiffness of the tissue located closer to the TM area. The details of these additional studies are provided in the S1 Text and in the S3 and S4 Figs. These studies emphasize that, even if the tumor can break into separated masses, the profile steepness is preserved.
The mathematical model, based on the description of the specie dynamics involved in front progression, captures and reproduces the development of patterns that connect tumor, protease, and integrin evolution. These patterns characterize the dynamics of the GB front and have a predictive value on the directionality of the tumor migration. The mathematical model is based on a non-linear system of evolution equations (one for each of the involved agents) in which the mechanisms of production, chemotaxis, haptotaxis, and front dynamics compete with movement induced by the saturated flux in a porous media context. The contribution of mathematical models to tackle health problems is of great relevance. Usually, GB patients undergo surgery to remove the solid mass of the GB and a security resection margin, which can include a minimal amount of healthy tissue. However, the infiltrative nature of GB contributes to a possible re-appearance of the tumor. The results of the proposed multidisciplinary approach bring novelties that can help in the prediction of the behavior of the GB cells infiltrated in the brain. For instance, data from GB resection samples can be used to perform measurements of active/inactive integrins, MMP1s, and GB membrane density after surgery. In light of the model results on the relative location of these proteins, the outcomes of these measurements can contribute to predicting the invasive migratory behavior of GB cells, which may have relevant clinical consequences.
In conclusion, the feedback between experiments and modeling find common ground in our analysis to advance knowledge on the dynamics of the GB front. The results confirm that Drosophila model represents a suitable platform to analyze the molecular and cellular mechanisms implicated in GB progression. Our results point out the role of some of the main agents involved in the GB battlefront and the importance of incorporating this information into mathematical models in order to attempt to predict GB cell invasion dynamics.
Drosophila glioblastoma model. To reproduce glioblastoma in Drosophila we took advantage of the binary expression system Gal4/UAS [54]. We over-expressed constitutively active forms of EGFR and PI3K under the control of UAS sequences. To restrict this expression to glial cells, we used the specific enhancer repo (repo-Gal4>UAS-EGFRλ, UAS-dp110).
Imaging. Drosophila brain images were mounted in Vectashield mounting media with DAPI (Vector Laboratories) and analyzed by Confocal microscopy (LEICA TCS SP5) with a 63x oil immersion objective and 3x magnification. Images were processed using Image J 1.52t.
Quantifications. FAK, MMP1, and Talin signals were determined from images taken at the same confocal settings. Pixel intensity was measured using the plot profile tool from Fiji 1.52t.
Survival. Animal survival is represented as the percentage of flies as compared to control siblings that reach adulthood upon GB induction (repo-Gal4>UAS-EGFRλ, UAS-dp110). n>100.
Numerical solution of the model. The whole system of partial differential equations (see Fig 6), coupled with no-flux boundary conditions, was numerically solved with a self-developed code in Matlab (MathWorks Inc., Natick, MA). For the spatial discretization, we considered the Galerkin method on a spatial grid of 500 points. Precisely, the flux-limited terms were discretized using an IMEX version of the Galerkin scheme. For the time discretization, an implicit Euler scheme was used for proteases and tumor equations, while a fourth-order Runge-Kutta method for the other involved populations over a total of 15 � 10 6 time points.
Analysis of the biological data. For the analysis of the experimental data concerning the distributions of MMP1s, active and inactive integrins, we used a self-developed code in Matlab and the curve fitting toolbox, considering a weighted smoothing spline interpolation. In particular, for the analysis of MMP1 distribution, we analyzed 26 measurements on 8 different images taken from as many animals; then, for the distribution of Talin and FAK, we analyzed 14 measurements on 6 different images taken from as many animals; finally, for the combined distributions of FAK and MMP1, we analyzed 23 measurements on 13 different images taken from as many animals. The biological measurement data underlying the figures and the results presented in the text are provided in the S1 Data.  [8] (see Supplementary References in S2 Text), the minimum value for the velocity relates to a value of the porosity of 50%, while the maximum occurs around the value of 66%. The red curve shows how, as the velocity changes due to the ECM degradation process increase the medium porosity, cells closer to the front start moving faster than inner cells. This determines a heterogeneous modification of the invasion front, which slightly exceeds the homogenous front related to the constant velocity case. Eventually, the entire main tumor mass feels the changes in the velocity and a unique front is recovered. If there is heterogeneity in the growth of the front, the profile might not unify and a new front might emerge from this disturbance, as in the S4C Fig). (TIF)

S4 Fig. Effects of chemotaxis strength and heterogenous proliferation on tumor profile.
A), B) and C) show the dynamics of the tumor profile after 5 hours in different cases. In A), the effect of changes in the chemotactic sensitivity a 1 is presented. Since the proteolytic activity and, consequently, the concentration of MMP1 is enhanced in the front area, the stronger the parameter a 1 , the more evident the localization of the tactic effect. Tumor cells closer to the front acquire an increased overall velocity (due to both J flux-sat and J chemo fluxes) that leads to heterogenous fronts and, eventually, it might leads to a break of the tumor in two separated masses, as it can be observed in the experimental Fig D). In D), cell nuclei are marked in blue with DAPI, while GB cell membrane is marked with mystoylated-RFP in red. In particular, in D) higher intensity of the GB membrane marker indicates areas of tumor invasion. B) and C) show the comparison of our model with classical homogenous proliferation with the two possible models for heterogeneous proliferation, in the case of a 1 = 0.005. (TIF) S1 Data. Excel file with the original biological data. (XLSX)