Strategies of Eradicating Glioma Cells: A Multi-Scale Mathematical Model with MiR-451-AMPK-mTOR Control

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 and biomechanical cues in the microenvironment. Recent studies have shown that a particular microRNA, miR-451, regulates downstream molecules including AMPK and mTOR to determine the balance between rapid proliferation and invasion in response to metabolic stress in the harsh tumor microenvironment. Surgical removal of 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 multi-scale nature of glioblastoma proliferation and invasion and its response to conventional treatment, we propose a hybrid model of glioblastoma that analyses spatio-temporal dynamics at the cellular level, linking individual tumor cells with the macroscopic behaviour of cell organization and the microenvironment, and with the intracellular dynamics of miR-451-AMPK-mTOR signaling within a tumour cell. The model identifies a key mechanism underlying the molecular switches between proliferative phase and migratory phase in response to metabolic stress and biophysical interaction between cells in response to fluctuating glucose levels in the presence of blood vessels (BVs). The model predicts that cell migration, therefore efficacy of the treatment, not only depends on oxygen and glucose availability but also on the relative balance between random motility and strength of chemoattractants. Effective control of growing cells near BV sites in addition to relocalization of invisible migratory cells back to the resection site was suggested as a way of eradicating these migratory cells.

where g represents the signaling pathways from glucose to miR-451, 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. Λ 2 , Λ 4 , Λ 8 are the Hill-type inhibition saturation parameters from the counterpart 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. Finally µ 1 , µ 2 , µ 3 are microRNA/protein degradation rates of miR-451, AMPK complex, and mTOR, respectively. As indicated in equation (1), the signal g increases the rate of miR-451 activation through the function f (g), while AMPK-dependent inhibition of miR-451 is through the function F (a) in the denominator. A 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 equation (2)  . Proposed role of miR-451 in the regulation of LKB1/AMPK-mTOR signaling in response to high and low glucose levels [18]. miR-451 levels determine cell migration or proliferation in response to glucose (red triangle on the left) via the AMPK-mTOR network [19]. (A) Normal glucose levels up-regulate miR-451, which leads to down-regulation of the AMPK complex (CAB39/LKB1/AMPK). In turn, up-regulated AMPK levels induce up-regulation of mTOR, which leads to increased proliferation and decreased cell migration. (B) Low glucose levels reduce miR-451 levels, resulting in up-regulation of AMPK activity and down-regulation of mTOR. This leads to reduced cell proliferation and enhanced cell motility. Schematic components of miR-451, CAB39/LKB1/AMPK complex, and mTOR are represented by modules 'M ' (box with brown dotted line), 'A' (box with blue dotted line), and 'R' (box with brown dotted line on the right) in our theoretical framework.
activity induced by the signal s, and miR-451-dependent suppression of AMPK activity 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 equation (3) represents an increase in mTOR activity induced by the signal s 2 , and AMPK-dependent suppression of mTOR activity 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 ( Figure A), we assume that The governing equations (1)-(3) with the simple assumption (4) in a dimensional form become where all parameters (Λ G , g, s 1 , s 2 , Λ i (i = 1, . . . , 9), µ i (i = 1, 2, 3)) are non-negative constants. Figure B. A schematic of the core control model. (A) Conceptual model of regulation of miR451, AMPK complex, and mTOR in GBM cell migration and proliferation [19]. (B) Cartoon model (extended from [18]): miR-451 level and activity of its target complex (CAB39/LKB1/AMPK), and mTOR levels were represented by 'm', 'a', and 'r', respectively.

Mechanical effects on tumor growth: The cell-based component
Cell-mechanics is believed to play a large role in tumor growth and invasion. Mechanical stresses and the signaling pathways influence growth in a phenomenologically-specified manner in our model. We begin in the next section with a description of the cell-based component of the hybrid model.

The forces acting on individual cells
The mechanical behavior of individual cells is based on the model developed by Dallon and Othmer [1] (hereafter the paper and model is denoted DO) and Kim et al. [2] (hereafter the paper and model is denoted KSO). The new aspect that is needed in the present context is the core control system (miR-451-AMPK-mTOR) for cell proliferation and migration. The forces on a cell in the DO model include (i) the active forces exerted on neighboring cells or the substrate, (ii) the dynamic drag forces that arise as a moving cell forms and breaks adhesive bonds with neighboring cells, (iii) a static frictional force that exists when cells are rigidly attached to each other or to the substrate, and (iv) a reactive force due to forces exerted by other cells on it. The active force on cell i is denoted T i,j , wherein j = 0 denotes the substrate, and the reaction force to this is denoted M j,i . In the context of glioma migration via miR-451-AMPK-mTOR regulation, the specific form of active forces will be introduced in Section 2.1.3. The static force, which is denoted S j,i , is the binding force on the ith cell when bound to the jth. Since S i,j = −S j,i , the cell-cell forces cancel on all but those cells attached to the substrate. A more detailed discussion of all forces involved can be found in [1,2]. The total force on the ith cell is then summarized by where N a i denotes the neighbors of i, including the substrate, upon which it can exert traction, N d i is the set of cells (which includes substrate and extracellular matrix) that interact with i via a frictional force, and N s i denotes the set of cells that statically bind to cell i.

