Modeling invasion patterns in the glioblastoma battlefield

RESEARCH ARTICLE Modeling invasion patterns in the glioblastoma battlefield Martina ConteID 1 , Sergio Casas-TintòID 2 *, Juan SolerID 3 * 1BCAM - Basque Center for Applied Mathematics, Bilbao, Spain, 2Department of Molecular, Cellular and Developmental Neurobiology, Instituto Cajal CSIC, Madrid, Spain, 3Departamento de Matema ´tica Aplicada and Research Unit “Modeling Nature” (MNat), Universidad de Granada, Granada, Spain *[email protected] (SCT); [email protected] (JS) Abstract 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. Author summary Glioblastoma (GB) is a type of brain cancer that originated from glial cells. The infiltrative nature of GB cells is a key feature for understanding its aggressiveness and resistance to current treatments. Cellular protrusions, named as Tumor Microtubes (TMs) in GB, mediate the interaction between tumor and healthy tissue and the processes leading GB invasion. These protrusions are also responsible for several cell communication pathways (e.g. Hedgehog or WNT). We have developed a multidisciplinary approach, which combined biological biomarker measurements performed in Drosophila GB with a novel PLOS COMPUTATIONAL BIOLOGY PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008632 January 29, 2021 1 / 20 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 OPEN ACCESS Citation: Conte M, Casas-TintòS, Soler J (2021) Modeling invasion patterns in the glioblastoma battlefield. PLoS Comput Biol 17(1): e1008632. https://doi.org/10.1371/journal.pcbi.1008632 Editor: Philip K. Maini, Oxford, UNITED KINGDOM Received: August 6, 2020 Accepted: December 14, 2020 Published: January 29, 2021 Copyright: ©2021 Conte et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: All relevant data are within the manuscript and its Supporting information files. Funding: This work has been partially supported by the MINECO-Feder (Spain) research grant number RTI2018-098850-B-I00 (JS), the Junta de Andaluca (Spain) Project PY18-RT-2422 & A-FQM311-UGR18 (MC, SCT, JS), and by the MINECOFeder (Spain by the project PID2019-110116GBI00 (SCT). This research has been partially supported by the Basque Government through the BERC 20182021 program and by the Spanish State Research Agency through BCAM Severo Ochoa excellence accreditation SEV 2017-0718 mathematical model, to determine the interactions between proteases, integrins, and TM dynamics. The resulting model is able to predict the formation and infiltration of GB fronts, and, therefore, the directionality, aggressiveness, and progression of the tumor. 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–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. PLOS COMPUTATIONAL BIOLOGY Modeling invasion patterns in the glioblastoma battlefield PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008632 January 29, 2021 2 / 20 (MC). This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No. 713673 (MC). The project that gave rise to these results received the support of a fellowship from “la Caixa” Foundation (ID 100010434), fellowship code is LCF/BQ/IN17/11620056 (MC). This research was also funded by the Consejerı ´a de Economı ´a, Conocimiento, Empresas y Universidad and European Regional Development Fund (ERDF), ref. SOMM17/6109/UGR (JS). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist. 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–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–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. Fig 1. Diagram of GB and protein dynamics. Diagram of the interactions between the proteins involved in GB progression and that give rise to the mathematical model. A): GB cells produce and release in the extracellular space Matrix Metalloproteases (MMPs), which proteolyze the Extra Cellular Matrix (ECM) components. B): Magnification of Tumor Microtubes (TMs). Integrins are activated in the GB tumor microtubes upon interaction with ECM proteins. Active integrins, interacting with Actin filaments and the Talin adaptor protein, activate the Focal Adhesion Kinase (FAK) protein to promote cytoplasm dynamics. https://doi.org/10.1371/journal.pcbi.1008632.g001 PLOS COMPUTATIONAL BIOLOGY Modeling invasion patterns in the glioblastoma battlefield PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008632 January 29, 2021 3 / 20 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 antiTalin 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 Fig 2. MMP1 and FAK proteins in GB. Fluorescent confocal images of Drosophila 3rd instar larvae brain with GB. A) 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 Modeling invasion patterns in the glioblastoma battlefield PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008632 January 29, 2021 4 / 20 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–40]. The process of growth and retraction of TMs plays a fundamental role in cell-to-cell communication routing [17,18, 20–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, PLOS COMPUTATIONAL BIOLOGY Modeling invasion patterns in the glioblastoma battlefield PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008632 January 29, 2021 5 / 20 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. Results 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 antiMMP1 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 consider (t,x)2[0, T]×O, with T>0 and O¼ ½0;bO� � R. First, we characterize the TM region (L TM ) and the heterogeneity of the tumor activity in the tumor domain using the functional FðNÞ. Defined the tumor support Sup(N) = [0, b N ]�O, the TM region is described as LTM ≔fx2O:9a2 ½0;hp�such that xa¼bNg;ð1Þ where h p is the maximum length of a microtube, while the functional for the tumor activity is given by FðNÞ≔1 ðN�I½ hp;hp�þ�ÞaF:ð2Þ Here, �indicates the convolution operator, I½ hp;hp�represents the identity function on the interval of semi-amplitude h p , and �and aFare parameters used for modulating the tumor activity. PLOS COMPUTATIONAL BIOLOGY Modeling invasion patterns in the glioblastoma battlefield PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008632 January 29, 2021 6 / 20 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: @N @t¼  @ @xJNðN;P;AÞþPNðNÞ:ð3Þ The total flux of tumor cells (denoted by JNðN;P;AÞ) is described by a combination of three main factors. First, the dynamics that JNexerts 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, JNcollects 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 PNð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 @t¼nN @ @x Nm ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N2þnN vN � �2@N @x �������� 2 s@N @x 0 B B B B @ 1 C C C C A @ @xNa1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1þ@P @x � �2 s@P @xþNa2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1þ@A @x � �2 s@A @x 0 B B B B @ 1 C C C C A þa3N1N KN � � ð4Þ 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. PLOS COMPUTATIONAL BIOLOGY Modeling invasion patterns in the glioblastoma battlefield PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008632 January 29, 2021 7 / 20 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 Fig 3. MMP1 accumulation in GB front and inner GB mass. Fluorescent confocal images of Drosophila 3rd instar larvae brain with GB marked in red (the white line in A 1 and B 1 indicated the location of the measurement at the front and in the inner GB mass, respectively), and stained with anti-MMP1 in green (A 2 and B 2 ). In A 3 and B 3 , we show the quantification and graphical representations of the fluorescent intensity for GB and MMP1 signals along the white lines in A 1 ,A 2 ,B 1 and B 2 . Dots represent the data and lines represent the fitting. The healthy tissue is the only component not marked in the immunostaining, thus it is represented by all the black regions in the images A 1 and B 1 . https://doi.org/10.1371/journal.pcbi.1008632.g003 PLOS COMPUTATIONAL BIOLOGY Modeling invasion patterns in the glioblastoma battlefield PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008632 January 29, 2021 8 / 20 levels in the inner tumor region. The overall governing equation for MMP1s is given by: @P @t¼nP @ @x Pm ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P2þnP vP � �2@P @x �������� 2 s@P @x 0 B B B B @ 1 C C C C Aþa4EFðNÞwLTM a5P N :ð5Þ MMP1 dynamics determines a dynamical degradation of the ECM E(t,x), described as @E @t ¼  a6E P :ð6Þ 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. PLOS COMPUTATIONAL BIOLOGY Modeling invasion patterns in the glioblastoma battlefield PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008632 January 29, 2021 9 / 20 S3 Text. Additional comments on the S1–S4 Figs. (PDF) S1 Table. Parameter estimation. (PDF) S1 Fig. Functional integrins are required for GB progression. Low magnification confocal images of Drosophila GB larvae brain in A) and rhea knockdown (UAS-rhea RNAi) in B), or mys knockdown (UAS-mys RNAi) in C). GB cell membrane is marked with myristoylatedRFP (red) and nuclei are marked with DAPI (blue). The images show that rhea or mys knockdown prevents the expansion of GB and rescues the lethality (percentage of survival in D). (TIF) S2 Fig. Talin and FAK dynamics: Comparison between the inner GB mass and the GB front. Fluorescent confocal images of Drosophila 3rd instar larvae brain with GB marked in red (A 1 and B 1 ), and stained with anti-Talin (green in A 2 and B 2 ) and anti-FAK (magenta in A 3 and B 3 ). In A 4 and B 4 , we show the quantification of the fluorescent signals and the graphical representation of the fluorescent intensity for GB, Talin and FAK signals along the white lines in A 1 -A 3 and B 1 -B 3 that indicate the location of the measurements. Dots represent the data and lines represent the fitting. The healthy tissue is the only component not marked in the immunostaining, thus it is represented by all the black regions in the images A 1 and B 1 . (TIF) S3 Fig. Effects of porosity changes on tumor profile. A), B), and C) show the comparison between the tumor density profile in the case of flux saturated model with constant velocity v N (in black) and with v N =v N (�) (in red) at the time steps referring to 3, 6 and 9 hours of tumor evolution, respectively. For this comparison we consider a tumor cell equation given by @N @t¼  Jflux-satðNÞ. In �/v), the profile of v N =v N (�) is shown. Specifically, accordingly with [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 Jflux-sat and Jchemo 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) PLOS COMPUTATIONAL BIOLOGY Modeling invasion patterns in the glioblastoma battlefield PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008632 January 29, 2021 16 / 20 S1 Data. Excel file with the original biological data. (XLSX) Acknowledgments The authors thank Isabel Guerrero for the reagents, her valuable advice during the development of this work, and for comments on the manuscript. The authors thank also Professors Luca Gerardo-Giorda, Jose ´L. Lo ´pez, and Juan J. Nieto for their helpful discussions. Author Contributions Conceptualization: Martina Conte, Sergio Casas-Tintò, Juan Soler. Data curation: Martina Conte, Sergio Casas-Tintò, Juan Soler. Formal analysis: Martina Conte, Sergio Casas-Tintò, Juan Soler. Funding acquisition: Martina Conte, Sergio Casas-Tintò, Juan Soler. Investigation: Martina Conte, Sergio Casas-Tintò, Juan Soler. Methodology: Martina Conte, Sergio Casas-Tintò, Juan Soler. Project administration: Martina Conte, Sergio Casas-Tintò, Juan Soler. Resources: Martina Conte, Sergio Casas-Tintò. Software: Martina Conte, Sergio Casas-Tintò, Juan Soler. Supervision: Juan Soler. Validation: Martina Conte, Sergio Casas-Tintò, Juan Soler. Visualization: Martina Conte, Sergio Casas-Tintò, Juan Soler. Writing – original draft: Martina Conte, Sergio Casas-Tintò, Juan Soler. Writing – review & editing: Martina Conte, Sergio Casas-Tintò, Juan Soler. References 1. Holland EC. Glioblastoma multiforme: the terminator. Proceedings of the National Academy of Sciences. 2000; 97(12):6242–6244. https://doi.org/10.1073/pnas.97.12.6242 PMID: 10841526 2. Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie SM, Wakimoto H, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science. 2014; 344(6190):1396–1401. https://doi.org/10.1126/science.1254257 PMID: 24925914 3. Johnson BE, Mazor T, Hong C, Barnes M, Aihara K, McLean CY, et al. Mutational analysis reveals the origin and therapy-driven evolution of recurrent glioma. Science. 2014; 343(6167):189–193. https://doi. org/10.1126/science.1239947 PMID: 24336570 4. Barthel FP, Johnson KC, Varn FS, Moskalik AD, Tanner G, Kocakavuk E, et al. Longitudinal molecular trajectories of diffuse glioma in adults. Nature. 2019; 576(7785):112–120. https://doi.org/10.1038/ s41586-019-1775-1 PMID: 31748746 5. Quail DF, Joyce JA. The microenvironmental landscape of brain tumors. Cancer cell. 2017; 31(3):326– 341. https://doi.org/10.1016/j.ccell.2017.02.009 6. Picariello HS, Kenchappa RS, Rai V, Crish JF, Dovas A, Pogoda K, et al. Myosin IIA suppresses glioblastoma development in a mechanically sensitive manner. Proceedings of the National Academy of Sciences. 2019; 116(31):15550–15559. https://doi.org/10.1073/pnas.1902847116 7. Osswald M, Jung E, Sahm F, Solecki G, Venkataramani V, Blaes J, et al. Brain tumour cells interconnect to a functional and resistant network. Nature. 2015; 528(7580):93–98. https://doi.org/10.1038/ nature16071 PMID: 26536111 PLOS COMPUTATIONAL BIOLOGY Modeling invasion patterns in the glioblastoma battlefield PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008632 January 29, 2021 17 / 20 8. Venkatesh HS, Morishita W, Geraghty AC, Silverbush D, Gillespie SM, Arzt M, et al. Electrical and synaptic integration of glioma into neural circuits. Nature. 2019; 573(7775):539–545. https://doi.org/10. 1038/s41586-019-1563-y PMID: 31534222 9. Venkataramani V, Tanev DI, Strahle C, Studier–Fischer A, Fankhauser L, Kessler T, et al. Glutamatergic synaptic input to glioma cells drives brain tumour progression. Nature. 2019; 573(7775):532–538. https://doi.org/10.1038/s41586-019-1564-x PMID: 31534219 10. Puchalski RB, Shah N, Miller J, Dalley R, Nomura SR, Yoon JG, et al. An anatomic transcriptional atlas of human glioblastoma. Science. 2018; 360(6389):660–663. https://doi.org/10.1126/science.aaf2666 PMID: 29748285 11. Read RD, Cavenee WK, Furnari FB, Thomas JB. A drosophila model for EGFR-Ras and PI3K-dependent human glioma. PLoS Genetics. 2009; 5(2):e1000374. https://doi.org/10.1371/journal.pgen. 1000374 12. Portela M, Venkataramani V, Fahey-Lozano N, Seco E, Losada–Perez M, Winkler F, et al. Glioblastoma cells vampirize WNT from neurons and trigger a JNK/MMP signaling loop that enhances glioblastoma progression and neurodegeneration. PLoS biology. 2019; 17(12):e3000545. https://doi.org/10. 1371/journal.pbio.3000545 PMID: 31846454 13. Portela M, SeguraCollar B, Argudo I, Sa ´iz A, Gargini R, Sa ´nchez-Gòmez P, et al. Oncogenic dependence of glioma cells on kish/TMEM167A regulation of vesicular trafficking. Glia. 2019; 67(2):404–417. https://doi.org/10.1002/glia.23551 PMID: 30506943 14. Chen X, Wanggou S, Bodalia A, Zhu M, Dong W, Fan JJ, et al. A feedforward mechanism mediated by mechanosensitive ion channel PIEZO1 and tissue mechanics promotes glioma aggression. Neuron. 2018; 100(4):799–815. https://doi.org/10.1016/j.neuron.2018.09.046 PMID: 30344046 15. Jacquemet G, Hamidi H, Ivaska J. Filopodia in cell adhesion, 3D migration and cancer cell invasion. Current opinion in cell biology. 2015; 36:23–31. https://doi.org/10.1016/j.ceb.2015.06.007 16. Ramı ´rez–Weber FA, Kornberg TB. Cytonemes: cellular processes that project to the principal signaling center in Drosophila imaginal discs. Cell. 1999; 97(5):599–607. https://doi.org/10.1016/S0092-8674 (00)80771-0 17. Verbeni M, Sa ´nchez O, Mollica E, Siegl–Cachedenier I, Carleton A, Guerrero I, et al. Morphogenetic action through flux-limited spreading. Physics of Life Reviews. 2013; 10(4):457–475. https://doi.org/10. 1016/j.plrev.2013.10.005 PMID: 23831049 18. Gonza ´lez–Me ´ndez L, Seijo-Barandiara ´n I, Guerrero I. Cytoneme-mediated cell-cell contacts for Hedgehog reception. Elife. 2017; 6:e24045. https://doi.org/10.7554/eLife.24045 19. Gonza ´lez–Me ´ndez L, Gradilla AC, Guerrero I. The cytoneme connection: direct long-distance signal transfer during development. Development. 2019; 146(9):dev174607. https://doi.org/10.1242/dev. 174607 20. Cardozo MJ, Sa ´nchez–Arrones L, Sandonis A, Sa ´nchez–Camacho C, Gestri G, Wilson SW, et al. Cdon acts as a Hedgehog decoy receptor during proximal-distal patterning of the optic vesicle. Nature communications. 2014; 5(1):1–3. https://doi.org/10.1038/ncomms5272 PMID: 25001599 21. Sanders TA, Llagostera E, Barna M. Specialized filopodia direct long-range transport of SHH during vertebrate tissue patterning. Nature. 2013; 497(7451):628–32. https://doi.org/10.1038/nature12157 22. Huang H, Liu S, Kornberg TB. Glutamate signaling at cytoneme synapses. Science. 2019; 363 (6430):948–955. https://doi.org/10.1126/science.aat5053 23. Kornberg TB, Gilboa L. Nanotubes in the niche. Nature. 2015; 523(7560):292–294. https://doi.org/10. 1038/nature14631 24. Rosenbauer J, Zhang C, Mattes B, Reinartz I, Wedgwood K, Schindler S, et al. Modeling of Wnt-mediated tissue patterning in vertebrate embryogenesis. PLOS Computational Biology. 2020; 16(6): e1007417. https://doi.org/10.1371/journal.pcbi.1007417 PMID: 32579554 25. Junyent S, Garcin CL, Szczerkowski JL, Trieu TJ, Reeves J, Habib SJ. Specialized cytonemes induce self-organization of stem cells. Proceedings of the National Academy of Sciences. 2020; 117(13):7236– 7244. https://doi.org/10.1073/pnas.1920837117 26. Mattes B, Dang Y, Greicius G, Kaufmann LT, Prunsche B, Rosenbauer J, et al. Wnt/PCP controls spreading of Wnt/β-catenin signals by cytonemes in vertebrates. Elife. 2018; 7:e36953. https://doi.org/ 10.7554/eLife.36953 PMID: 30060804 27. Kress H, Stelzer EH, Holzer D, Buss F, Griffiths G, Rohrbach A. Filopodia act as phagocytic tentacles and pull with discrete steps and a load-dependent velocity. Proceedings of the National Academy of Sciences. 2007; 104(28):11633–11638. https://doi.org/10.1073/pnas.0702449104 28. Tamkun JW, DeSimone DW, Fonda D, Patel RS, Buck C, Horwitz AF, et al. Structure of integrin, a glycoprotein involved in the transmembrane linkage between fibronectin and actin. Cell. 1986; 46(2):271– 282. https://doi.org/10.1016/0092-8674(86)90744-0 PMID: 3487386 PLOS COMPUTATIONAL BIOLOGY Modeling invasion patterns in the glioblastoma battlefield PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008632 January 29, 2021 18 / 20 29. Brown NH, Gregory SL, Rickoll WL, Fessler LI, Prout M, White RA, et al. Talin is essential for integrin function in Drosophila. Developmental cell. 2002; 3(4):569–579. https://doi.org/10.1016/S1534-5807 (02)00290-3 PMID: 12408808 30. Ellert–Miklaszewska A, Poleszak K, Pasierbinska M, Kaminska B. Integrin Signaling in Glioma Pathogenesis: From Biology to Therapy. International Journal of Molecular Sciences. 2020; 21(3):888. https://doi.org/10.3390/ijms21030888 31. Huttenlocher A, Horwitz AR. Integrins in cell migration. Cold Spring Harbor perspectives in biology. 2011; 3(9):a005074. 32. Demuth T, Berens ME. Molecular mechanisms of glioma cell migration and invasion. Journal of neurooncology. 2004; 70(2):217–228. https://doi.org/10.1007/s11060-004-2751-6 33. Yu Q, Stamenkovic I. Cell surface-localized matrix metalloproteinase-9 proteolytically activates TGF-β and promotes tumor invasion and angiogenesis. Genes & development. 2000; 14(2):163–176. 34. Rao JS. Molecular mechanisms of glioma invasiveness: the role of proteases. Nature Reviews Cancer. 2003; 3(7):489–501. https://doi.org/10.1038/nrc1121 35. Egeblad M, Werb Z. New functions for the matrix metalloproteinases in cancer progression. Nature reviews cancer. 2002; 2(3):161–174. https://doi.org/10.1038/nrc745 36. Matı ´as–Roma ´n S, Ga ´lvez BG, Genı ´s L, Ya ´ñez-Mo ´M, de la Rosa G, Sa ´nchez–Mateos P, et al. Membrane type 1–matrix metalloproteinase is involved in migration of human monocytes and is regulated through their interaction with fibronectin or endothelium. Blood. 2005; 105(10):3956–3964. https://doi. org/10.1182/blood-2004-06-2382 PMID: 15665118 37. Ogier C, Bernard A, Chollet AM, Le Diguardher T, Hanessian S, Charton G, et al. Matrix metalloproteinase2 (MMP2) regulates astrocyte motility in connection with the actin cytoskeleton and integrins. Glia. 2006; 54(4):272–284. https://doi.org/10.1002/glia.20349 PMID: 16845676 38. Berges R, Denicolai E, Tchoghandjian A, Baeza–Kallee N, Honore S, Figarella–Branger D, et al. Proscillaridin A exerts anti-tumor effects through GSK3βactivation and alteration of microtubule dynamics in glioblastoma. Cell death & disease. 2018; 9(10):1–14. https://doi.org/10.1038/s41419-018-1018-7 PMID: 30250248 39. Berges R, Tchoghandjian A, Honore S, Estève MA, Figarella–Branger D, Bachmann F, et al. The novel tubulin-binding checkpoint activator BAL101553 inhibits EB1-dependent migration and invasion and promotes differentiation of glioblastoma stem-like cells. Molecular Cancer Therapeutics. 2016; 15 (11):2740–2749. https://doi.org/10.1158/1535-7163.MCT-16-0252 PMID: 27540016 40. Berges R, Baeza–Kallee N, Tabouret E, Chinot O, Petit M, Kruczynski A, et al. End-binding 1 protein overexpression correlates with glioblastoma progression and sensitizes to Vinca-alkaloids in vitro and in vivo. Oncotarget. 2014; 5(24):12769. https://doi.org/10.18632/oncotarget.2646 PMID: 25473893 41. Galjart N. Plus-end-tracking proteins and their interactions at microtubule ends. Current Biology. 2010; 20(12):R528–R537. https://doi.org/10.1016/j.cub.2010.05.022 42. Cherry AE, Vicente JJ, Xu C, Morrison RS, Ong SE, Wordeman L, et al. GPR124 regulates microtubule assembly, mitotic progression, and glioblastoma cell proliferation. Glia. 2019; 67(8):1558–1570. https:// doi.org/10.1002/glia.23628 PMID: 31058365 43. Kowalski–Chauvel A, Modesto A, Gouaze–Andersson V, Baricault L, Gilhodes J, Delmas C, et al. Alpha-6 integrin promotes radioresistance of glioblastoma by modulating DNA damage response and the transcription factor Zeb1. Cell death & disease. 2018; 9(9):1–12. https://doi.org/10.1038/s41419018-0853-x PMID: 30158599 44. Hoshino A, Costa–Silva B, Shen TL, Rodrigues G, Hashimoto A, Mark MT, et al. Tumour exosome integrins determine organotropic metastasis. Nature. 2015; 527(7578):329–335. https://doi.org/10. 1038/nature15756 PMID: 26524530 45. Mallet DG, Pettet GJ. A mathematical model of integrin-mediated haptotactic cell migration. Bulletin of mathematical biology. 2006; 68(2):231. https://doi.org/10.1007/s11538-005-9032-1 46. Kim Y, Jeon H, Othmer H. The role of the tumor microenvironment in glioblastoma: A mathematical model. IEEE Transactions on Biomedical Engineering. 2016; 64(3):519–527. 47. Calvo J, Campos J, Caselles V, Sa ´nchez O, Soler J. Pattern formation in a flux limited reaction–diffusion equation of porous media type. Inventiones mathematicae. 2016; 206(1):57–108. https://doi.org/ 10.1007/s00222-016-0649-5 48. Calvo J, Campos J, Caselles V, Sa ´nchez O, Soler J. Flux-saturated porous media equations and applications. EMS Surveys in Mathematical Sciences. 2015; 2(1):131–218. https://doi.org/10.4171/EMSS/ 11 49. Klapholz B, Herbert SL, Wellmann J, Johnson R, Parsons M, Brown NH. Alternative mechanisms for talin to mediate integrin function. Current biology. 2015; 25(7):847–857. https://doi.org/10.1016/j.cub. 2015.01.043 PLOS COMPUTATIONAL BIOLOGY Modeling invasion patterns in the glioblastoma battlefield PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008632 January 29, 2021 19 / 20 50. Macagno JP, Vera JD, Yu Y, MacPherson I, Sandilands E, Palmer R, et al. FAK acts as a suppressor of RTK-MAP kinase signalling in Drosophila melanogaster epithelia and human cancer cells. PLoS Genetocs. 2014; 10(3):e1004262. https://doi.org/10.1371/journal.pgen.1004262 PMID: 24676055 51. Fox GL, Rebay I, Hynes RO. Expression of DFak56, a Drosophila homolog of vertebrate focal adhesion kinase, supports a role in cell migration in vivo. Proceedings of the National Academy of Sciences. 1999; 96(26):14978–14983. https://doi.org/10.1073/pnas.96.26.14978 52. Hsia DA, Mitra SK, Hauck CR, Streblow DN, Nelson JA, Ilic D, et al. Differential regulation of cell motility and invasion by FAK. The Journal of cell biology. 2003; 160(5):753–767. https://doi.org/10.1083/jcb. 200212114 PMID: 12615911 53. Das A, Monteiro M, Barai A, Kumar S, Sen S. MMP proteolytic activity regulates cancer invasiveness by modulating integrins. Scientific reports. 2017; 7(1):1–13. 54. Brand AH, Perrimon N. Targeted gene expression as a means of altering cell fates and generating dominant phenotypes. development. 1993; 118(2):401–415. PLOS COMPUTATIONAL BIOLOGY Modeling invasion patterns in the glioblastoma battlefield PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1008632 January 29, 2021 20 / 20