Role of extracellular matrix and microenvironment in regulation of tumor growth and LAR-mediated invasion in glioblastoma

The cellular dispersion and therapeutic control of glioblastoma, the most aggressive type of primary brain cancer, depends critically on the migration patterns after surgery and intracellular responses of the individual cancer cells in response to external biochemical cues in the microenvironment. Recent studies have shown that miR-451 regulates downstream molecules including AMPK/CAB39/MARK and mTOR to determine the balance between rapid proliferation and invasion in response to metabolic stress in the harsh tumor microenvironment. Surgical removal of the main tumor is inevitably followed by recurrence of the tumor due to inaccessibility of dispersed tumor cells in normal brain tissue. In order to address this complex process of cell proliferation and invasion and its response to conventional treatment, we propose a mathematical model that analyzes the intracellular dynamics of the miR-451-AMPK- mTOR-cell cycle signaling pathway within a cell. The model identifies a key mechanism underlying the molecular switches between proliferative phase and migratory phase in response to metabolic stress in response to fluctuating glucose levels. We show how up- or down-regulation of components in these pathways affects the key cellular decision to infiltrate or proliferate in a complex microenvironment in the absence and presence of time delays and stochastic noise. Glycosylated chondroitin sulfate proteoglycans (CSPGs), a major component of the extracellular matrix (ECM) in the brain, contribute to the physical structure of the local brain microenvironment but also induce or inhibit glioma invasion by regulating the dynamics of the CSPG receptor LAR as well as the spatiotemporal activation status of resident astrocytes and tumor-associated microglia. Using a multi-scale mathematical model, we investigate a CSPG-induced switch between invasive and non-invasive tumors through the coordination of ECM-cell adhesion and dynamic changes in stromal cells. We show that the CSPG-rich microenvironment is associated with non-invasive tumor lesions through LAR-CSGAG binding while the absence of glycosylated CSPGs induce the critical glioma invasion. We illustrate how high molecular weight CSPGs can regulate the exodus of local reactive astrocytes from the main tumor lesion, leading to encapsulation of non-invasive tumor and inhibition of tumor invasion. These different CSPG conditions also change the spatial profiles of ramified and activated microglia. The complex distribution of CSPGs in the tumor microenvironment can determine the nonlinear invasion behaviors of glioma cells, which suggests the need for careful therapeutic strategies.

