Strategies in regulating glioblastoma signaling pathways and anti-invasion therapy

Glioblastoma multiforme is one of the most invasive type of glial tumors, which rapidly grows and commonly spreads into nearby brain tissue. It is a devastating brain cancer that often results in death within approximately 12 to 15 months after diagnosis. In this work, optimal control theory was applied to regulate intracellular signaling pathways of miR-451–AMPK–mTOR–cell cycle dynamics via glucose and drug intravenous administration infusions. Glucose level is controlled to activate miR-451 in the up-stream pathway of the model. A potential drug blocking the inhibitory pathway of mTOR by AMPK complex is incorporated to explore regulation of the down-stream pathway to the cell cycle. Both miR-451 and mTOR levels are up-regulated inducing cell proliferation and reducing invasion in the neighboring tissues. Concomitant and alternating glucose and drug infusions are explored under various circumstances to predict best clinical outcomes with least administration costs.


Introduction
Glioblastoma multiforme (GBM) is the most common and the most aggressive type of brain cancer. The median length of survival time is approximately 12 to 15 months following diagnosis. GBM is characterized by anaplasia, nuclear atypia, cellular pleomorphism, mitotic activity, and more importantly, alternating phases of rapid proliferation and aggressive invasion into the surrounding brain tissue. This leads to an inevitably critical recurrence even after the surgical resection of the main tumor mass [1,2]. The mainstay of treatment for GBM is surgery, followed by radiotherapy and chemotherapy. Despite advances in these approaches, glioma cells can still invade the neighboring tissues beyond detection leading to tumor recurrence. High probability of main treatment failure also encourages researchers to investigate the use of innovative treatments when the first line of therapy has failed, in order to improve clinical outcomes [3].
In the tumor microenvironment (TME), glioma cells encounter many challenges including hypoxia, acidity, and limited nutrient availability. To maintain rapid growth, tumor cells need to adapt to these biochemical changes and modify their metabolic activity by increasing glycolysis even in the presence of oxygen. This process is called the Warburg effect which requires consuming considerable amounts of glucose [4]. The tricarboxylic acid cycle, or Krebs cycle, plays an important role in the breakdown of organic fuel molecules and the survival in non-a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 hypoxic normal differentiated cells. These molecules include glucose, fatty acids, and some amino acids. While differentiated cells favor this type of metabolism, which is very efficient in terms of ATP production, tumor cells adopt the inefficient aerobic glycolysis producing relatively large amounts of waste product in the form of lactic acid [5]. This may provide cancer cells the advantage of not having to depend on oxygen as an energy source especially in a hostile tumor microenvironment, thus leading to longer survival [6]. Inhibition of glycolysis may also prevent drug resistance thus a better understanding of this metabolic pathway may lead to better treatment options and clinical outcomes [7]. Developing strategies of metabolic adaptation, angiogenesis, and migration is critical for cancer cells in order to survive metabolic stress and ensure enough nutrient supply as tumor mass accumulates where glucose supply may fluctuate due to heterogeneous biochemical and biophysical conditions [8]. Therefore, adequate cellular responses to glucose withdrawal are critical for cancer cell survival. Cancer cells then activate the 5 0 -adenosine monophosphate activated protein kinase (AMPK) pathway under metabolic stress. It is the master cellular sensor of energy availability which enhances glucose uptake and conserve energy, thus avoiding cell death [9]. MicroRNAs, also abbreviated as miRNA, are approximately 22 nucleotide single-stranded non-coding ribonucleic acids (RNAs) that are known to regulate gene expression [10]. Dysregulation of microRNA expression has been linked to oncogenic and tumor suppressor activities in several types of cancer, including GBM [11,12], where altered miRNA expression contributes to tumorigenesis [13,14].
Godlewski et al. [8] identified an interesting mechanism of glioma cell migration and proliferation wherein a particular microRNA, miR-451, and its counterpart, AMPK complex (CAB39/ LKB1/STRAD/AMPK), determine whether the cell favors growth at the expense of invasion or conversely. Moreover, they also identified a potential feedback loop between LKB1 and miR-451 allowing for a sustained and robust response to glucose withdrawal [15]. It was found out that (i) under high (normal) glucose conditions, up-regulation of miR-451 leads to the down-regulation of AMPK complex, which then leads to elevated proliferation and decreased migration of glioma cells and (ii) glucose withdrawal induces down-regulation of miR-451 and up-regulation of AMPK, which promotes cell migration with reduced proliferation.
The mathematical models developed by Kim et al. [14,16,17] describe the effects of the miR-451-AMPK core control system on cell proliferation and migration in glioblastoma. It explains the response of miR-451 to high and low glucose levels as well as the mutual antagonism between miR-451 and AMPK complex concentrations. Kim et al. [18] then extended the model to include the dynamics of mammalian target of rapamycin (mTOR), which is a protein kinase that links with other proteins as well as to the cell cycle dynamics, to form the miR-451-AMPK-mTOR core control system. This mutual antagonism between miR-451 and AMPK, which was predicted in these mathematical models [14,[16][17][18], was confirmed by recent experiments. For example, Ansari et al. [19] recently found that (i) the miR-451 transcription in GBM cells is induced by unrestricted activity of its transcription factor OCT1 (official gene symbol POU2F1) in the presence of abundant glucose, resulting in AMPK inhibition through direct targeting of CAB39 in the LKB1 complex; and (ii) the miR-451 level is inhibited through the phosphorylation and inactivation of OCT1 at S335 by AMPK in response to glucose depletion-induced metabolic stress, leading to a reciprocal negative feedback loop between miR-451 and AMPK. In this case, suppression of miR-451 in turn leads to sustained AMPK activities and a robust response to glucose withdrawal in GBM cells. The multiscale mathematical models [14,[16][17][18] also predict the growth-invasion cycling patterns of glioma cells in response to fluctuating glucose uptake in the tumor microenvironment. The core control system predicts bistability and hysteresis bifurcation when delayed down-regulation of miR-451 activities along certain molecular pathways would induce glioma cells to stay longer in the proliferative phase despite relatively low glucose concentrations, making this mechanism a therapeutic target.
The cell cycle represents an integrated series of events that regulates complex processes including cell proliferation, cell division and DNA replication, regulated by a complex hierarchy of genetic and metabolic networks which involves several transition states of varied lengths and checkpoints [20]. The stages of the cell cycle are as follows: (i) synthesis phase (S), a period where DNA replication occurs; (ii) gap phase 2 (G2), during which proteins required for mitosis are produced; (iii) mitosis phase (M), a period where chromatin condensation, nuclear envelope breakdown (NEBD), chromatid separation, and cytokinesis happens; (iv) gap phase 1 (G1), in which genes necessary for DNA replication are activated and the protein agents of S phase progression are accumulated; and (v) resting phase (G0), a state in which cells can exit the cell cycle and enter a phase of quiescence or relative inactivity [21]. The progression of mammalian cell cycle is tightly regulated by coordinated activation of cyclin-dependent kinases (CDKs) family [22]. The CDKs are positively regulated by cyclins and negatively by CDK inhibitors (CDKIs) such as the proteins p15, p16, p21 and p27. In cancer cells, cyclins are over-expressed while CDKIs are under-expressed which results in the dysregulation of the cell cycle, and promoting uncontrolled cell growth [20]. Tyson and Novak [23,24] identified that the transition between two stable steady states, G1 and S-G2-M cell-cycle phases, are described using the kinetic relations of the model that is controlled by changes in cell mass.
A standard GBM treatment is surgery followed by chemotherapy and radiotherapy. However, even under best circumstances, the mean survival of this disease is about a year. Poor outcomes of standard care treatments are due to the topographically diffuse nature of the disease [25]. By the time of diagnosis, typical GBM cells may have widely spread throughout the brain tissue [26][27][28], increasing the potential of recurrence. Thus, annihilation of distant tumor satellites is implausible despite surgically removing all the essential tumor seen on enhanced MRI scan [29]. Knowing the exact margins of a tumor mass in real patients is indeed a daunting task. In this study, it is assumed that a major tumour mass has been surgically removed and that the infiltrative tumor cells are near the surgical site. The objective is to prevent the glioma cells from further diffusing into the surrounding brain tissue. Localization approach will be utilized, that is, glioma cell invasion will be blocked keeping them in a proliferative phase while also attempting to limit excessive growth before a second surgery [17]. Our analytical tool is primarily based on the framework of optimal control theory, which has been successfully used to make informed decisions involving biological models such as optimal treatment strategies in human immunodeficiency virus (HIV) models [30][31][32][33], tuberculosis [34][35][36], and cardiopulmonary resuscitation (CPR) techniques [37,38]. Optimal control theory is also applied to the miR-451-AMPK core control system to determine the intravenous glucose and/or drug infusion protocols with least possible cost under various circumstances in de los Reyes, et al. [39].
Recently, Kim et al. [40] developed an intracellular signaling pathway model that extends the miR-451-AMPK-mTOR core control system including the cell cycle dynamics. In this work, a potential control problem is formulated in order to maintain high levels of miR-451 and mTOR (low levels of AMPK) inducing cell proliferation prohibiting cell motility and invasion to the neighboring tissues. With glucose levels as a key regulator of miR-451 activity which also activates mTOR in the downstream signaling pathway, glucose intravenous infusion is considered to up-regulate miR-451 and mTOR concentrations above certain threshold values. In addition, a drug suppressing the inhibitory effect of mTOR by the AMPK complex is incorporated. This drug can be administered concomitantly or alternately with glucose as a secondary intravenous infusion. The controls are then given by dose rates of glucose and drug intravenous administrations regulating upstream and downstream signaling pathway, respectively. Solution of the optimal control problem aims to determine infusion protocols with minimal glucose and drug amount, and least administration costs. Hence, glucose and drug levels are regulated to prevent rapid tumour growth, hyperglycemia, and further drug complications.
The results propose plausible glucose and drug intravenous infusion controls which indicate the time of administration, frequency, number of administrations, and dosages of glucose and drug per infusion.
In the current study, we will first present the miR-451-AMPK-mTOR-cell cycle intracellular signalling pathway developed by Kim et al. [40]. A drug module is incorporated to up-regulate mTOR activities inducing cell proliferation. Then an optimal control problem is formulated with the goal of activating miR-451 and mTOR levels through glucose and drug intravenous infusions. Two different infusion protocols will be explored, namely, concomitant and alternating glucose and drug intravenous administrations. Optimal solutions for two strategies are presented and results on frequency, dosage per infusion, total glucose and drug amount, and relative cost incurred in the administrations are compared. The conclusion section discusses and summarizes the optimal control results, and provides outlook for future research directions.

Intracellular signaling dynamics model
In this section, we present the basic components of intracellular signalling pathway of tumor cell containing the core control miR-451-AMPK-mTOR and cell cycle pathway developed in Kim et al. [40], incorporating a drug which blocks the inhibitory pathway of mTOR by AMPK complex. The core control model identifies a key mechanism which determines the molecular switches between the proliferative and migratory phases in response to fluctuating glucose and drug levels. The simplified signaling pathways consists of five key determinants of the intracellular structure, namely, glucose level G, miR-451 level M, AMPK complex activity A, mTOR concentration R, and drug level D. The intracellular cell cycle dynamics are developed by Tyson and Novak [23,24] including only the essential interactions for its regulation and control. The model captures the kinetics of chemical processes within the cell such as production, destruction, and different molecule interactions. The transition between two main steady states, G1 and S-G2-M phases, of the cell cycle are described using the kinetic relations of the model that is controlled by changes in cell mass. . Also included are the effects of the changes in oxygen dynamics at the macroscopic level through the activation and inactivation of HIF-1 α. This results in changes in cell cycle length. In addition, the cell cycle inhibitory effect of p21 or p27 genes of HIF-1 α is incorporated. Kim et al. [40] proposed to link both the miR-451-AMPK-mTOR control system and the cell cycle dynamics to provide a mechanism driving the cell cycle to undergo the quiescent stage G0-phase depending on the concentration level of mTOR. Fig 1A depicts the detailed schematic diagram of the intracellular signalling networks including miR-451, AMPK complex, mTOR, and key players in the cell cycle module (CycB, Cdh1, p55cdcT, p55cdcA, Plk1). Kinetic interpretation of arrows and hammerheads represent induction and inhibition in the signaling network, respectively. The dimensionless diagram of the core control miR-451, AMPK complex, and mTOR linking to the cell cycle is depicted in Fig 1B. It should be noted that u G and u D are the sources of glucose and drug with decay rates μ G and μ D , respectively, which can be controlled exogenously. S 1 and S 2 are signaling sources to AMPK complex and mTOR, respectively, while α, β and γ are inhibition strengths, and ϕ denotes decay.
It has been shown that high (normal) glucose concentration yields over expression of miR-451 levels (down-regulation of AMPK complex and up-regulation of mTOR) leading to elevated cell proliferation and reduced migration, while low glucose levels leads to down-regulation of miR-451 activities (up-regulation of AMPK complex and down-regulation of mTOR) reducing cell proliferation and inducing migration into the surrounding brain tissue [8,15,18,40]. The effect of various glucose levels in the regulation of the core control is depicted in Fig 2. Kim et al. [14] developed a core control system of glioma cell migration and proliferation by using a regulatory network of key molecules (miR-451 (M), AMPK (A)) as follows: As observed in experiments [8,15,19], the regulatory system in Eq (1) includes a mutually antagonistic loop between miR-451 and AMPK complex in response to high and low glucose levels (G). This genetic toggle switch induces a monostable and bistable system, characterizing proliferation and critical cell infiltration of GBM cells in brain [14]. In the follow-up studies [14,[16][17][18] including mTOR (R) and cell cycle modules, this mutual antagonism and bistable system played a critical role in developing anti-invasion strategies. A mutually antagonistic feedback loop has been well-studied for its bistable properties by use of mathematical models ([42-45] and other references in [45]). For instance, in a study on a genetic toggle switch in Escherichia coli [42], a mutually inhibitory loop between repressor 1 and repressor 2 was  In this study, we consider a drug D suppressing the inhibitory effect of mTOR by AMPK complex where the inhibition strength is given by z(D) = e −D . When D is large, γe −D is small which makes dR/dt to be large. Thus, the presence of drug up-regulates mTOR activity that could eventually lead to cell proliferation. Fig 3 illustrates the mTOR bifurcation curve. Observe that increasing drug concentration shifts the hysteresis curve upwards keeping the same bistability window. Hence, with higher drug levels for the same glucose concentrations, mTOR is activated prompting elevated cell growth.
In the miR-451-AMPK-mTOR core control system developed by Kim and colleagues [14,[16][17][18]40], the bistability regime of main variables emerges in response to glucose levels. As it was shown in [40], the existence and size of the bistability window (|W b |) depend on other essential parameters and may disappear under perturbations of parameters. As it will be shown later (Figs 4 and 5), some of key parameters are sensitive in creation or destroying the bistability while other parameters are not. After achieving the equilibrium, continuation of the curve is computed by varying the glucose concentration G and bifurcation points are detected labeled LP and CP for limit and cusp points, respectively. Fig 6A illustrates the hysteresis diagram (G, R)−curves for different S 1 values. In order to obtain the cusp point (CP), fold continuation is computed starting at a limit point where two parameters G and S 1 are activated resulting to codim 2 bifurcations. Both parameters (G, S 1 ) are varied along the curve where each point is a limit point for the equilibrium curve at the corresponding value of S 1 . This is depicted in Fig 6B. The cusp point gives the threshold values th M , th A , and th R . A Matlab software MatCont was used for numerical continuation and bifurcation study of continuous and discrete dynamical systems [46].
The governing model equations for the dimensionless intracellular signaling dynamics are then described by the following ordinary differential equations d½Plk1� dt ¼ k 9 ½mass s �½CycB�ð1 À ½Plk1�Þ À k 10 ½Plk1�; Strategies in regulating glioblastoma signaling pathways and anti-invasion therapy The cell cycle dynamics and the regulatory core control system are linked by a variable called pseudo-mass ([mass s ]) given by The oxygen dynamics through HIF-1 α is described as and the growth rate μ is expressed as wherem is the probability density function with uniform distribution between −1 and 1. This growth rate formulation introduces cell cycle heterogeneity of length between 20 and 30 hrs to account for the natural variability between cell growth rates and to have a non-synchronous population [41]. The model parameters in system (2) are listed in Tables 1 and 2. In the current model formulation, mTOR (R) is activated/inactivated in the same way as miR-451 (M). In addition, R links the core control system and the cell cycle dynamics via the pseudo-mass ([mass s ]) which influences the cell mass [mass] and the intracellular proteins (refer to Eq (3)). Therefore, when glucose supply is high, R is up-regulated (proliferative phase; up-regulated M, down-regulated A) and [mass s ] � [mass] yielding a typical cell cycle (see Fig  7). On the other hand, when glucose supply is low, R is down-regulated (migratory phase; down-regulated M; up-regulated A) and [mass s ] exceeds [mass] influencing the cells to enter into resting phase G0 [40].
Let us assume that glucose (G) and drug (D) can be regulated through intravenous infusions and can be periodically administered. For illustrative purposes, consider 3h infusions every 12h with maximum dosage of 1 unit, and refer this as regular infusion. We then have The intracellular dynamics under regular infusion can be seen in Fig 8. Since glucose and drug levels oscillate, it follows that M and R also periodically fluctuates around the threshold (th M � 1.87, th R � 2.76) as can be seen in Fig 8A. Note that when M and R crosses the threshold value from above, respectively, peak in pseudo-mass is generated. It can thus be inferred that these peaks indicate cell migration since M and R levels are below their respective threshold values. Strategies in regulating glioblastoma signaling pathways and anti-invasion therapy As shown in Fig 8B, a trajectory of core control concentrations in mTOR-miR-451-AMPK space switches between proliferation and migration region. A closer look at the intracellular cell cycle proteins, mass, and mass s dynamics are illustrated in Fig 8C. Note that regular infusion significantly perturb the cell cycle dynamics stimulating several cell divisions (see mass profiles) and cell migration (see mass s ). Indeed, fluctuating glucose levels in the microenvironment leads to the dichotomy of grow and go dynamics of glioblastoma cells yielding bigger tumor mass as reported in [14], even in the presence of (fluctuating) drug concentrations to keep the cells in proliferation phase.

Optimal control problem
In the current investigation, we aim to regulate the amount of glucose and drug infusions to up-regulate miR-451 and mTOR above its threshold values inducing cell proliferation Strategies in regulating glioblastoma signaling pathways and anti-invasion therapy avoiding migration to neighboring tissues. The modeling approach utilizes optimal control theory to identify infusion administration protocols. Let u G (t) and u D (t) be the controls of the system representing dose rates of glucose and drug intravenous administrations, respectively. Two different administration protocols are examined, namely, (1) concomitant, and (2) alternating glucose and drug infusions, with the following objective functionals respectively. Here, M(t) and R(t) denote the level of miR-451 and mTOR concentrations, respectively. Parameters B 1 and B 2 are weight factors measuring the relative cost based on maximizing miR-451 M(t) and mTOR R(t), and administering glucose and drug intravenous infusions over a specified time interval, respectively. The control costs are modeled by the linear combination of quadratic terms u 2 G ðtÞ and u 2 D ðtÞ. Our objective is to find optimal infusion regimen for glucose and drug administrations, denoted by u � G ðtÞ and u � D ðtÞ, such that the objective functionals are satisfied, that is, where The bounds for the controls represent the limits on dose rates for glucose and drug administrations. Assuming that miR-451, AMPK, and mTOR responses are regulated by glucose levels and influenced by drug levels, our control strategies deal with finding optimal control regimens for both glucose and drug intravenous infusions. We also note that the existence of optimal controls is guaranteed by standard results in control theory [51]. In this maximization problem, the necessary convexity of the integrand in the objective functional holds. Therefore, we can proceed with applying Pontryagin's Maximum Principle [52]. An iterative method is used for solving the optimality system which is a two-point boundary value problem having initial conditions for the state variables and terminal conditions for the adjoints. Numerical simulations are obtained using a fourth-order iterative Runge-Kutta method. Given the initial conditions and guess for the controls, state equations are solved using the forward scheme while the corresponding adjoint equations are solved using the backward scheme with the transversality conditions. The controls are updated by using a convex combination of the previous controls and the value from the characterizations. This is commonly referred as the Forward-Backward Sweep Method (FBSM) which is shown to be convergent [53]. Further details on the optimal control under study can be found in the S1 Appendix.  (2)) was performed in order to determine which parameters have the most/least influence on the reference output (main variables). All the core control parameters were considered to assess their corresponding influence on the miR-451, AMPK, and mTOR activity. It was inferred that miR-451 activity will be enhanced by increasing glucose signal (G) and autocatalytic production rates (ℓ 1 , ℓ 2 ) of miR-451. These parameters were negatively correlated with AMPK activity due to the mutual antagonistic mechanism between miR-451 and AMPK complex. In addition, AMPK activity will be up-regulated by an increase in signaling source S 1 and autocatalytic rates (ℓ 3 , ℓ 4 ) of the AMPK complex in the model. It has been also noted that increasing S 1 and the inhibition strength of mTOR by AMPK (γ) down-regulates mTOR level. In a similar fashion, S 2 is strongly positively correlated with mTOR levels but little correlated with G, ℓ 1 , ℓ 2 , ℓ 3 , ℓ 4 , α, β, � 1 , � 2 at time 100.

Sensitivity analysis on bistability of the miR-451-AMPK-mTOR system
In this work sensitivity analysis is carried out to determine which model parameters have consequential effect in achieving or inhibiting bistability in the (G, R)−curve. As in [40], a method from [54] is adapted where a range for each parameter is selected and divided into N intervals of uniform length. For each parameter of interest, a partial rank The PRCC values range between -1 and 1 with the sign determining whether a change in the parameter value will promote (+) or suppress (−) bistability. The following algorithm is performed: 1. Assign a probability distribution D½m j;min ; m j;max � to each parameter μ j and let N be the number of samples to be selected. Divide the interval [μ j,min , μ j,max ] into N equiprobable subintervals, and draw an independent sample from each subinterval.
2. Assemble the LHS matrix L, wherein each row of L represents a unique combination of parameters sampled without replacement.
3. For each row of the LHS matrix L, solve for mTOR R in terms of glucose G and check for bistability window W b . If bistability exists assign 1 to output variable W b in matrix Y, else assign −1.
4. Rank-transform the matrices L and Y to obtain L R and Y R . By rank-transform, we mean to replace the value by its rank when the data are sorted from lowest to highest, e.g. the smallest value is assigned a rank 1. Tied values are assigned an average rank. Strategies in regulating glioblastoma signaling pathways and anti-invasion therapy 5. Fix a parameter μ j , which is encoded in the j th column in the matrix L R . Form the following linear regression models using the data matrices L R and Y R for μ j and y, respectively: Compute ðm j Àm j Þ and ðy ÀŷÞ, the residuals in the input parameter and the output after removing the linear effects of the other input parameters.
6. Obtain the PRCC of μ j using PRCC m j ;y ¼ Covðm j Àm j ; y ÀŷÞ 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 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Varðm j Àm j ÞVarðy ÀŷÞ 7. Repeat Steps 5 and 6 for the remaining parameters. It illustrates that a change in the miR-451 autocatalytic production rate (ℓ 1 ) or inhibition strength of AMPK complex by miR-451 (β) enhances bistability. On the contrary, the Hill-type coefficient ℓ 4 in the regulation of AMPK is responsible for losing bistability. It should be noted that due to the model's structure, ℓ 1 up-regulates miR-451 which in turn increases mTOR and suppress AMPK activity. The parameter ℓ 4 promotes AMPK complex and down-regulates miR-451 and mTOR. However, other parameters ℓ 2 , α, ℓ 3 , S 1 , ℓ 5 , ℓ 6 , γ, S 2 are little sensitive in emergence or destruction of the bistability.
The bistability of miR-451-AMPK-mTOR core control system depends on the geometric structure of its nullclines. In particular, bistability arises when M-and A-nullclines (i.e., dM dt ¼ 0 and dA dt ¼ 0) intersect at three distinct points, producing one unstable and two stable steady states. The nullclines intersect three times due to their sigmoidicity influenced by catalytic rates ℓ 1 and ℓ 4 . These rates must be proportionate for a given glucose G level. Otherwise, the nullclines will intersect only once. This bistability condition has been shown for a genetic toggle switch in E. coli [42]. Fig 5A-5C depict the M-and A-nullclines under different G levels with various ℓ 1 and ℓ 4 values. The bifurcation curves for ℓ 1 , ℓ 4 and region of bistability for specific G levels are depicted in Fig 5D. Under given circumstances, ℓ 1 and ℓ 4 showed significant sensitivities in the bistability of the system. On the contrary, both ℓ 5 and ℓ 6 affect the R-nullcline only. It is illustrated in Fig 5E that for several ℓ 5 and ℓ 6 values, R-and A-nullclines intersect at only one point, leading to a single steady state. In effect, M-and A-nullclines will also intersect once, producing monostability. Hence, ℓ 5 and ℓ 6 show insignificant consequence in achieving bistability (see Fig 4).
In the following numerical experiments, it is assumed that initially glioma cells are in growth phase (a probable occurrence of post primary tumour surgery) with M > th M , A < th A , and R > th R . It is also considered that glucose and drug can be administered exogenously as intravenous infusions. Further, the weight parameters B 1 = 1 and B 2 = 1 are used as default values unless specified.

Concomitant infusion control
In this strategy, glucose and drug infusions are administered simultaneously, in particular, both controls u G (t) and u D (t) are infused at the same. Thus, this is referred to as concomitant control. Here, it is assumed that the drug in consideration which blocks the inhibitory pathway of mTOR by AMPK complex had negligible side effects and had inconsequential chemical reactions with glucose. In order to determine an efficient strategy of concomitant infusion protocol, glucose and drug are administered concurrently at initial time t = 0 for 3 hours using the numerical scheme FBSM. Infusion spontaneously increases glucose and drug concentrations as depicted in Fig 9(A) and 9(B). Consequently, miR-451 and mTOR levels are up-regulated while AMPK complex is down-regulated as shown in Fig 9C. A closer look at the control curves shows that both u G (t) and u D (t) decrease from 0 < t i < 3 suggesting that both glucose and drug dose rates should be decreased from time t i . This leads to the decrease of glucose and drug concentrations due to consumption and decay. Accordingly, miR-451 and mTOR levels decrease while AMPK complex increases. Before miR-451 crosses the threshold value th M , FBSM is again applied for the next 3 hours. This suggests a time for the next concomitant administrations in which miR-451 profiles are monitored subsequently. The procedures in tracking miR-451 profiles and applying FBSM for 3 hours are repeated over a specified time duration. Thus, the number of concomitant administrations are determined. It is important to note that keeping miR-451 level above its threshold value consequently confine AMPK and mTOR levels, below and above its corresponding threshold values, respectively (see Fig 9C).

Alternating infusion control
Suppose that concurrent glucose and drug administrations is not plausible due to unwanted chemical reactions. We propose another control strategy that administers glucose u G (t) and drug u D (t) infusions alternately. At initial time t = 0, glucose infusion is obtained by solving the optimal control problem using the numerical scheme FBSM for 3 hours. Next, drug infusion is attained for the next 3 hours in a similar manner. Hence, the controls u G (t) and u D (t) are applied one after the other. Subsequently, miR-451 levels are then tracked and before it crosses the threshold value, glucose infusion u G (t) and drug infusion u D (t) are again administered alternately in a similar fashion above, over the specified duration of administration. Thus, the time for the next alternating infusion is determined by tracking miR-451 levels. This strategy is referred simply as alternating control. This infusion protocol is shown in Fig 11(A) and 11(B). Note that glucose infusion increases miR-451 levels and drug infusion activates mTOR activities keeping AMPK down-regulated as depicted in Fig 11C. The intracellular dynamics under alternating control is illustrated in Fig 12. As shown in Fig 12A,

Comparison between control strategies
Both concomitant and alternating control strategies are able to sustain elevated miR-451 and mTOR levels above their threshold values and AMPK levels below its threshold value. It was also shown that mTOR-miR-451-AMPK trajectory is restrained in the proliferation region prohibiting cell migration. Further, number of cell divisions is reduced with slightly higher mass concentration before division compared to regular infusions but lower compared to constant (intermediate) high glucose concentrations. In this section, we compare the cost efficiency of the proposed strategies in terms of frequency of administration, dose per infusion, total glucose and drug amount, relative cost per infusion, and total cost incurred in the administrations. For the following results, the time for simulation duration is considered to be 168 hours (7days).
Recall that parameters B 1 and B 2 are the weight factors associated in our objective functional which represent the measure of costs involved in the administration of glucose u G (t) Strategies in regulating glioblastoma signaling pathways and anti-invasion therapy and drug u D (t) infusions, respectively, which also includes dosage, type, brand, medical fee for administration, etc. Fig 13A shows that as the cost of glucose administration becomes expensive (increasing B 1 values) with fixed drug administration cost (B 2 = 1.0), frequency of concomitant and alternating infusion increases. However, it should be observed that as frequency of administration increases, the optimal glucose dose per infusion decreases with drug dose per infusion increases, see Fig 13B. Increasing drug dosage compensates the decreasing glucose dosage in order to keep mTOR activities up-regulated leading to cell proliferation. On the contrary, if drug administration cost increases (increasing B 2 values) with fixed glucose administration cost (B 1 = 1.0), frequency of administration remains almost constant. This is depicted in Fig 13C. The glucose dose per infusion is almost the same but the drug dose per infusion decreases with increasing B 2 as illustrated in Fig 13D. Further, note that frequency of concomitant infusion control is always lower than that of alternating infusion control. In addition, drug (glucose) dose per infusion of concomitant control is always more (slightly less) than that of alternating control.
For increasing glucose administration cost (increasing B 1 ) with fixed drug administration cost (B 2 = 1.0), total amount of glucose used for both control strategies generally decrease (except for high B 1 > 3.0 values), as depicted in Fig 14A, while total amount of drug increases, see Fig 14B. On the other hand, when drug administration becomes expensive (increasing B 2 Strategies in regulating glioblastoma signaling pathways and anti-invasion therapy values) with fixed glucose administration cost (B 1 = 1.0), the total amount of glucose dosage is almost constant (just slightly increasing for concomitant control) as shown in Fig 14C, but the total drug amount is decreasing as can be seen in Fig 14D. As illustrated in Fig 14, concomitant control use less total amount of glucose than that of alternating control. Contrarily, concomitant administration consumes more total drug amount as compared to that of alternating control (except possibly for high glucose administration cost, B 1 > 3.0, see Fig 14B).
Figs 15 and 16 reflect the relative cost per infusion and total administration cost of glucose and drug infusions incurred under concomitant and alternating controls. With fixed drug administration cost (B 2 = 1.0), relative glucose cost per infusion and total administration cost increases as B 1 increases. Note that alternating control cost for glucose is always higher than concomitant infusions, see Fig 15(A) and 15(B). The relative and total administration cost for concomitant infusion slightly increase as compared to alternating control as B 1 increases. Again, alternating control incur more cost for higher glucose administration cost as illustrated in Fig 15(C) and 15(D). On the other hand, when B 1 = 1.0 and drug administration cost (B 2 ) increases, the relative glucose cost per infusion and total administration cost is almost constant. This can be seen in Fig 16(A) and 16(B). Again, relative total cost for alternating control is higher than that of concomitant control. Further, as B 2 increases, relative drug cost per infusion and total administration cost decreases (Fig 16(C) and 16(D)), since total drug dosage also decreases as shown in Fig 14D. In this case, it can be seen that drug administration cost for alternating control is generally lower than that of concomitant control.
Observe that in Fig 17, both concomitant (blue) and alternating (orange) control trajectories in glucose-mTOR-drug space are restricted in a smaller region avoiding aggressive invasion, rapid proliferation, and unwanted drug complications. Both strategies achieved the goal of keeping mTOR (miR-451) up-regulated inducing cell proliferation and thus avoiding aggressive cell migration. It is important to note that under these proposed optimal control infusions, glucose and drug levels are regulated to prevent excessive cell division and tumor growth. As a consequence, these strategies suggest safer infusion administration preventing hyperglycemia for diabetic patients and risk of drug complications. Table 3 provides the average frequency, dosage and relative cost of concomitant and alternating control strategies where B 1 = B 2 = 1.0 and simulation time is 168h (7d).

Conclusion
The periodic switching behavior of glioblastoma cells between proliferation and invasion phases is highly influenced by fluctuating glucose levels [14,17]. In response to high glucose supply, miR-451 and mTOR are up-regulated and AMPK complex is down-regulated inducing cell growth. On the contrary, low glucose level up-regulates AMPK complex, down-regulating miR-451 and mTOR, promoting cell migration [15]. The mutual antagonistic mechanism Strategies in regulating glioblastoma signaling pathways and anti-invasion therapy between miR-451 (mTOR) and AMPK complex and the cell's strategic metabolic adaptation support the survival of cancer cells even in a nutrient-deprived microenvironment [14,55].
In addition to rapid proliferation of glioblastoma cells, aggressive invasion to the surrounding tissue is a major cause of treatment failure. Despite advances in medical imaging technology such as MRI and PET, glioblastoma cells can spread beyond detection leading to tumor recurrence within 2 to 3 cm of the resection cavity even after surgical removal of a malignant glioblastoma [29]. These glioma cells are capable to deform cell membrane and nucleus for cell infiltration through a narrow gap between normal glial cells in brain tissue by upregulation of myosin II along with actin bundles [26,56]. While exact migratory patterns are still poorly understood, these invasive glioma cells prefer white matter and blood vessels [26,57] with a wide range of speeds in the range of 5-80 μm/h [58][59][60][61] showing sometimes saltatory patterns in the migration direction in brain [58]. Prediction of tumour invasion directions in nearby tissue may help define exact boundaries of focal treatments (surgery or radiosurgery), preventing future growth and recurrence [57]. For example, medical doctors could determine specific locations for radiation target volumes, not just using a rough estimate of 2-3 cm in all directions as a guidance as commonly done today [57].
Assuming that migratory cells are localized near the surgery site [16], one possible approach is to keep the cells in its proliferative phase preventing them from invading brain tissue in a combination with transport of therapeutic drugs near blood vessels [18]. As a result, the tumor Strategies in regulating glioblastoma signaling pathways and anti-invasion therapy mass will be visible for a succeeding surgery while killing proliferative tumor cells at the blood sites.
In this study, we considered the intracellular dynamics of the miR-451-AMPK-mTOR-cell cycle signaling pathway model developed recently by Kim et al. [40]. Incorporated in the model is a drug component which blocks the inhibitory pathway of mTOR by AMPK complex. This drug targets up-regulation of mTOR activities enhancing cell proliferation. The focus of the current work is to regulate up-stream signaling pathway via glucose infusion activating  Strategies in regulating glioblastoma signaling pathways and anti-invasion therapy miR-451, and control the down-stream pathway to cell cycle via drug infusion enhancing mTOR activities. Optimal control problem is formulated with the goal of keeping high levels of miR-451 and mTOR to induce cell growth and reduce invasion to the surrounding tissues.
In the framework of optimal control theory, two administration strategies are explored to achieve the goal with minimal cost incurred in glucose and drug administrations. The control strategies investigated in this study are (1) concomitant infusion control and, and (2) alternating glucose and drug infusion control. Both strategies are able to switch on the proliferative mode of glioblastoma cells and turn off its migratory mode. Cell cycle is regulated with fewer mass divisions restricting rapid growth. Numerical results show that concomitant control had fewer infusions, lesser glucose dosage and cheaper administration cost. However, when glucose and drug poses unwanted chemical reactions during concurrent administration, alternating control would be beneficial with lower drug amount usage. The mathematical models of the miR-451-AMPK-mTOR core control [14,[16][17][18]40] are based on a 'go-or-grow' hypothesis and supporting experiments in GBM cells [8,15,19]. However, the range of bistable behavior indicates a 'go-and-grow' program which has been proposed for other cancers too [62]. In the glioma cases, 'go-or-grow' hypothesis has been long suggested [29,55,[63][64][65][66] in addition to miR-451-induced 'go-or-grow' mechanism [8,15,19]. Glioma consists of a bulky, proliferative core in the center and highly invasive individual cells in the outer rim [55,63,64,67] and sequential transition between proliferative and invasive phenotypes characterizes tumor progression and may lead to faster growth [14]. At least at the microscopic level, these migration/proliferation phenotypes appear to be mutually exclusive characteristics at different time frames [64], as suggested by in vivo imaging data of glioma cells migrating in a saltatory fashion [58]. Glioma cells was shown to pause for a short period of time (�an hour) for cell division before the daughter cells begin to move again [58]. Gal et al. [68] further experimentally observed that this phenotypic switch can happen under different extracellular environment: (i) proliferative type in soft agar via activation of Myc signaling pathways in response to hepatocyte growth factor (HGF); and (ii) infiltrative type in Matrigel through Ras signaling pathways in response to the same HGF. Recently, Dhruv et al. [69] also affirmed the role of activation of key molecules in dichotomy between proliferation and invasion: c-Myc and NFκB in the proliferative core and radially dispersed, invasive region of GBM tumors, respectively. Glial cells can also interact with GBM tumor cells for phenotypic switch of GBM cells between cell migration and active proliferation [66]. Although these experiment and mathematical modeling support the 'go-or-grow' hypothesis, it is not certain if alternative mechanism such as 'go-andgrow' occurs in such a heterogenous TME in real patients [64]. Mansury et al. [70] for instance illustrated that individual glioma cell in a mixed group of proliferative and invasive phenotypes can depend on genotype of counter part and tumor microenvironment. In a similar conceptual studies on epithelial-to-mesenchymal transition (EMT) and mesenchymalto-epithelial transition (MET) [45] can give a hint on the complexities of this complex regulation of those two phenotypes and glioma in TME might not posses the simple 'go-or-grow' dogma. Increasing number of evidences now suggest that the 'go-or-grow' model has similar molecular basis with EMT [64]. For instance, HGF and TGFβ are major regulators of EMT, and these also provide strong stimuli of GBM invasion [71,72] and, upregulation of MET and CD44 activities, as well as an activated NFκB signaling pathway, was reported in both mesenchymal GBM and metastatic cancer [73][74][75][76]. All these experimental observations suggest that metastatic epithelial cancers and mesenchymal GBM drive common mechanisms that regulate the phenotypic transition between invasion and proliferation, suggesting the possibility of the 'go-and-grow' mechanism. The EMT paradigm in GBM has not been extensively studied as relevant to the progression of the disease due to different origin and rare metastasis of glioma [77,78]. Further studies and experimentation need to be done for better understanding of the fundamental principles.
In this paper, we did not take into account key microenvironmental factors such as endogenous immune dynamics including NK cells [79] and M1/M2 macrophages [72], other major signaling networks [80,81] such as E2F and Myc [11,82], angiogenesis [80,83], biophysical interaction between tumor cells and blood vessels [80], ECM remodeling for therapy [81,[84][85][86], or growth factors [87,88] such as epidermal growth factors [72,89,90], fibroblast growth factors [91], transforming growth factor-β [72,92], and CSF-1 [72,93,94], that may play critical roles in proliferation, progress, aggressive invasion of gliomas and development of anticancer strategies [95]. For example, endogenous NK cells may interfere oncolytic virus combination therapy in GBM while exogenous NK cells increase anti-tumor efficacy [79], illustrating the complex nature of tumor microenvironment. Recently, mTOR was considered to be a master regulator of cell growth and recognized as a good therapeutic target for therapies in glioblastoma [96]. A recent study found that withaferin A and temozolomide can induce apoptosis and reduce drug resistance by G2/M cell cycle arrest through intrinsic and extrinsic apoptotic signaling pathways [97]. More detailed analysis and experiments on Akt/mTOR/PI3K are necessary. Interestingly, radiation was shown to indirectly promote the export of lactate into the extracellular space and inhibition of AMPK/NFκB signaling pathways were involved in radiation-induced invasion of cancer cells [98]. On the other hand, M2 microglia/macrophages induce matrix remodeling and glioma cell invasion [72,[99][100][101]. Since PI3K signaling was shown to contribute to M2-polarization of these tumor associated macrophages (TAMs) [102] and PI3K binding to CSF1R was shown to enhance spreading of macrophages, thus promoting glioma cell invasion [103]. Therefore, better understanding of these PI3K-mTOR-CSF1R signaling networks in macrophages would lead to development of blocking aggressive infiltration of cancer cells. Signaling pathways of apoptosis and necroptosis are important parts of oncolytic virus (OVs) therapy [104,105]. Tumor extracellular matrix (ECM) plays a significant role in regulation of glioma invasion in brain tissue as well as OV spread [106]. For example, CSPGs, one of major tumor ECM in brain can characterize the invasive and non-invasive gliomas in a complex TME where microglia and astrocytes mechanically interact with CSPGs and tumor cells [40, [107][108][109]. Hybrid models [56,[110][111][112] and its associated optimal control strategies can be adapted to take into account this important intracellular signaling as well as cell population dynamics and tissue dynamics. Better understanding of various roles of these components in tumor microenvironment may provide better anti-invasion strategies of glioma cells.
However, our mathematical model in this work is a first step toward further experimental/ clinical investigation and more optimal anti-invasion strategies of GBM by incorporating these microenvironmental components. We will address these issues in future work.
Supporting information S1 Appendix. Optimal control problem. Formulation of the optimal control problem for concomitant and alternating glucose and drug infusion. (PDF)