Cell growth and the rheology of the cytoplasm
There are two different kinds of cells involved in the system: proliferative cells and motile cells. Proliferative cells are modeled as in the KSO model [2] with intracellular dynamics (miR-451-AMPK-mTOR module) and an algorithm for migratory cells is introduced as in [1,3]. The KSO model, which is based on an earlier model for the mechanical behavior of cells and tissues under stress [1], addresses four major aspects (i) how an individual cell reacts to forces on it, (ii) how stress affects growth and inhibits growth, (iii) how cells interact mechanically with their surrounding microenvironment, and (iv) how growth and division can be described. The cells are treated as oriented ellipsoids and their cytoplasm considered as an incompressible, viscoelastic solid. In the absence of growth their volume is preserved (constant under all deformations), and growth is included in series with the active response and the passive forces. Changes in the length of an axis a of a cell consist of the passive change due to the length of the viscoelastic component (springdashpot module; u 0 a ) and the change due to the growth (u g a ) in response to given mechanical forces and biochemical signals (oxygen and glucose). The governing equations of the length of the i-th axis, i = a, b, c, of a cell are where u i is the change in the length of the ith axis, u 0 i and u g i are the changes in the length of the ith axis due to a change in the passive and growth element, respectively, f i is the magnitude of the force applied at each end, f 2 is the nonlinear spring force from the spring in parallel, k i is the spring constant for the spring in the Maxwell element, µ i is the viscous coefficient of the dashpot,p is the force due to pressure. See [1] for specific form of the function f 2 and details of how these equations are established. It is assumed that the passive response is incompressible [2,4], leading to the volume constraint for u 0 a , u 0 b , and u 0 c : where a * 0 , b * 0 and c * 0 are the lengths of three axes a, b, c after growth ( where a 0 , b 0 and c 0 are the initial lengths of three axes). Specific form of the growth termu g i is introduced below. In the KSO model [2] a cell can grow and divide into two daughter cells following the developed algorithm. This schematic algorithm incorporates the following major aspects (i) how to define an intrinsic cell-cycle time τ c ; (ii) how to define an intrinsic cell volume V 0 that cells attain immediately after division; (iii) how a cell divides under a given stress; (iv) how mechanical stress can inhibit cell growth. A cell relaxes to a spherical shape, whatever their initial shape, in the absence of external forces. In the present work we need to refine our assumptions on growth. In a previous work, Kim et al. [2,4] assumed that, in the absence of nutrient limitation, mechanical stress limitations [2,4], or other anti-growth signaling factors [4] cells grow to the volume 2V 0 and then instantly divide into equal two daughter cells. The orientation of cell division is determined by the direction of the net force exerted on the cell as in [2]. Several control mechanisms may delay or inhibit growth or cell cycle. Many different biomechanical and biochemical factors may affect the growth of normal cells and its regulation in cancerous cells may be very complex. (See the detailed discussion for regulation of cell growth below.) In the current model, we assume that the control mechanism of growth inhibition has three components. The first component controls the growth along each axis, which is zero when the force acting on the cell in that direction is too compressive or tensile. The second one is switch-like component, which is specific to the model used here, miR-451-AMPK-mTOR-dependent module. Finally, the third control component consists in cell cycle arrest and apoptosis due to chemotherapeutic that target either G 1 -or S-G 2 -M -phase in a cell cycle. We assume that the growth rate of proliferative cells depends on biomechanical signals (the mechanical stress acting on the cells) and biochemical signals (the up-or down-regulation of miR-451, AMPK, mTOR in the core control system in response to glucose) independently, and relate the growth of the axes of the cell to the tensions along major three axes. Thus we use the multiplicative form [3][4][5] of the growth rate function for the i-th axis given by where σ is the force acting on the cell and P is a function of the levels of miR-451 (M ), AMPK complex (A), and mTOR (R), and quiescent status ([G 0 ]) in the cell cycle. See Figure C. The growth function f (σ) is defined so that cells do not grow if forces are too large, but can grow under sufficiently small tensile and compressive forces. We also assume that the growth rate due to stress decays linearly with increasing compressive or tensile stress. Thus we define f (·) as follows (same as KSO with α = 0 there): where c + , c − are positive constants, As a result growth is either on or off in the context of mechanical stress acting on the cell. See Figure 4 in [2] for more detail. Cell proliferation depends on levels of many intracellular variables that control internal biological clocks of cell proliferation and migration. In the present work, we assume that the core control system and conditions of G 0 phase ([G 0 ] = 0 or 1) in the cell cycle module also determine cell proliferation, i.e., when the cell i is getting the right proliferation signal from the core control system (miR-451-AMPK-mTOR) and the cell is in the regular cell-cycle ([G 0 ] = 0). Thus, the function P (M, A, R, [G 0 ]) in the equation (12) is defined as where M, A, R is the levels of miR-451, AMPK complex, and mTOR, respectively, at the cell site and th M , th A , th R are the threshold values of M, A, R, respectively. See Figure C(B). Many factors including temporal changes [6] in chemicals or signaling levels, spatial differences [6,7], growth factors [8], or mechanical stress [8] can affect cell growth even in a simple system such as that in flies. Mechanical stresses can be processed through integrins-cytoskeleton and affect biochemical interactions [8][9][10][11][12][13], gene expression and overall growth [14,15]. Both compression and tension [8,[13][14][15] may promote cell division through the formation of mitotic spindles along the main axis of tension. This phenomenon of growth inhibition under mechanical constraints was also modeled with a simplification of those complicated relationships between growth and stresses in the context of growth of multilayer spheroids [2,4,16] and compared to experiments [14,15]. Mechanical interaction between a cell and substrate was also shown to govern growth behavior in epithelial cell growth [17]. The same hybrid model approach was proven to be useful in reproducing growth patterns of tumor epithelial cells under mechanical stresses in ductal carcinoma DCIS in situ [4,5].

Active force and equations of motion
For the force balance equations for each cell we take into consideration all the forces acting on the cell: the active forces, adhesion forces, internal pressure, and other forces as discussed in Section 2.1.1. We first describe the active force here. In DO model, the active force is generated by pulling on the neighboring cell whose center is closest to the line in the direction of desired motion when the cell is not in contact with the substrate. In the current model, we assume that migratory cells transmit the active force directly to the substrate and generate the active force by pulling the substrate. The model can easily be extended to incorporate active forces generated by all cells, but the cells in the interior of an aggregate need to connect to the surrounding tissue in order to do meaningful work toward moving the aggregate [1]. Furthermore, as we mentioned earlier, glioma cells in the interior of the growing tumor mass may grow but do not migrate. We assume that a glioma cell becomes motile (which is the case when the AMPK level stays above the threshold value (A > th A )), this generates an active force in the direction of a gradient consisting of a combination of glucose and chemoattractants, with a probability as well as in the direction from random motion as long as it is not surrounded by neighboring cells. The active force T i for migratory cell i is given by where d r is a unit vector of the moving direction from random motion, G, C are the concentrations of glucose and a chemoattractant, respectively (described in Materials and Methods Section in the main text), ψ 1 , ψ 2 , ψ 3 are scaling factors of weight distribution favoring random motion, glucose and other chemoattractants, respectively (ψ 1 , ψ 2 , ψ 3 ∈ [0, 1]; ψ 1 + ψ 2 + ψ 3 = 1), A is the level of AMPK complex at the ith cell site. We allow a small randomness in the magnitude of active force. Then, the indicator function φ(A) is given by where F 0 is the basal magnitude of the active force (0 ≤ |T i | ≤ F 0 ) and φ r is a random number (φ r ∈ [0.8, 1.2]). Therefore, the active force is completely turned off for cells in the proliferative phase (M > th A , A < th A , R > th R ) or cells under physical constraints. For example, when a cell is completely surrounded by neighboring cells, the active force is also turned off. From the formulation in the equation (15), it would be turned off when gradients of glucose (∇G) and chemoattractants (∇C) are zero in the absence of random motility (ψ 1 = 0), i.e., no active force is generated in the absence of chemotactic signals.
In summary, Newton's law for the ith cell reduces to where v i is the velocity of cell i, N i is the neighborhood of cell i, µ cell (resp., µ s , µ f ) is the degree of adhesiveness between the cells (resp., between the substrate and the cells, and the fluid viscosity), and  Table 2 in the main text. These mechanical components of the cell-based model, equations (9)- (17), are integrated with the reaction diffusion equations and the core control system in the hybrid framework in the main text. Figure D shows the overall interactions and flow of three major model components: miR-451-AMPK-mTOR system at the intracellular level, cell mechanics at the cellular level and microenvironmental factors such as glucose and oxygen.
(v) l K c (oxygen consumption rate) and µ K : Powanthil et al. [31] took the oxygen consumption rate of 0.2 s −1 based on the relation L = D K /φ where φ is rate of oxygen consumption by a cell at a cell site in the exponential decay term. In a similar fashion, we find (v) l G c (glucose consumption rate) : The measured level α = 1.6 pg/cell/min for piecewise increasing linear consumption term α(G)n was used in a glioma cell migration study [32] with a threshold value G th 1 = 2.0 × 10 −4 g/cm 3 , taken from the work of Li et al. [33]. We take l G c = 0.8 pg/cell/min and use a sphere with a radius r = 8 − 20 µm for estimation of cell volume.
(v) l 1 (ECM degradation rate by MMPs): It is hard to measure the decay rate of MMPs due to it's fast decay. Eisenberg et al. used the ECM degradation rate, 1.17×10 4 cm 3 g −1 s −1 in a myoferlin-mediated invasion model in breast cancer. We take the adjusted value l 1 =3.0×10 4 cm 3 g −1 s −1 from [34].
(iii) µ D (Decay rate of chemo-drugs): In a study [23] of effects of cell-cycle heterogeneity on the response of a tumour to chemotherapy, various decay rates of chemo-drugs were used from [30] Reference values: (i) G * : Sander and Deisboeck [32] used the characteristic concentration 2 × 10 −4 g/cm 3 of glucose with high values 6 × 10 −4 g/cm 3 in the far field (see also [46,47]). In in vitro experimental study by Godlewski et al. [19], low (0.3 g/l) and high (4.5 g/l) levels of glucose induced the up-and down-regulated expressions of miR-451 and mTOR (down-and up-regulation of AMPK complex), respectively. We take the high glucose level, G * = 4.5 × 10 −3 g/cm 3 , as a reference value.
(ii) C * : EGF is known to be chemoattractants in glioma cell lines [39]. We take the reference value C * = 1.0×10 −8 g/cm 3 from experimental studies of the role of EGF in the regulation of cancer progression [48,49].
(iv) P * : MMP concentrations were measured to be 1.6 × 10 −9 g/cm 3 in a breast cancer cell invasion study [51]. On the other hand, it was observed that PCK3145 may down-regulate MMP-9 expression in prostate cancer patients with up-regulated MMP-9 level of > 1.0 × 10 −7 g/cm 3 [52]. We take the bit high levels of MMPs: P * = 1.0 × 10 −7 g/cm 3 . Table 3 in the main text lists reference values and all parameter values above.
3.2 Core control system (miR-451, AMPK, mTOR) miRNAs are typically more stable than their targets [53,54] and the typical half-life of a miRNA is much larger than one of AMPK [55,56] and mTOR, leading to the small values of the fractions 1 = µ1 µ2 = 0.02, 2 = µ1 µ3 = 0.02 [3,18,57]. From the estimated concentrations of miRNAs (80 pM − 2.2 µM ) in an animal cell (assuming 1,000-25,000 µm 3 volume) [58], we take m * = 1.0 µM . The normal (high; 4.5 g/l) and low (0.3 g/l) glucose levels were used in [19] 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 [19] 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 [59], 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 [3,18]. By observing the same patterns of up-and down-regulation of mTOR as ones 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 [3,18]. 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 [3,18]. Once the parameters above determined, we fitted the data in the level of LKB1/AMPK activity in [19] in response to negative control and over expressed levels of miR-451 (Fig 5B in [19]) 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 equation (18) 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 [19]. 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 bit larger than the inhibition strength (β = 1.0) of the AMPK complex by miR-451 [3,18] while keeping the same value for the inhibition strength (γ = 1.0) of mTOR by the AMPK complex.
5 Dynamics of the system Figure E shows spatial profiles of concentrations of oxygen (K(x, t) in (A)), glucose (G(x, t) in (B)), ECM (ρ(x, t) in (C)), and MMP (P (x, t) in (D)) at corresponding times (t = 0, 80, 160, 240 h) as that of the spatial patterns of the tumour cells that are given in Figures 5A-5D in the main text. While invasive glioma cells after surgery undergo random motion (ψ 1 = 0.6) in the brain tissue, they are also attracted to BV sites via biochemical signals, in this case glucose and oxygen levels, which have highest values at fixed source points (BV sites; Figure