where Λ g is the signaling rate from glucose to miR-451, g is the glucose level, s 1 , s 2 are the signalling pathways to AMPK complex and mTOR, respectively, Λ 1 , Λ 3 , Λ 7 are the autocatalytic enhancement parameters for miR-451, AMPK complex, mTOR, respectively and Λ 2 , Λ 4 , Λ 8 are the Hill-type inhibition saturation parameters from the counter part of miR-451, AMPK complex, mTOR, respectively, Λ 5 is the inhibition strength of miR-451 by the AMPK complex, Λ 6 is the inhibition strength of the AMPK complex by miR-451, Λ 9 is the inhibition strength of the mTOR by the AMPK complex, and finally µ 1 , µ 2 , µ 3 are microRNA/protein degradation rates of miR-451, AMPK complex, and mTOR, respectively. As indicated in Eq (1), the signal g increases the rate of miR-451 activation through the function f (g) in the first term, while AMPK-dependent inhibition of miR-451 is expressed through the function F (a) in the denominator of the second term. A simple requirement of these functions is that ∂f ∂g > 0 for all non-negative g, and ∂F ∂a > 0 for all non-negative a. Similarly, the first term h 1 (s 1 ) on the right-hand side of Eq (2) represents an increase in AMPK activity induced by the signal s, and miR-451-dependent suppression of AMPK activity is explored through the function H 1 (m) in the denominator. In the same manner, the first term h 2 (s 2 ) on the right-hand side of Eq (3) represents an increase in mTOR activity induced by the signal s 2 , and AMPK-dependent suppression of mTOR activity is accounted for through the function H 2 (a) in the denominator. One must also have ∂h1 ∂s1 > 0, ∂h2 ∂s2 > 0, ∂H1 ∂m > 0, ∂H2 ∂a > 0 for all non-negative s 1 , s 2 , m, a. Based on biological observations (Fig 1A), we assume that f (g) = Λ g g, F (a) = a 2 , h 1 (s 1 ) = s 1 , H 1 (m) = m 2 , h 2 (s 2 ) = s 2 , H 2 (a) = a 2 .
Furthermore, performing now the following non-dimensionalization from the dimensional Eqs (1)-(4) we obtain the associated dimensionless equations with a smaller set of essential control parameters, namely: 2 dR dt = S 2 + λ 5 λ 2 6 λ 2 6 + γA 2 − R. (8) Parameter estimation miRNAs are typically more stable than their targets [2,3] and the typical half-life of a miRNA is much larger than that of AMPK [4,5] or mTOR, leading to the small values of the fractions 1 = µ1 µ2 = 0.02, 2 = µ1 µ3 = 0.02 [6][7][8]. From the estimated concentrations of miRNAs (80 pM − 2.2 µM ) in an animal cell (assuming 1,000-25,000 µm 3 volume) [9], we take m * = 1.0 µM . The normal (high; 4.5 g/l) and low (0.3 g/l) glucose levels were used in [1] in order to test the effect of glucose on the miR451 expression and AMPK activities, and proliferation/migration patterns of glioma cells. Using this range of glucose levels (0.3 − 4.5 g/l) in [1] and the miR-451 reference value (m * ), we determined the glucose signaling strength, (2.4 × 10 −5 − 2.4 × 10 −3 ) µM/h, resulting in a range of dimensionless glucose input level λ g G = Λgg µ1m * = 0.01 − 1.0. Using the reference value of AMPK, a * = 100nM taken from the measured values (35-150 nM ) in rat liver [10], and taking the signal source of the AMPK complex, s 1 = 2.4 nM/h, we estimate the dimensionless value of signaling strength to the AMPK complex to be S 1 = s1 µ2a * = 0.2 [6,7]. By observing the same patterns of up-and down-regulation of mTOR as those of miR-451 in response to various glucose levels, we take slightly larger values of S 2 = s2 µ3r * = 1.2 in the work. We take the dimensionless parameter value λ 1 = Λ1 µ1m * = 4.0 by assuming that the autocatalytic rate (Λ 1 ) of miR-451 is at least 4-fold larger than its negative contribution (µ 1 m * ) from its decay/consumption in the absence of the inhibitory pathway from the AMPK complex [6,7]. Relying on the observed mutual antagonism between miR-451 and AMPK complex, we also take λ 3 = Λ3 µ2a * = 4.0 in a similar fashion. The Hill-type coefficients Λ 2 , Λ 4 , Λ 6 (and their corresponding dimensionless parameters λ 2 , λ 4 , λ 6 ) are fixed to be unity [6,7]. Once the parameters above were determined, we fitted the data in the level of LKB1/AMPK activity in [1] in response to negative control and over expressed levels of miR-451 (Fig 5B in [1]) in the following way: The steady state of AMPK complex (A s ) in terms of miR-451 level (M s ) and other parameters in the miR-451-AMPK model can be rewritten by The up-regulated AMPK complex level (∼500 pmol of phosphate incorporated in a dimensional form) for negative control of miR-451 was down-regulated (∼100 pmol of phosphate incorporated in a dimensional form). Using Eq (9) above, we estimated the parameter value of β to be 1.0 which gives the low AMPK complex activity (∼ A s = 0.9) in response to a high dimensionless miR-451 level (M s = 4.2) and the high AMPK complex activity (A s = 4.2) in response to negative control of the miR-451 level (M s = 0.0), resulting in a reasonable ∼4.7-fold difference in the AMPK complex activities as in the experiments in [1]. Finally, by observing the behavior of the system and experimental data, we assume that the inhibition strength (α = 1.6) of miR-451 by the AMPK complex is slightly larger than the inhibition strength (β = 1.0) of the AMPK complex by miR-451 [6,7] while keeping the same value for the inhibition strength (γ = 1.0) of mTOR by the AMPK complex.
Analysis of the core control system Characterization of migration and proliferation under perturbation of key parameters In Fig S1, we show sensitivity of variables in the core control system in response to glucose levels (G) to changes in inhibition strength of miR-451 by AMPK complex (α=1.0, 1.6, 2.0, 2.5). As α is increased, the bifurcation curves of all variables (miR-451, AMPK, mTOR) shift to the right (Figs S1A-S1C). Overall, this increase in α also moves the bi-stability window (W b ) to the right in all cases of miR-451, AMPK, and mTOR and leads to a decrease in the size of the window (|W b |), which is also illustrated in Fig S1D. This indicates that the probability of switching to the migration phase is increased when the inhibition strength of miR-451 by AMPK (α) is increased. For example, the cell is already in the migratory phase (M < th M , A > th A , R < th R ) when G = 0.7 in the case of higher α (α = 2.5) while it is still in the proliferative phase (M > th M , A < th A , R > th R ) in the base case (α * = 1.6). On the other hand, glioma cells would grow even under glucose withdrawal conditions when α is decreased. For example, the cell is still in the proliferative phase (M > th M , A < th A , R > th R ) for a low glucose level (G = 0.2) when α is lowered (α = 1.0) while it should be in the migratory phase (M < th M , A > th A , R < th R ) in the base case (α * = 1.6). Fig S1D shows a phase diagram of G − α for switching behaviors including a monostable region and bistable and one-way switches. An increase in the inhibition strength of miR-451 by AMPK (α) induces a transition from a one-way switch to a bistable switch and to a mono-stability. A decrease in α increases the size of the bi-stability window (|W b |). We recall that, for the base case (α * = 1.6), the existence of bi-stability regimes illustrates that glucose can regulate the forward transition from the migratory phenotype (T m ) to the proliferative phenotype (T p ) and the recurrent process from the T pstate to the T m -state. However, when miR-451 activity is highly enhanced with the decreased inhibition signal from AMPK complex (a decrease in α), glucose may regulate only forward transition from T m -state to T p -state (one way switch) and a monotonic decrease in glucose does not push cells to the migratory Other parameters are fixed as in Table 1. The base case (α * = 1.6) was marked in star in (A-C). (E) Phenotypic changes in response to glucose levels (G) when the key inhibition parameters (α) is varied.
phase. Therefore, this mechanism may provide a way of keeping glioma cells in a proliferative phase for a certain period of time in order to prevent the invasion process for second surgery [7]. The phenotypic changes between T m -state and T p -state for increased or decreased α is illustrated in Fig S1E. An increase (or decrease) in α generates higher probability of switching to the migratory (or proliferative) phenotype for given glucose level by moving the bifurcation curve and bistable window (W b ) to the right or left. In Fig S2(A-D) we show the effect of the inhibition strength of AMPK complex by miR-451 (β = 0.07, 1.0, 2.5) on the dynamics of the core control system. As β is increased, the bifurcation curves of all variables (miR-451, AMPK, mTOR) shift to the left (Figs S2(A-C)). In contrast to the case of changes in α, this increase in β reverses its shifting direction i.e., it moves the bi-stability window (W b ) to the left in all cases of miR-451, AMPK, and mTOR and leads to an increase in the size of the window (|W b |) despite disappearance of the bistable switch, which is also illustrated in Fig S2D. This indicates that the probability of switching to the proliferative phase is increased in response to fluctuating glucose levels (0 ≤ G ≤ 1.0) when the inhibition strength of miR-451 by AMPK (β) is increased. For example, the cell is already in the proliferative phase (M > th M , A < th A , R > th R ) when G = 0.2 in the case of higher β (β = 2.5) while it should still be in the migratory phase (M < th M , A > th A , R < th R ) in the base case (β * = 1.0). On the other hand, glioma cells would initiate the motility machinery even under normal glucose conditions when β is decreased. For example, the cell would increase probability of switching to the migratory phase by lowering mTOR levels for a normal glucose level (G = 1.0) when β is lowered (β = 0.07) while it should be in the proliferative phase (M > th M , A < th A , R > th R ) in the base case (β * = 1.0). Fig S2D shows a phase diagram of G − β for switching behaviors including a monostable region and bistable and one-way switches. An increase in the inhibition strength of AMPK by miR-451 (β) induces a transition from a bistable switch to a one-way switch and to a mono-stability. A decrease in β reduces the size of the bi-stability window (|W b |) but an increase in β reduces the size of the one-way switch window. When miR-451 activity is highly enhanced with the increased inhibition signal of AMPK by miR-451 (an increase in β) due to mutual antagonism between miR-451 and AMPK, glucose may also regulate only forward transition from T m -state to T p -state (one way switch) and a monotonic decrease in glucose does not push cells to the migratory phase. Therefore, this mechanism has a similar effect on keeping high miR-451 levels as that of a decrease in α, providing a way of keeping glioma cells in a proliferative phase. In Fig S2E we summarize the phenotypic changes between T m -and T p -groups as β is varied. An increase (or decrease) in β leads to weakening (or strengthening) of AMPK module and up-regulation (or down-regulation) of miR-451 due to mutual antagonism between miR-451 (M ) and AMPK (A), resulting in cell proliferation (or migration). Fig S2F illustrates the effect of the increased and decreased inhibition strength of mTOR activity by AMPK complex (γ = 0.08, 1.0, 9.0) on the regulation of mTOR expression. Recall that cell migration and proliferation heavily depend on glucose levels in the base parameter set, i.e., up-and down-regulation of mTOR expression levels in response to high and low glucose levels (G). However, we found that the dynamics of proliferation and migration are changed when γ is varied. As γ increases from the base value (γ * = 1.0), mTOR activity is significantly reduced, leading to underexpression of mTOR regardless of high and low glucose levels. Therefore, enhanced inhibition of mTOR by AMPK complex (γ > γ * = 1.0 induces complete blocking of mTOR activity (R < th R , ∀G ∈ [0, 1]) followed by cell migration even in the case of high glucose injections (γ = 9.0). On the other hand, cell migration may be inhibited when γ is decreased from the base value. For example, mTOR levels stay above the threshold value (th R = 3.0) even at low glucose levels (G = 0.2 for example) when γ = 0.08. This means glioma cells may stay in the proliferative phase in response to high and low glucose levels when the inhibition of mTOR by AMPK complex is weakened. We also note that bifurcation curves of two other key molecules (miR-451 and AMPK), thus bi-stability windows, are not affected by the changes in γ due to the network structure of the signaling pathway in the signaling pathway (Fig 1). See Fig S2G. The corresponding phenotypic changes between T m -and T p -groups are shown in S2H. As in the case of α, an increase (or decrease) in γ generates higher probability of switching to the migratory (or proliferative) phenotype for a given glucose level by lifting the bifurcation curve down (or up) but keeping the bistable window (W b ).
We recall that the system exhibits a switching behavior, that is, an increase (or decrease) in G induces a transition from T m -status (or T p -status) to T p -status (or T m -status). Steady-state mTOR bifurcation diagrams as a function of the parameter G for various values of λ 1 are presented in the top panels of Fig  S3. It shows the pushing and pulling properties of mTOR (see Figs S3(A-C)). The mTOR level (R s ) decreases with increasing α for the fixed G and λ 1 and the bifurcation curve moves to the right for the fixed λ 1 as seen in Fig S1. While the system exhibits the hysteresis (property of bistable systems) for the base parameter set and the increased λ 1 and various perturbed α's considered (λ 1 = 5.0, α = 1.0, 1.6, 2.0, 2.5). However, when λ 1 is decreased, for larger values of α, the response of the switch is slightly different. For example, the core control system loses the bi-stability property for larger α's (α = 1.6, 2.0, 2.5) when λ 1 = 3.0 (Figs S3A and S3D). Therefore, even for the base value of α (α * = 1.6) in this case, an increase (or decrease) in glucose leads to a smooth transition from T m -status (or T p -status) to T p -status (or T mstatus) along the branch. However, the system also illustrates an extreme case of hysteresis, a one-way switch, when λ 1 is increased (λ 1 =5.0 in Fig S3C) where the upper branches (T p -regime) stretches into the larger domain but actually need to be eliminated due to physiologically irrelevant regimes (G < 0) of glucose (for example, α=1.0 in Fig S3B and α=1.0,1.6 (base) in Fig S3C). In this oneway switching case, glucose withdrawal from the normal level would not induce cell migration from the main tumor mass but the slow supply of glucose from the complete withdrawal condition would put the glioma cell in both T mand T p -status depending on the initial condition of variables (M, A, R) in the core control system.  Table 1.
corresponding profiles of the bifurcation curves of R s as a function of glucose levels (G) for different λ 1 's (λ 1 =3, 4,5). A bistable area is enclosed by two saddle-node bifurcation points (the right and left knees of the curves in the top panels of Fig S3). For any glucose levels in the bi-stability window , the system has two stable solutions and one unstable solution as we saw before for the basic parameter set (Fig 3J). The bottom panels of Fig S3 show phase diagrams of switching behaviors including a monostable region and bistable and one-way switches for different λ 1 values. In all three cases, an increase in the inhibition strength of the miR-451 (α) induces a transition from a one-way switch to a bistable switch, and to a mono-stability. An increase in the autocatalytic production rate of miR-451 (λ 1 ) results in an increase in the size of both oneway and bistable regions in α − G planes. It also moves up the bistable switch region. Therefore, the increased value of λ 1 creates a rich dynamic mixture of oneway switch, bistable switch, mono-stability of the core control system, thus, complex phenotypic switching between T p -and T m -types.
In the top panels of Fig S4 we present steady-state mTOR bifurcation diagrams as a function of the glucose level (G) for various values of β (β=0.8,1.0,2.0). The mTOR level (R s ) decreases with increasing λ 3 for the fixed G and β. An increased autocatalytic production rate of AMPK (λ 3 : 2.5 → 4.0 → 5.5 in Figs S4A-S4C) pushes the bifurcation curve to the right for all fixed β's, increasing probability of establishing the T m -phenotype (R > th R = 3.0). As shown in Fig S2, an increase in β causes an overall shift to the left of those curves, leading to higher chances of acquiring the T p -phenotype. The core control system loses the bi-stability property, leading to a oneway switch when λ 3 is decreased (λ 3 = 2.5). Therefore, this has the same effect on cell fate as a decrease in α above (Figs S3B and S3C), i.e., glucose withdrawal does not guarantee the initiation of glioma cell migration from the main tumor core but the glucose supply from the zero level may generate both T m -and T p -phenotype depending on initial conditions of those key molecules (miR-451, AMPK, mTOR) in the core control system. The bottom panels of Figs S4D-S4F show the phase diagrams of switching behaviors including a monostable region and bistable and one-way switches for different β's (β=0.8, 1.0, 2.0). In all three cases, an increase in λ 3 induces a transition from a one-way switch to a bistable switch, and to mono-stability. An increase in β leads to an increase in the size of both oneway and bistable regions in λ 3 − G planes. It also moves up the bistable switch region. Therefore, the increased value of β creates richer dynamics for larger ranges of β 3 : a mixture of a oneway switch, bistable switch and mono-stability of the core control system, thus, complex phenotypic switching between T p -and T m -types.

Sensitivity Analysis
In the model developed in this paper, there are a number of parameters for which no experimental data are known and/or they may have a significant effect on the results. We consider all parameters (G, S 1 , S 2 , λ 1 , λ 2 , λ 3 , λ 4 , λ 5 , λ 6 , α, β, γ, 1 , 2 ) for sensitivity analysis of the core control model. In order to determine the sensitivity of the main variables (miR-451 level, AMPK activity, and mTOR activity) in the model at certain time to these parameters, we have performed sensitivity analysis. We have chosen a range for each of these parameters and divided each range into 10,000 intervals of uniform length. For each of the fourteen parameters of interest, a partial rank correlation coefficient (PRCC) value is calculated. PRCC values range between -1 and 1 with the sign determining whether an increase in the parameter value will decrease (-) or increase (+) the miR-451 level, AMPK activity, and mTOR activity at a given time. The sensitivity analysis to be described below was carried out using the method from [11] and Matlab files available from the website of Denise Kirschner's Lab: http://malthus.micro.med.umich.edu/lab/usadata/. Fig S5 shows the sensitivity analysis results for the miR-451 level, AMPK activity, and mTOR activity at t = 10 (in (A)), 100 (in (B)). We calculate PRCC values and associated p-values for fourteen perturbed parameters (G, S 1 , S 2 , λ g , λ 1 , λ 2 , λ 3 , λ 4 , λ 5 , λ 6 , α, β, γ, 1 , 2 ) for the model (6)-(8), i.e., we compute PRCC values for the miR-451 level, AMPK activity, and mTOR activity at time t = 10, 100. We conclude that the miR-451 level at t = 100 is positively correlated to the parameters G, λ 1 , λ 2 but is poorly correlated (and thus not sensitive) to S 2 , λ 3 , λ 4 , λ 5 , λ 6 , β, γ, 1 , 2 . Thus, in particular, the miR-   Table 1. indicates PRCC values of the miR-451 level, AMPK activity, and mTOR activity for model parameters (G, S 1 , S 2 , λ g , λ 1 , λ 2 , λ 3 , λ 4 , λ 5 , λ 6 , α, β, γ, 1 , 2 ) at t = 10 (in (A)), 100 (in (B)).
The analysis was carried out using the method of [11] with sample size 10000.
451 activity will increase significantly if the glucose signal to miR-451 (G) and miR-451 autocatalytic production rates (λ 1 , λ 2 ) are increased. On the other hand, the miR-451 activity is negatively correlated to the parameters S 1 , α, i.e., there would be significant decrease in the miR-451 level if the signal to the AMPK complex (S 1 ) and the inhibition strength of miR-451 by the AMPK complex (α) are increased. Due to mutual antagonism between the miR-451 level (M ) and AMPK activity (A), it is expected to see opposite correlations between AMPK activity and those parameters. Indeed, AMPK activity is positively correlated to the parameters S 1 , λ 3 , λ 4 but negatively correlated to glucose signal (G) and inhibition parameter of AMPK (β). AMPK activity is not sensitive to parameters λ 1 , λ 2 , α, 1 . Due to the structure of the model, it might be expected to see sensitivity of mTOR activity to parameters that are sensitive to miR-451 levels. However, mTOR activity is positively correlated to the parameters (S 2 , λ 5 , λ 6 ) but is poorly correlated to G, λ 1 , λ 2 , λ 3 , λ 4 , α, β, 1 , 2 . On the other hand, mTOR activity is negatively correlated to the signaling strength of the AMPK complex (S 1 ) and inhibition strength of mTOR by AMPK complex (γ). This is due to a strong inhibitory effect of LKB-AMPK on mTOR [1,12].

Therapeutic approaches: regulation of migration in glioblastoma
We test the effect of intervention of three pathways on regulation of cell migration and proliferation.
Here, we assume that glucose is delivered to the glioma cell from blood vessels nearby, (ii) drugs can be injected either intravenously or via simple injection directly to the site. For this purpose, we introduce two equations for glucose (G) and drugs (D): where µ G is the consumption rate of glucose by brain tissue, µ D is the natural decay rate of drugs. Here, S G is the external periodic source of glucose from blood vessels on the time interval [t i , t i + h G ], i = 1, . . . , N G , for the time duration h d and a period τ G (= t i+1 − t i , i = 1, . . . , N G − 1; h G < τ G ) between those intervals, N G is the number of glucose injections, I[·] is the indicator function. In a similar fashion, S D is the external periodic source of drugs on the time interval [t j , t j + h D ], j = 1, . . . , N D , for the time duration h D and a period τ D (= t j+1 − t j , j = 1, . . . , N D − 1; h D < τ D ) between those intervals, N D is the number of drug injections. In Figs S6A-S6F, we show the effect of inhibition of miR-451 by AMPK complex by a drug in Eq (11). So, we consider the modified equation where δ is the effectiveness of the inhibitory drug on the pathway. Eqs (10)- (12) are coupled with the core control system (6) (Fig S6A). The same effect of dominant proliferation modes can be obtained when we block the inhibition pathways of mTOR by AMPK complex, i.e., γ → γe −δD (Fig S6C). These results indicate possible strategies of either increasing or decreasing cell migration at the molecular level. For instance, Kim et al. [7,13,14] suggested therapeutic strategies to use miR-451-AMPK core control system in terms of localization of invasive cells immediately after first surgery, followed by secondary surgery for removal of the localized invasive cells [15] on the periphery of original surgical site. An optimal control theory [16] was also applied for the strategy in order to obtain the optimal solution of maintaining 'non'-invasive status given glucose constraints since too much glucose could lead to clinical complications especially for diabetic patients with glioblastoma. AMPK stimulation or inhibition could lead to better clinical outcomes in different types and stages of cancer [12]. For instance, some authors reported that the stimulation of AMPK with pharmacological compounds may induce anti-tumoral effects in experiments [17,18] especially for diabetic type II patients with anti-diabetic drug metformin, toxic for malignant tumor cells [19]. However, activation of AMPK provides better chance of survival in many types of cancer. For instance, AMPK plays a central role in overcoming anoikis [20], programmed cell death under a loss of attachment to the basement membrane. AMPK is also the central coordinator of the metabolic reprogramming known as Warburg effect [21], a hallmark of malignant tumor cells, and AMPK inhibition rather than promotion may represent a potential therapeutic strategy in glioblastoma where the active infiltration of undetected cancer cells leads to final recurrence and poor clinical outcomes [22]. For example, compound C, a reversible and ATP-competitive AMPK inhibitor [23], was used to inhibit AMPK [24]. Glioma cells are exposed to highly heterogenous micro-environment such as fluctuating glucose levels. Experimentally it is not clear how many pathways or how fast action is necessary in order to induce the inhibitory action of miR-451 by AMPK complex. We consider two time delays in inhibition of pathways: (i) inhibition of miR-451 by AMPK complex (τ 1 ) (ii) inhibition of AMPK complex by miR-451 (τ 2 ). In Fig S7 we test the effect of time delays in the inhibition pathways on regulation of cell migration or  proliferation. We use the following modified delay differential equations (DDEs) where the time t was rescaled usingt = 10t for the computational convenience. While miR-451 directly acts on the AMPK complex, inhibition of miR-451 by AMPK complex involves indirect regulations via many other pathways [1], leading to delayed inhibition of AMPK complex by miR-451 (τ 2 ). Intervention of direct action on the AMPK complex by miR-451 will also induce time delays (τ 1 ). Different combinations of those two time delays in inhibition strengths lead to diverse behaviors in terms of cell migration and proliferation. In Fig S7 we investigate the effect of time delays on regulation of cell migration and proliferation in response to an intermediate glucose level (G = 0.5). The model predicts that tumor cells will adapt to the proliferative mode when τ 1 > τ 2 , τ 2 < τ † 2 for a threshold τ † 2 (Fig S7D) while the cells will switch to the migratory phase when τ 1 < τ 2 , τ 1 < τ † 1 for a threshold τ † 1 ( Figure S7A). When the two time delays are small, the system leads to either proliferative or migratory phases depending on the relative strength of τ 1 and τ 2 (Fig S7C). Note that the cell is in the migratory phase for the same initial condition (M (0) = 0.01, A(0) = 0.05, R(0) = 4.0) in the absence of time delays, i.e., when τ 1 = τ 2 = 0 (Fig 3). On the other hand, the system leads to the oscillations between T p and T m when τ 1 , τ 2 are relatively large, for example, when (τ 1 , τ 2 ) = (1.2, 3) in Fig S7B. These results indicate the difficulty of treatment of GBM in the presence of those time delays (τ 1 , τ 2 ) in a microenvironment where glucose may fluctuate, i.e., possible drugs targeting AMPK inhibition may not work for heterogeneous glioma populations and other alternative strategies may need to be developed.
Stochastic fluctuations on signaling networks of miRNAs may play a significant role in regulation of cell fates [25]. In Fig S8 we illustrate the stochastic effect of the miR-451 and AMPK complex on the downstream mTOR activity, and regulation of cell migration in response to fluctuating glucose (G(t) = 0.5 + 0.2 sin(πt/24)) levels. Figs S8A-S8C show one realization for perturbed activities of miR-451 ( Fig S8A) and AMPK complex (Fig S8B) in G − M and G − A planes, respectively, and time course of mTOR levels ( Fig S8C). Gray S -curves in Figs S8A and S8B represent the bifurcation curve in the absence of noise (Fig 3J). While the system shows periodic oscillations in mTOR levels between T p and T m phases in the absence of noise (green solid line), it maintains the T p -phase in the first cycle (interval between arrows in Fig S8C) under the noise. Chances of infiltration of a glioma cell will be positively correlated to the time spent in the T m -phase, increasing the risk of tumor recurrence. Figs S8D and S8E show risk maps without ( Fig S8D) and with (Fig S8E) noise for various combinations of α and β in response to fluctuating glucose and noise in Fig S8(A-C). While the overall pattern of high risk (smaller α and larger β) in Fig S8D is similar to the case with noise in Fig S8E, the intermediate boundary area between T m and T p is larger in Fig S8E. In other words, cell migration is more subject to stochastic noise in miR-451 and AMPK pathways. Finally, Fig S8F shows time courses of mTOR activities with the noise in the inhibition of mTOR activities by AMPK complex for various γ's (γ = γ 0 + σB(t); σ = 0.5, γ 0 = 0.5 (red), 1.0 (green), 2.0 (blue)). We found that the core control system is highly adaptive in terms of selecting the T p and T m phases even under this perturbation on the mTOR inhibition.

Discussion
miRNAs recently emerged as one of the key regulators of cellular process such as the cell cycle [8,26]. Despite recent breakthroughs in our knowledge of these miRNAs, it is not clear yet how many miRNAs Figure S8. Stochastic effect on behaviors of the core system in response to fluctuating glucose levels (G(t) = 0.5 + 0.2 sin(πt/24)). The miR-451 and AMPK levels were perturbed by random Gaussian input (we added σdB(t) on the right-hand side of Eqs (6)- (7) where B(t) is a standard Brownian motion). (A-C) One realization of the stochastic simulations for trajectories in G − M (A) and G − A planes (B), and the corresponding time course of mTOR activity (C) are shown in red. Gray bifurcation curves without noise are shown in (A,B). Green curves represent the trajectories (A,B) and time course (C) in the absence of noise. In response to glucose withdrawal, the system stays in the T p -phase with noise (between arrows) while it switches to the T m phase; shaded light blue box) without noise (green curve). Other parameters are same as in Table 1 are involved in regulation of cancer progression or how these miRNAs exactly control the cellular process in tumor growth overall. Some upregulated miRNAs such as miR-21 in brain tumors are believed to play a pivotal role in pro-oncogenic activities supporting growth, proliferation, migration, and survival of cancer cells while expression of other miRNA is suppressed in gliomas. These miRNAs harbor a clinical significance as therapeutic components in anti-cancer therapy [27][28][29][30][31]. AMPK and mTOR are one of master players in metabolic reprogramming in glioma [32]. Godlewski et al. [1,33] identified a key molecular control system (miR-451, AMPK complex, and mTOR) that provides a critical switch between cell proliferation and migration. They [1,33] found that glucose can regulate glioma cell proliferation and migration via the core control system: (i) low glucose ⇒ down-regulation of miR-451, up-regulation of the AMPK complex, and down-regulation of mTOR ⇒ cell migration; (ii) normal (high) glucose ⇒ upregulation of miR-451, down-regulation of AMPK complex, and up-regulation of mTOR ⇒ proliferation.
We developed a mathematical model of the core control system (miR-451-AMPK-mTOR) based on the experimental observations [1,33] and analyzed the model's behavior in response to high and low glucose levels. The responses of the core control system for various glucose levels are in good agreement with experimental results [1,33]: In particular, (i) the up-and down-regulation of miR-451 and mTOR in response to high and low glucose levels, respectively, and (ii) the down-and up-regulation of the AMPK complex in response to normal and low glucose. This allowed us to define the migratory and proliferative phase based on the status of these molecules and the hysteresis system generates a bistable window (W b in Fig 3J) as predicted in the previous smaller miR-451-AMPK model [6,7]. The model system predicts oneway-, bistable-, and mono-stability switches under the perturbations of the key inhibition parameters α (Fig S1), β (Fig S2), γ (Fig S2)). Phenotypic changes under the perturbation of these parameters (α, β, γ) enabled us to predict the T m /T p -status of glioma cells in the glucose fluctuating microenvironment when specific components of the signaling system are perturbed, for example by miR-451-suppressing drugs (Figs S1,S2). The changes in the size and location of the bistable and one-way regions in G − α plane under perturbations of the inhibition (α) and enhancing (λ 1 ) parameters of miR-451 also generate a rich dynamics where the cell fate of migration and proliferation critically depends on the sign of dG/dt, a decrease in glucose (steady glucose withdrawal by the consumption) or increase in glucose (supply from blood vessels (BVs) or glucose intravenous infusion), and initial conditions (Fig S3). We also investigated the system responses under the perturbation of the inhibition (β) and enhancing (λ 3 ) parameter of the AMPK complex ( Fig S4). Onset of T p and T m phenotypes in response to drugs targeting the inhibition pathways of the AMPK complex in the presence of fluctuating glucose are sensitivity to duration and strength of the drug while only T p -phenotype is selected for drugs targeting the inhibition pathways of miR-451 or mTOR (Fig S6). We also illustrated the difficulty of treating GBM in the presence of time delays (τ 1 , τ 2 ) in the inhibition signaling network due to diverse responses of GBM cells (proliferation, migration, or periodic changes in T p and T m phases; Fig S7). In order to eradicate all heterogeneous glioma populations, alternative strategies may need to be developed. For instance, pre-conditioning chemotaxis-driven localization of cells [13,14] after conventional surgery might be a way of homogenization of the tumor population in order to increase the efficacy of the drugs or clinical outcomes [34]. In this case, the optimized solutions need to be obtained in the controlling the core control system [16].
The growth pattern of a GBM can be determined based on perivascular accumulation of cancer cells, perineuronal satellitosis, subpial growth, and intrafascicular growth in the corpus callosum [22]. While miR-451 levels are up-regulated near blood vessels due to high levels of glucose supply [1,33], glucose levels and stiffness in brain tissue are irregular due to the heterogeneous microenvironment. Many of these microenvironmental factors in addition to intrinsic noise in vast signaling networks including the miR-45-AMPK-mTOR system may introduce biochemical perturbations within the core control system and the cell may select the proliferative phase instead of the migratory phase even under glucose conditions favoring migration (Fig S8). Cancer cells within tumors exist in distinct phenotypic states and stochasticity in individual behaviors may promote phenotypic equilibrium [35] and affect transitions between