Next Article in Journal
The Management of Peutz–Jeghers Syndrome: European Hereditary Tumour Group (EHTG) Guideline
Next Article in Special Issue
Evolved Resistance to Placental Invasion Secondarily Confers Increased Survival in Melanoma Patients
Previous Article in Journal
Breathing Re-Education and Phenotypes of Sleep Apnea: A Review
Previous Article in Special Issue
Erratum: Iyer, A., et al. Integrative Analysis and Machine Learning Based Characterization of Single Circulating Tumor Cells. J. Clin. Med. 2020, 9, 1206
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Multi-Stability and Consequent Phenotypic Plasticity in AMPK-Akt Double Negative Feedback Loop in Cancer Cells

1
Department of Molecular Reproduction, Development, and Genetics, Indian Institute of Science, Bangalore 560012, India
2
Centre for BioSystems Science and Engineering, Indian Institute of Science, Bangalore 560012, India
3
Basic Sciences Division, Fred Hutchinson Cancer Research Center, Seattle, WA 98109, USA
*
Authors to whom correspondence should be addressed.
J. Clin. Med. 2021, 10(3), 472; https://doi.org/10.3390/jcm10030472
Submission received: 21 December 2020 / Revised: 7 January 2021 / Accepted: 21 January 2021 / Published: 26 January 2021

Abstract

:
Adaptation and survival of cancer cells to various stress and growth factor conditions is crucial for successful metastasis. A double-negative feedback loop between two serine/threonine kinases AMPK (AMP-activated protein kinase) and Akt can regulate the adaptation of breast cancer cells to matrix-deprivation stress. This feedback loop can significantly generate two phenotypes or cell states: matrix detachment-triggered pAMPKhigh/ pAktlow state, and matrix (re)attachment-triggered pAkthigh/ pAMPKlow state. However, whether these two cell states can exhibit phenotypic plasticity and heterogeneity in a given cell population, i.e., whether they can co-exist and undergo spontaneous switching to generate the other subpopulation, remains unclear. Here, we develop a mechanism-based mathematical model that captures the set of experimentally reported interactions among AMPK and Akt. Our simulations suggest that the AMPK-Akt feedback loop can give rise to two co-existing phenotypes (pAkthigh/ pAMPKlow and pAMPKhigh/pAktlow) in specific parameter regimes. Next, to test the model predictions, we segregated these two subpopulations in MDA-MB-231 cells and observed that each of them was capable of switching to another in adherent conditions. Finally, the predicted trends are supported by clinical data analysis of The Cancer Genome Atlas (TCGA) breast cancer and pan-cancer cohorts that revealed negatively correlated pAMPK and pAkt protein levels. Overall, our integrated computational-experimental approach unravels that AMPK-Akt feedback loop can generate multi-stability and drive phenotypic switching and heterogeneity in a cancer cell population.

1. Introduction

Despite major advances in cancer research, metastasis remains clinically unsolved and claims the vast majority of cancer-related deaths. Metastasis is a highly inefficient process with extremely high (>99.8%) rates of attrition. A hallmark of cancer cells that can successfully metastasize is their ability to dynamically adapt to their changing microenvironments, called phenotypic plasticity or switching [1,2]. Thus, understanding the rules of phenotypic plasticity and identifying therapeutic perturbations to reduce the fitness of metastasizing cells can be crucial for restricting the disease aggressiveness.
Metastasizing cells can often display phenotypic plasticity along multiple interconnected axes. Two of the most well-characterized axes are epithelial-mesenchymal plasticity (EMP) and cancer stem cell (CSC) plasticity [3,4]. More recently, the axes of metabolic plasticity (i.e., switching between more glycolytic vs. more oxidative phosphorylating states) [5,6] and drug resistance (i.e., switching between drug-resistant and drug-sensitive states) [7,8] are being investigated. A hallmark of networks involved in plasticity along these axes is the presence of double negative or mutually inhibitory feedback loops that can generate two (or more) phenotypes that cells can acquire. In the context of EMP, ZEB1 forms such loops with miR-200 and GRHL2, thus giving rise to multiple states–epithelial (high miR-200 and GRHL2, low ZEB), mesenchymal (low miR-200 and GRHL2, high ZEB), and hybrid epithelial/ mesenchymal (medium levels of miR-200, GRHL2 and ZEB) [9]. Similarly, for CSC plasticity, LIN28 and let-7 inhibit each other [10], and for metabolic plasticity, HIF-1α and AMPK (AMP-activated protein kinase) can inhibit each other [5]. The emergent dynamics of the above mentioned loops has been thoroughly investigated, and they have been shown as capable of exhibiting multi-stability (i.e., co-existence of multiple phenotypes). Such multi-stability has been posited to underlie phenotypic switching; disrupting such loops can restrict phenotypic switching, as witnessed for feedback loops of ZEB1 with GRHL2 and miR-200 [11,12].
Recent work from our laboratory has uncovered a double-negative feedback loop between two serine/threonine kinases AMPK and Akt operating in the adaptation of breast cancer cells to matrix-deprivation [13]. In epithelial cells, matrix-deprivation usually drives programmed cell death known as “anoikis” [14], but few detached epithelial cells can develop resistance to anoikis [15,16]. AMPK (AMP-activated protein kinase) is activated in cells facing bioenergetic or metabolic stress and can switch on energy-generating catabolic processes such as glycolysis and inhibit energy-consuming anabolic processes such as lipid and protein synthesis [17]. Conversely, upon growth factor stimulation, Akt becomes activated, promoting anabolic processes of lipid and protein synthesis, driving cell growth and proliferation [18]. Upon matrix deprivation, AMPK is activated which drives upregulation of PHLPP2 protein levels which can inactivate Akt. On the other hand, upon matrix (re)attachment, Akt is activated which can repress AMPK activity through PP2Cα [13]. Thus, while adherent cells showed a (re)attachment-triggered pAkthigh/pAMPKlow state, matrix-deprived cells demonstrated a detachment-triggered pAMPKhigh/pAktlow state. However, because this analysis was done at a population (or bulk) level, two questions remain to be answered: (a) can these cell states/phenotypes co-exist in the same cell population?, and (b) can these two subpopulations ‘spontaneously’ switch between themselves to give rise to one another?
Here, we adopt an integrative computational-experimental approach to answer these questions. First, we develop a mechanism-based mathematical model that captures the set of experimentally reported interactions among AMPK, Akt, PHLPP2 and PP2Cα. Simulations reveal that the AMPK–Akt feedback loop can give rise to two phenotypes–pAkthigh/ pAMPKlow, and pAMPKhigh/pAktlow–that can co-exist in specific parameter regimes, and switch between one another under the influence of biological noise. Next, we segregated the two subpopulations in MDA-MB-231 cells and observed that, under adherent conditions, each of them was capable of giving rise to another, thus validating our model prediction. Finally, clinical data analysis revealed a negative correlation between pAMPK and pAkt protein levels in The Cancer Genome Atlas (TCGA) breast cancer and pan-cancer cohorts. Overall, our results suggest that the AMPK–Akt feedback loop can be bistable, and therefore drive phenotypic switching and non-genetic heterogeneity in a cancer cell population.

2. Materials and Methods

2.1. ODE Model of The AMPK-Akt Network

The dynamics of the species in the regulatory network (Figure 2A) are represented using a system of Ordinary differential equations (ODEs) given below:
d d t A M P K     =   k a c A M P K   × t o t a l A M P K   A M P K k d a c A M P K × H n P P 2 C a , λ P P 2 c a , P P 2 c a 0 , P P 2 c α × A M P K )
d d t A K T   =   k a c A K T   ×   t o t a l A K T   A K T k d a c A K T   ×   H n P H L P P 2 , λ P H L P P 2 , P H L P P 2 0 , P H L P P 2   ×   A K T )
d d t P H L P P 2   =   k a c P H L P P 2   ×   t o t a l P H L P P 2   P H L P P 2   ×   H n A M P K , λ A M P K , A M P K 0 , A M P K k d a c P H L P P 2   ×   P H L P P 2 )
d d t P P 2 C a   =   k a c P P 2 C a   ×   t o t a l P P 2 C a   P P 2 C a   ×   H n A K T , λ A K T , A K T 0 , A K T k d a c P P 2 C a   ×   P P 2 C a )
Here, each equation represents the rate of change of the active levels of the entity given in the left-hand side of the equation. The total levels of the entities are constant and only switch between active and inactive states. The activation and deactivation rates of the species, given as kac(X) and kdac(X), (X ∈{AMPK,AKT,PHLPP2,PPCA}), respectively, are taken in units of t-1; all other parameters are in arbitrary units. The regulatory interactions are represented using shifted hill function form [19] as given below:
H n ,   λ , A 0 , A   =   1 +   λ ×   A A 0 n 1 + A A 0 n
where A is the effector species, n is the Hill’s coefficient that represents the non-linearity of the regulatory interaction, λ is the maximum fold change in activation/deactivation rate caused by A. λ > 1 indicates activation and λ < 1 indicates inhibition, A0 is the Activation/Inhibition threshold level of A. The parameter value ranges, with the corresponding dimensions, are given in Table 1 and Table S1. Assumptions made in the model are given below:
  • Total levels of each molecule are taken as a constant value of 100 arbitrary units (A.U.) and do not change throughout the simulation. However, the concentration of active and inactive molecular species can change.
  • Only the active state of the molecule affects another molecule’s conversions between its active and inactive forms.
  • Each molecule has its intrinsic activation and deactivation rate. The influence of interaction with other protein causing state changes is accounted by multiplying the corresponding rate term with a hill function.

2.2. Temporal Profiles and Steady State Estimation

10,000 Parameter sets were randomly sampled from the ranges mentioned in Table S1. For each parameter set, the ODE system was simulated using 1000 randomly generated initial conditions between the range (0–total molecules). The temporal profiles for each initial condition were computed numerically using the ode23 function of Matlab R2019A until 1000 time-steps. Initial conditions across distinct parameter sets reach steady state values by 1000 time-steps (Figure S1), thus the final state reported after 1000 time steps was considered to be the steady state value. For each parameter set, unique steady states of the system are taken from 1000 temporal profiles. Z-score calculation was performed across all the parameter sets for all the molecules. Based on the Z-score values of AMPK and Akt, parameter sets are grouped into four categories (LL, LH, HL, HH). Positive Z-score is considered as high state (H) and negative Z-score is considered as low state (L). A heatmap was generated based on the Z-scores obtained for all parameter sets using the ComplexHeatmap package in R [20].

2.3. Nullcline, Bifurcation and Phase Plane Analysis

For a given parameter set for the ODE system, nullclines for AMPK and Akt were obtained by setting their respective rates of change to zero and thus using constant value for the active level of that species. The corresponding steady state of the other three variables was obtained by temporal simulations of the modified ODE system as described above. Bifurcation analysis was done using the software package MATCONT [21]. Phase planes were obtained by combining multiple bifurcation diagrams.

2.4. Noise Induction

We used a simple formalism to induce noise which introduces a white noise sampled from a normal random distribution of mean 0 and variance η, η ∈ {20,30,40}. This noise is added to the species levels at fixed intervals (tstep = 100) in the simulation, representing the additive stochasticity to the species levels in the time period of simulation (ttotal = 5000). The fixed time interval is based on the observation that the mean time taken by different initial conditions across diverse parameter sets to attain equilibrium is around 100 timesteps (Figure S1B).

2.5. Clinical Data

Reverse Phase Protein Array (RPPA) dataset for TCGA- Pan Cancer 32, Breast Cancer and Sarcoma were downloaded from https://www.tcpaportal.org/tcpa/download.html. Pearson’s correlation analysis between phosphorylated AMPK (pT172) and phosphorylated AKT (pT308 and pS473) levels were calculated using cor.test() function from stats package and scatterplot was generated using plot() function in R 3.6.1.

2.6. Cell Line and Culture Condition, Fluorescence Activated Cell Sorting (FACS) Sorting and Analysis of The Plasticity

Breast cancer cell line MDA-MB-231 was procured from The American Type Culture Collection (ATCC) and validated subsequently by STR analysis. These cells were cultured in DMEM (Sigma-Aldrich) supplemented with 10% FBS (fetal bovine serum) containing penicillin and streptomycin, at 37 °C in 5% CO2 incubator. Cells were trypsinised and counted before seeding for every experiment. For FACS sorting purposes, MDA-MB-231 cells stably expressing EGR1promoter-TurboRFP (readout for the ERK activity; RFP: Red Fluorescent Protein) were sorted into high and low RFP cells after culturing on 90 mm tissue culture dish. These sorted cells were cultured in attached condition for the indicated time-points. Analysis of the high and low RFP sorted cells for their ability to attend the original heterogeneity was performed after the indicated time point of the culture in attached condition. Representative data were analysed using Summit software V5.2.1.12465 (Beckman Coulter, Miami, FL, USA).

2.7. Markov Chain Modelling and Simulations

To study the transition dynamics between the phenotypes observed in FACS data, the phenotypic transition is modelled as a Markov chain. In this formalism, a population of cells consists of 3 phenotypes, namely, the green phenotype (low EGR), red phenotype (high EGR) and black phenotype. Therefore, at any time t, the population is represented as:
F t = f R t ,   f G t , f B t
where fRt (fGt, fBt) is the fraction of the red (green, black) phenotype in the cell population at time t. Because the population size is large, these fractions are approximated to be the probabilities of a cell in the population showing the corresponding phenotype. The transition rates are defined as the conditional probabilities of a cell switching to a phenotype at time t+1 given its phenotype at time t. These transition rates build a transition matrix as follows:
T m =   P R | R P G | R P B | R P R | G P G | G P B | G P R | B P G | B P B | B
The transition matrix is assumed to be constant over time. Following the Markov chain property, the population at t+1 is obtained as follows:
F t + 1 = F t   ×   T m
The R package CellTrans is used to infer the transition matrix from the FACS data. Briefly, the package makes use of the above mentioned property of Markov chains and back calculates Tm using the information of Ft and F(t+1) from the FACS data. More details of the methodology are provided in [22].
Using the transition matrix, we simulated the trajectories of sorted cell populations. In each instance of the simulations, we start with a sorted population of size 10,000 cells. At each time step, the transition of a given cell from its current phenotype to a new phenotype is decided using a uniform random number and the row of the transition matrix corresponding to the current phenotype of the cell. 1000 such simulations were performed to obtain a distribution of the population fractions of cell phenotypes at each time point.

3. Results

3.1. AMPK-Akt Feedback Loop Can Give Rise To Two States: pAkthigh/ pAMPKlow and pAMPKhigh/pAktlo

First, we gathered experimentally curated information about interconnections among AMPK and Akt. AMPK and Akt can antagonistically regulate common downstream effectors such as mTOR signaling and FOXO signaling through differential phosphorylation [23]. Moreover, they can affect the activation status of one another. For instance, AMPK activating agents such as AICAR and phenformin can reduce the phosphorylation of Akt [24]. Adiponectin-activated AMPK can dephosphorylate Akt by increasing the activity of protein phosphatase 2A (PP2A) through dephosphorylating PP2Ac at Tyr307 [25]. On the other hand, upon insulin treatment, Akt is activated and it phosphorylates Ser487 of the serine/threonine rich loop (ST loop) in AMPK-α1 subunit, thus reducing subsequent phosphorylation and LKB1- or CAMKKβ-dependent AMPK activation at Thr172. Also, GSK3, another substrate of Akt, can phosphorylate the AMPK-α1 subunit at Thr481 and Ser477, further inhibiting AMPK activation by Thr172 phosphorylation [26]. Thus, AMPK and Akt pathways seem to inhibit the activity of one another.
Such a double negative association between AMPK and Akt was also reported in breast cancer cells during matrix attachment and detachment. The activation of AMPK upon matrix-deprivation drove upregulation of PHLPP2 protein levels, which can inactivate Akt. On the other hand, when cells were (re)attached to the matrix, Akt was activated which repressed AMPK activity through PP2Cα [13] (Figure 1). Put together, the abovementioned interactions reveal a mutually inhibitory feedback loop between the AMPK and Akt. This loop is reminiscent of “toggle switches” formed by various ‘master regulators’ of two (or more) diverse cell states, seen during embryonic development and disease progression [27]. Such feedback loops can occur at transcriptional [28,29], post-transcriptional [30,31], and cell-cell communication levels [32,33]. Here, a ‘toggle switch’ is observed between two kinases. For further analysis, we have focused on interactions known in the context of matrix-deprivation in breast cancer cells (Figure 2A).
Next, we investigated the emergent dynamics of the AMPK-Akt feedback loop. For the sake of simplicity, we considered the set of interactions reported for breast cancer cells during matrix-deprivation (Figure 2A): (a) AMPK and Akt can switch back and forth between their phosphorylated (active) and dephosphorylated (inactive) forms, (b) phosphorylated AMPK (pAMPK) can upregulate the levels of PHLPP2 which promote the dephosphorylation of Akt, and (c) phosphorylated Akt (pAkt) can upregulate the levels of PP2Cα associated with AMPK, thus enhancing AMPK dephosphorylation. These interactions are represented via a set of four coupled ordinary differential equations (ODEs). Each ODE tracks temporal evolution of the levels of AMPK, Akt, PHLPP2, PP2Cα, and the set of ODEs is solved numerically to obtain the steady state values for each of these four variables. To identify the robust dynamic features of this set of experimentally identified interactions, the kinetic parameters were chosen from a biologically relevant range of values (Table 1; see Section 2 (Materials and Methods)). We chose 10,000 such unique parameter sets to represent the effects of cell-to-cell heterogeneity and 1000 initial conditions for each parameter set to characterize all the possible phenotypes across parameter sets.
We collated the levels of AMPK, Akt, PHLPP2 and PP2Cα obtained from all parameter combinations and plotted them using a heatmap. We observed the emergence of four major clusters – pAMPKlow/pAktlow, pAMPKhigh/pAktlow, pAMPKlow/pAkthigh and pAMPKhigh/pAkthigh state with pAMPKhigh/ pAkthigh state being the least frequent and pAMPKlow/ pAktlow state being the most frequent (Figure 2B and Figure S2). A scatter plot between the pAMPK and pAKT levels revealed a significantly negative correlation (Figure 2C, Figure S2), suggesting that pAMPKhigh/pAktlow and pAMPKlow/pAkthigh states can be the dominant outputs of the network. The other nodes of the network also showed expected correlations trends (Table S2). However, the model also suggests the possible existence of pAMPKlow/pAktlow and pAMPKhigh/ pAkthigh states; the biological evidence for their existence still remains inconclusive. The occurrence of these four states associated with respective ratios of activation and deactivation of AMPK and Akt as identified from parameter sampling. The pAMPKlow/pAktlow state was observed in parameter cases where the net activation rate of AMPK or Akt was quite low, independent of the AMPK or Akt activities (i.e., the effects of PHLPP2 or PP2Cα) (Figure S3). Relatively low levels of pAMPK and/or pAKT have been reported in certain experimental reports [34,35]. However, here in the context of breast cancer, we focused our attention to pAMPKhigh/ pAktlow and pAMPKlow/pAkthigh states [13], and the parameter sets that converged to them. In particular, we focused on parameter sets that enabled a co-existence of these two states (bistability).

3.2. The Two States (pAkthigh/ pAMPKlow and pAMPKhigh/pAktlow) Can Co-Exist and Stochastically Switch Between One Another

To better understand bistability in this system, we performed nullcline and bifurcation analysis on the parameter sets showing bistability. First, for each parameter set, we constructed nullclines for pAMPK and pAKT. For a two-component system such as this, a nullcline represents the steady state levels of one component obtained over a range of values of the second component with its (i.e., the second component’s) rate of change over time being set to zero. The intersections of nullclines are, therefore, the points at which the rate of change for both the components is zero, i.e., a steady state of the system. The advantage of this approach over dynamical ODE-based simulations is that we can also identify unstable steady states that act as tipping points for transition from one stable steady state to another. In other words, for these parameter states, the two identified stable states can co-exist and can transition from one to another.
Here, for a given parameter set, we first calculated the steady state levels of active Akt (pAKT) for different constant values of active AMPK (pAMPK) (green curve in Figure 3A(I). Next, we calculated the steady state levels of pAMPK for different fixed levels of pAkt (red curve in Figure 3A(I)). These two curves are called as nullclines and their intersections identify the steady states of the system, two of which are stable (filled circles), and one unstable (hollow circle). The two stable states are pAkthigh/pAMPKlow and pAMPKhigh/pAktlow. The unstable state acts as a ‘tipping point’ beyond which a perturbation can allow switch from one state to another. Similar dynamical behavior was seen for other parametric combinations too, suggesting underlying bistability of the AMPK-Akt feedback loop (Figure 3A(II-III) and Figure S4).
To investigate how this feature of bistability may depend on various kinetic parameters in the model, we plotted a bifurcation diagram that tracks the levels of pAKT at varied values of activation rate of AMPK (kac AMPK). We chose kac AMPK as the bifurcation parameter to reflect the cases where the activation rate of AMPK can be altered by cell-intrinsic or cell-extrinsic factors. At low kac AMPK values (<0.05), pAkt levels are relatively higher, because the dephosphorylation of Akt by AMPK is weaker. Similarly, at high kac AMPK values (>0.13), pAkt levels are relatively lower, because of strong inactivation of Akt by AMPK. However, at intermediate range of these values (area bounded between dotted blue lines), a cell can exhibit bistability in terms of pAMPK levels, i.e., it can exist in either a pAkthigh/pAMPKlow or pAMPKhigh/pAktlow state (Figure 3B(I)). Similar trends are seen for other parameter sets shown earlier (Figure 3B(II-III)).
To quantify the range of bistability in the bifurcation diagram and its dependence on other kinetic parameters besides kac AMPK, we plotted this bifurcation at different values of activation rate of Akt (kac Akt) and observed that these curves (blue and magenta curves in Figure 3C(I) largely overlapped with that seen earlier (green curve in Figure 3C(I)). Next, we varied each parameter, one at a time, by +/− 10% and calculated the percentage change in the range of kac AMPK values enabling bistability, i.e., the distance between the dotted vertical lines. This sensitivity analysis suggested that the percentage change in this range was above 10% for only a few parameters such as those defining the effect of PP2Ca on AMPK (Figure 3C(II), Figures S5 and S6). These results underscore that the existence of bistability in AMPK-Akt feedback loop can be considered as largely robust to small parameter variations. Consequently, cells in an isogenic population may exist in two distinct states: pAkthigh/pAMPKlow and pAMPKhigh/pAktlow.
Next, we varied two parameters simultaneously, to map the two-dimensional parameter region in terms of (co-)existence of the two states. At high kac AMPK and low kac Akt, only the pAMPKhigh/pAktlow state was observed. Similarly, low kac AMPK and high kac Akt allowed only the pAkthigh/pAMPKlow state. Both these states co-existed in a parameter region lying between these extremes (yellow region in Figure 4A(I)). These trends were reinforced based on bifurcation diagrams drawn for an intermediate value of kac Akt (= 0.15), with kac AMPK as the bifurcation parameter. Bistability existed for intermediate values of kac AMPK; while higher or lower values led to monostable regions where only one state existed (Figure 4A(II)). Bifurcation diagram with kac Akt as the parameter confirmed the trend (Figure 4A(III)). Similar characteristics dynamics was seen for other parameter sets (Figure S7).
The co-existence of two (or more) states in a bifurcation diagram indicates that they may ‘spontaneously’ switch among one another. Thus, we performed stochastic simulations to examine for state switching under the influence of biological noise. These simulations revealed that cells can switch back and forth between the pAkthigh/ pAMPKlow and pAMPKhigh/pAktlow states for varying strengths of noise parameter (Figure 4B and Figure S8), highlighting that spontaneous state switching may be an outcome of the AMPK-Akt loop.
Such state switching implies that, when a cell population is sorted into subpopulations, it is possible for a subpopulation to give rise to another and potentially generate a population distribution similar to that seen in the parental population [36,37]. The rates of switching to and from a subpopulation to/from another one may be unequal, depending on the relative stability of the two (or more) states, evidenced by different mean residence times in a given state [38].

3.3. Experimental and Clinical Data Supports The Model Predictions of Bistability in AMPK-Akt Loop

To experimentally interrogate our observations of spontaneous state switching between pAkthigh/pAMPKlow and pAMPKhigh/pAktlow states, we performed FACS based cell sorting experiments in MDA-MB-231 cells stable for EGR1(promoter)-TurboRFP (EGR1 promoter-reporter system). In this approach, EGR1 promoter is used as a readout for AMPK activity, where high EGR activity (assessed by RFP intensity) corresponds to low AMPK activity and vice-versa [39]. Cells were grown in attached condition and sorted into high and low RFP group by FACS. These cells were cultured again for the indicated times in attached condition and we observed that these cells regained their original heterogeneous nature with time (Figure 5A and Figure S9A).
To quantify the observed phenotypic transitions between low and high AMPK (red and green populations respectively) cells, we constructed a discrete time Markov chain model [31]. The model assumes that transition rates between these cell types are independent of each other and do not change with time. We further assumed that the non-expressing cells (black population in the FACS distributions) do not transition into any other cell type. Using the R package CellTrans [22], we obtained the transition matrix, an ordered collection of transition rates (Figure 5B and Figure S9B). Using these transition rates, we simulated the evolution of population composition starting from homogeneous populations (Figure 5C and Figure S9C). The model predicts a higher transition rate from pAkthigh/pAMPKlow to pAMPKhigh/pAktlow phenotype, than that of pAMPKhigh/pAktlow to pAkthigh/pAMPKlow, suggesting that the pAMPKhigh/pAktlow population may be relatively more stable.
Further, to support our finding in the clinical setting, we performed a correlation of pAkt and pAMPK levels in patient samples using RPPA data from the TCGA cohort. A significant negative correlation was observed between pAkt and pAMPK levels across various cancer types (Figure 6, Table S4). This correlation was stronger for one of the two different phosphorylated versions of pAKT (compare top row vs. bottom row in Figure 6). These data further support the existence of bistability in the dynamics of AMPK-Akt feedback loop. Put together, our integrated computational-experiment analysis shows that the AMPK-Akt crosstalk can give rise to bistability and can lead to switching of states, driving heterogeneity at population level.

4. Discussion

Phenotypic switching can play crucial roles during cancer progression, as seen across cancer types. Prostate cancer cells can switch to a neuroendocrine-like state that is refractory to various therapies [40]. Similarly, in small cell lung cancer, cells under therapy-induced stress can reversibly switch to a hybrid neuroendocrine/mesenchymal state [41]. Besides these specific examples, EMP and CSCs are archetypal examples of phenotypic switching reported in many carcinomas [42,43] and non-epithelial tumors [44,45,46,47]. Recent progress in collecting high-throughput spatiotemporal data and mapping the regulatory networks underlying these axes of plasticity has led to developing complex mechanism-based and population-based mathematical models to decode dynamical traits of phenotypic switching such as dose-dependence, reversibility, hysteresis and transition rates among the cell states [48,49,50,51,52,53,54,55,56].
A hallmark of regulatory networks enabling phenotypic switching in cancer cell populations is multi-stability, i.e., the ability of isogenic cells to reversibly acquire diverse phenotypes [5,11,12,57,58,59,60,61], as reported earlier also for bacterial [62] and viral [63] populations. Multi-stability can enable ‘spontaneous’ switching among cell phenotypes (different attractors in the Waddington’s landscape) due to biological noise (that can operate at multiple levels including transcriptional or conformational [64,65]), and thus facilitate non-genetic heterogeneity [8,66]. For instance, in the context of EMP, PMC42-LA cells showed a bimodal distribution of EpCAM, and either subpopulation (EpCAM-high or EpCAM-low) was capable of generating the other without any exogenous overt induction [53]. Similar observations have been reported in maintaining a dynamic equilibrium of CSCs and non-CSCs in breast cancer [67].
Here, we demonstrate breast cancer cells switching between pAkthigh/pAMPKlow and pAMPKhigh/pAktlow states in adherent cell populations. Together with reinforcing observations of switching among these subpopulations in matrix-detached conditions using the same reporter construct in MDA-MB-231 cells [39], our results strongly support the existence of multi-stability in AMPK-Akt double negative feedback loop, as predicted by our mechanism-based mathematical model. A limitation of our mathematical model is the limited set of interactions between AMPK and Akt that we incorporated for the purpose of characterizing their dynamics in the context of matrix-deprivation. Other context-dependent interactions between AMPK and Akt have been reported, for instance, AMPK can activate Akt in acute lymphoblastic leukemia [68]. Also, AMPK mediated phosphorylation of Skp2 at S256 can activate the Skp2 SCF complex, driving K63-linked ubiquitination and eventual activation of Akt [69]. However, these interactions have not been yet observed in matrix-deprivation conditions. Such context-specific differences may allow for other dynamics for AMPK and Akt, such as oscillations seen in MCF10A cells upon inhibition of glycolysis and mitochondrial ATPase [70]. Intriguingly, AMPK can form double negative feedback loops with other molecules such as mTORC1 [71], which is involved in a similar feedback loop with ULK1 [72]. Coupling of such “toggle switches” can influence the emergence of multi-stability [73].
Another salient feature of multi-stability is hysteresis, as observed for bacterial cells [74] and for EMP in cancer cells [11]. Future experiments should investigate the possibility and implications of hysteresis in AMPK-Akt feedback loop driving the adaptation of breast cancer cells to matrix-deprivation (in other words, anchorage-independence) stress. Also, anchorage-independence has been shown to be associated with other axes of phenotypic plasticity: EMP [75,76,77], CSCs [78], and metabolic reprogramming [79,80,81]. However, a systems-level understanding of coordination of cell states along these interconnected axes remains elusive. Future integrative computational-experimental efforts, similar to the approach taken here, can be critical in investigating such coupled dynamics of phenotypic plasticity during the challenging metastatic cascade and identify therapeutic targets that can impact multiple axes of plasticity simultaneously.

Supplementary Materials

The following are available online at https://www.mdpi.com/2077-0383/10/3/472/s1, Table S1: Shows the literature reported values of parameter values. Table S2: Shows the Pearson’s correlation values and corresponding p values for the comparisons between AMPK, Akt, PHLPP2 and PP2Cα for three replicates, Table S3: Shows the 10 representative parameter sets used in the study, Table S4: Shows the Pearson’s correlation values and p values between AMPK (pT172) and Akt (pT308 and pS473) for 32 different TCGA cancer cohorts.

Author Contributions

M.K.J. conceived research; M.K.J. and A.R. supervised research; A.C., K.H. and S.K. performed research; all authors analyzed data and participated in writing and editing of the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

M.K.J. was supported by Dr. Vijaya and Rajgopal Rao Biomedical Research Fund awarded to the Centre for BioSystems Science and Engineering (BSSE), Indian Institute of Science (IISc), Bangalore. A.R. acknowledges support from DBT-IISc Partnership Program Phase-II (BT/PR27952-INF /22/212/2018), and support from DST-FIST and UGC, Govt. of India, to the department of MRDG, IISc. A.R. is a recipient of the Welcome Trust/DBT India Alliance Fellowship. A.C. acknowledges the Junior Research Fellowship (JRF) awarded by the Council of Scientific & Industrial Research (CSIR), India. K.H. was supported by Prime Ministers’ Research Fellowship (PMRF) awarded by the Ministry of Human Resources and Development (MHRD), Government of India.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The corresponding codes are available at https://github.com/csbBSSE/AmpkAkt. Reverse Phase Protein Array (RPPA) dataset for TCGA- Pan Cancer 32, Breast Cancer and Sarcoma were downloaded from https://www.tcpaportal.org/tcpa/download.html.

Acknowledgments

The authors acknowledge MRDG and IISc FACS facility.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Celià-Terrassa, T.; Kang, Y. Distinctive properties of metastasis- initiating cells. Genes Dev. 2016, 30, 892–908. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Welch, D.R.; Hurst, D.R. Defining the Hallmarks of Metastasis. Cancer Res. 2019, 79, 3011–3027. [Google Scholar] [CrossRef] [PubMed]
  3. Archana, P.T.; Saxena, K.; Murali, R.; Jolly, M.K.; Nair, R. Cancer Stem Cell Plasticity—A deadly deal. Front. Mol. Biosci. 2020, 7, 79. [Google Scholar]
  4. Celià-Terrassa, T.; Jolly, M.K. Cancer Stem Cells and Epithelial-to-Mesenchymal Transition in Cancer Metastasis. Cold Spring Harb. Perspect. Med. 2020, 10, a036905. [Google Scholar] [CrossRef]
  5. Jia, D.; Lu, M.; Jung, K.H.; Park, J.H.; Yu, L.; Onuchic, J.N.; Kaipparettu, B.A.; Levine, H. Elucidating cancer metabolic plasticity by coupling gene regulation with metabolic pathways. Proc. Natl. Acad. Sci. USA 2019, 116, 3909–3918. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Hanahan, D.; Weinberg, R.A. Hallmarks of cancer: The next generation. Cell 2011, 144, 646–674. [Google Scholar] [CrossRef] [Green Version]
  7. Jolly, M.K.; Kulkarni, P.; Weninger, K.; Orban, J.; Levine, H. Phenotypic plasticity, bet-hedging, and androgen independence in prostate cancer: Role of non-genetic heterogeneity. Front. Oncol. 2018, 8, 50. [Google Scholar] [CrossRef]
  8. Bell, C.C.; Gilan, O. Principles and mechanisms of non-genetic resistance in cancer. Br. J. Cancer 2020, 122, 465–472. [Google Scholar] [CrossRef]
  9. Jolly, M.K.; Tripathi, S.C.; Somarelli, J.A.; Hanash, S.M.; Levine, H. Epithelial/mesenchymal plasticity: How have quantitative mathematical models helped improve our understanding? Mol. Oncol. 2017, 11, 739–754. [Google Scholar] [CrossRef] [Green Version]
  10. Bocci, F.; Jolly, M.K.; George, J.T.; Levine, H.; Onuchic, J.N. A mechanism-based computational model to capture the interconnections among epithelial-mesenchymal transition, cancer stem cells and Notch-Jagged signaling. Oncotarget 2018, 9, 29906–29920. [Google Scholar] [CrossRef] [Green Version]
  11. Celià-Terrassa, T.; Bastian, C.; Liu, D.D.; Ell, B.; Aiello, N.M.; Wei, Y.; Zamalloa, J.; Blanco, A.M.; Hang, X.; Kunisky, D.; et al. Hysteresis control of epithelial-mesenchymal transition dynamics conveys a distinct program with enhanced metastatic ability. Nat. Commun. 2018, 9, 5005. [Google Scholar] [CrossRef] [PubMed]
  12. Hari, K.; Sabuwala, B.; Subramani, B.V.; La Porta, C.A.; Zapperi, S.; Font-Clos, F.; Jolly, M.K. Identifying inhibitors of epithelial-mesenchymal plasticity using a network topology based approach. NPJ Syst. Biol. Appl. 2020, 6, 15. [Google Scholar] [CrossRef] [PubMed]
  13. Saha, M.; Kumar, S.; Bukhari, S.; Balaji, S.A.; Kumar, P.; Hindupur, S.K.; Rangarajan, A. AMPK—Akt double-negative feedback loop in breast cancer cells regulates their adaptation to matrix deprivation. Cancer Res. 2018, 78, 1497–1510. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Frisch, S.M.; Francis, H. Disruption of epithelial cell-matrix interactions induces apoptosis. J. Cell Biol. 1994, 124, 619. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Guadamillas, M.C.; Cerezo, A.; del Pozo, M.A. Overcoming anoikis—Pathways to anchorage-independent growth in cancer. J. Cell Sci. 2011, 124, 3189–3197. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  16. Buchheit, C.L.; Weigel, K.J.; Schafer, Z.T. Cancer cell survival during detachment from the ECM: Multiple barriers to tumour progression. Nat. Rev. Cancer 2014, 14, 632–641. [Google Scholar] [CrossRef]
  17. Hardie, D.G.; Schaffer, B.E.; Brunet, A. AMPK: An Energy-Sensing Pathway with Multiple Inputs and Outputs. Trends Cell Biol. 2016, 26, 190–201. [Google Scholar] [CrossRef] [Green Version]
  18. Hoxhaj, G.; Manning, B.D. The PI3K–AKT network at the interface of oncogenic signalling and cancer metabolism. Nat. Rev. Cancer 2020, 20, 74–88. [Google Scholar] [CrossRef]
  19. Lu, M.; Jolly, M.K.; Levine, H.; Onuchic, J.N.; Ben-Jacob, E. MicroRNA-based regulation of epithelial-hybrid-mesenchymal fate determination. Proc. Natl. Acad. Sci. USA 2013, 110, 18174–18179. [Google Scholar] [CrossRef] [Green Version]
  20. Gu, Z.; Eils, R.; Schlesner, M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 2016, 32, 2847–2849. [Google Scholar] [CrossRef] [Green Version]
  21. Dhooge, A.; Govaerts, W.; Kuznetsov, Y.A.; Meijer, H.G.E.; Sautois, B. New features of the software MatCont for bifurcation analysis of dynamical systems. Math. Comput. Model. Dyn. Syst. 2008, 14, 147–175. [Google Scholar] [CrossRef]
  22. Buder, T.; Deutsch, A.; Seifert, M.; Voss-Böhme, A. CellTrans: An R Package to Quantify Stochastic Cell State Transitions. Bioinform. Biol. Insights 2017, 11, 1177932217712241. [Google Scholar] [CrossRef] [PubMed]
  23. Zhao, Y.; Hu, X.; Liu, Y.; Dong, S.; Wen, Z.; He, W.; Zhang, S.; Huang, Q.; Shi, M. ROS signaling under metabolic stress: Cross-talk between AMPK and AKT pathway. Mol. Cancer 2017, 16, 79. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. King, T.D.; Song, L.; Jope, R.S. AMP-activated protein kinase (AMPK) activating agents cause dephosphorylation of Akt and glycogen synthase kinase-3. Biochem. Pharmacol. 2006, 71, 1637–1647. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Kim, K.Y.; Baek, A.; Hwang, J.E.; Choi, Y.A.; Jeong, J.; Lee, M.S.; Cho, D.H.; Lim, J.S.; Kim, K., II; Yang, Y. Adiponectin-activated AMPK stimulates dephosphorylation of AKT through protein phosphatase 2A activation. Cancer Res. 2009, 69, 4018–4026. [Google Scholar] [CrossRef] [Green Version]
  26. Hawley, S.A.; Ross, F.A.; Gowans, G.J.; Tibarewal, P.; Leslie, N.R.; Hardie, D.G. Phosphorylation by Akt within the ST loop of AMPK-α1 down-regulates its activation in tumour cells. Biochem. J. 2014, 459, 275–287. [Google Scholar] [CrossRef] [Green Version]
  27. Zhou, J.X.; Huang, S. Understanding gene circuits at cell-fate branch points for rational cell reprogramming. Trends Genet. 2011, 27, 55–62. [Google Scholar] [CrossRef]
  28. Guantes, R.; Poyatos, J.F. Multistable decision switches for flexible control of epigenetic differentiation. PLoS Comput. Biol. 2008, 4, 1000235. [Google Scholar] [CrossRef] [Green Version]
  29. Varankar, S.S.; Kamble, S.C.; Mali, A.M.; More, M.M.; Abraham, A.; Kumar, B.; Pansare, K.J.; Narayanan, N.J.; Sen, A.; Dhake, R.D.; et al. Functional balance between TCF21-Slug defines cellular plasticity and sub-classes in high-grade serous ovarian cancer. Carcinogenesis 2020, 17, 515–526. [Google Scholar] [CrossRef]
  30. Jolly, M.K.; Preca, B.-T.; Tripathi, S.C.; Jia, D.; George, J.T.; Hanash, S.M.; Brabletz, T.; Stemmler, M.P.; Maurer, J.; Levine, H. Interconnected feedback loops among ESRP1, HAS2, and CD44 regulate epithelial-mesenchymal plasticity in cancer. APL Bioeng. 2018, 2, 031908. [Google Scholar] [CrossRef]
  31. Yang, X.; Lin, X.; Zhong, X.; Kaur, S.; Li, N.; Liang, S.; Lassus, H.; Wang, L.; Katsaros, D.; Montone, K.; et al. Double-Negative Feedback Loop between Reprogramming Factor LIN28 and microRNA let-7 Regulates Aldehyde Dehydrogenase 1-Positive Cancer Stem Cells. Cancer Res. 2010, 70, 9463–9472. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  32. Jolly, M.K.; Boareto, M.; Lu, M.; Onuchic, J.N.; Clementi, C.; Ben-Jacob, E. Operating principles of Notch-Delta-Jagged module of cell-cell communication. N. J. Phys. 2015, 17, 055021. [Google Scholar] [CrossRef] [Green Version]
  33. Sprinzak, D.; Lakhanpal, A.; Lebon, L.; Santat, L.A.; Fontes, M.E.; Anderson, G.A.; Garcia-Ojalvo, J.; Elowitz, M.B. Cis-interactions between Notch and Delta generate mutually exclusive signalling states. Nature 2010, 465, 86–90. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Zhande, R.; Zhang, W.; Zheng, Y.; Pendleton, E.; Li, Y.; Polakiewicz, R.D.; Sun, X.J. Dephosphorylation by Default, a Potential Mechanism for Regulation of Insulin Receptor Substrate-1/2, Akt, and ERK1/2. J. Biol. Chem. 2006, 281, 39071–39080. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Willows, R.; Sanders, M.J.; Xiao, B.; Patel, B.R.; Martin, S.R.; Read, J.; Wilson, J.R.; Hubbard, J.; Gamblin, S.J.; Carling, D. Phosphorylation of AMPK by upstream kinases is required for activity in mammalian cells. Biochem. J. 2017, 474, 3059–3073. [Google Scholar] [CrossRef] [Green Version]
  36. Huang, S. Non-genetic heterogeneity of cells in development: More than just noise. Development 2009, 136, 3853–3862. [Google Scholar] [CrossRef] [Green Version]
  37. Gupta, P.B.; Fillmore, C.M.; Jiang, G.; Shapira, S.D.; Tao, K.; Kuperwasser, C.; Lander, E.S. Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells. Cell 2011, 146, 633–644. [Google Scholar] [CrossRef] [Green Version]
  38. Biswas, K.; Jolly, M.; Ghosh, A. Stability and mean residence times for hybrid epithelial/mesenchymal phenotype. Phys. Biol. 2019, 16, 025003. [Google Scholar] [CrossRef]
  39. Kumar, S.; Hari, K.; Jolly, M.K.; Rangarajan, A. Feedback loops involving ERK, AMPK and TFEB generate non-genetic heterogeneity that allows cells to evade anoikis. bioRxiv 2019, 736546. [Google Scholar] [CrossRef]
  40. Hu, C.-D.; Choo, R.; Huang, J. Neuroendocrine Differentiation in Prostate Cancer: A Mechanism of Radioresistance and Treatment Failure. Front. Oncol. 2015, 5, 90. [Google Scholar] [CrossRef] [Green Version]
  41. Udyavar, A.A.; Wooten, D.J.; Hoeksema, M.; Bansal, M.; Califano, M.; Estrada, L.; Schnell, S.; Irish, J.M.; Massion, P.P.; Quaranta, V. Novel Hybrid Phenotype Revealed in Small Cell Lung Cancer by a Transcription Factor Network Model That Can Explain Tumor Heterogeneity. Cancer Res. 2017, 77, 1063–1074. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Scheel, C.; Weinberg, R.A. Cancer stem cells and epithelial-mesenchymal transition: Concepts and molecular links. Semin. Cancer Biol. 2012, 22, 396–403. [Google Scholar] [CrossRef] [PubMed]
  43. Jolly, M.K.; Celià-Terrassa, T. Dynamics of Phenotypic Heterogeneity Associated with EMT and Stemness during Cancer Progression. J. Clin. Med. 2019, 8, 1542. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Kahlert, U.D.; Joseph, J.V.; Kruyt, F.A.E. EMT- and MET-related processes in nonepithelial tumors: Importance for disease progression, prognosis, and therapeutic opportunities. Mol. Oncol. 2017, 11, 860–877. [Google Scholar] [CrossRef] [PubMed]
  45. Somarelli, J.A.; Shetler, S.; Jolly, M.K.; Wang, X.; Bartholf Dewitt, S.; Hish, A.J.; Gilja, S.; Eward, W.C.; Ware, K.E.; Levine, H.; et al. Mesenchymal-epithelial transition in sarcomas is controlled by the combinatorial expression of microRNA 200s and GRHL2. Mol. Cell. Biol. 2016, 36, 2503–2513. [Google Scholar] [CrossRef] [Green Version]
  46. Genadry, K.C.; Pietrobono, S.; Rota, R.; Linardic, C.M. Soft tissue sarcoma cancer stem cells: An overview. Front. Oncol. 2018, 8, 475. [Google Scholar] [CrossRef] [PubMed]
  47. Bonnet, D.; Dick, J.E. Human acute myeloid leukemia is organized as a hierarchy that originates from a primitive hematopoietic cell. Nat. Med. 1997, 3, 730–737. [Google Scholar] [CrossRef]
  48. Karacosta, L.G.; Anchang, B.; Ignatiadis, N.; Kimmey, S.C.; Benson, J.A.; Shrager, J.B.; Tibshirani, R.; Bendall, S.C.; Plevritis, S.K. Mapping Lung Cancer Epithelial-Mesenchymal Transition States and Trajectories with Single-Cell Resolution. Nat. Commun. 2019, 10, 5587. [Google Scholar] [CrossRef] [Green Version]
  49. Enderling, H. Cancer stem cells: Small subpopulation or evolving fraction? Integr. Biol. 2015, 7, 14–23. [Google Scholar] [CrossRef] [Green Version]
  50. Devaraj, V.; Bose, B. Morphological State Transition Dynamics in EGF-Induced Epithelial to Mesenchymal Transition. J. Clin. Med. 2019, 8, 911. [Google Scholar] [CrossRef] [Green Version]
  51. Mandal, M.; Ghosh, B.; Anura, A.; Mitra, P.; Pathak, T.; Chatterjee, J. Modeling continuum of epithelial mesenchymal transition plasticity. Integr. Biol. 2016, 8, 167–176. [Google Scholar] [CrossRef] [PubMed]
  52. Zhang, J.; Tian, X.-J.; Zhang, H.; Teng, Y.; Li, R.; Bai, F.; Elankumaran, S.; Xing, J. TGF-β–induced epithelial-to-mesenchymal transition proceeds through stepwise activation of multiple feedback loops. Sci. Signal. 2014, 7, 91. [Google Scholar] [CrossRef] [PubMed]
  53. Bhatia, S.; Monkman, J.; Blick, T.; Pinto, C.; Waltham, A.; Nagaraj, S.H.; Thompson, E.W. Interrogation of Phenotypic Plasticity between Epithelial and Mesenchymal States in Breast Cancer. J. Clin. Med. 2019, 8, 893. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  54. Watanabe, K.; Panchy, N.; Noguchi, S.; Suzuki, H.; Hong, T. Combinatorial perturbation analysis reveals divergent regulations of mesenchymal genes during epithelial-to-mesenchymal transition. NPJ Syst. Biol. Appl. 2019, 5, 21. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  55. Steinway, S.N.; Zañudo, J.G.T.; Michel, P.J.; Feith, D.J.; Loughran, T.P.; Albert, R. Combinatorial interventions inhibit TGFβ-driven epithelial-to-mesenchymal transition and support hybrid cellular phenotypes. NPJ Syst. Biol. Appl. 2015, 1, 15014. [Google Scholar] [CrossRef] [Green Version]
  56. Wang, W.; Douglas, D.; Zhang, J.; Kumari, S.; Enuameh, M.S.; Dai, Y.; Wallace, C.T.; Watkins, S.C.; Shu, W.; Xing, J. Live-cell imaging and analysis reveal cell phenotypic transition dynamics inherently missing in snapshot data. Sci. Adv. 2020, 6. [Google Scholar] [CrossRef]
  57. Tripathi, S.; Chakraborty, P.; Levine, H.; Jolly, M.K. A mechanism for epithelial-mesenchymal heterogeneity in a population of cancer cells. PLoS Comput. Biol. 2020, 16, 1007619. [Google Scholar] [CrossRef] [Green Version]
  58. Goetz, H.; Melendez-Alvarez, J.R.; Chen, L.; Tian, X.-J. A plausible accelerating function of intermediate states in cancer metastasis. PLoS Comput. Biol. 2020, 16, 1007682. [Google Scholar] [CrossRef]
  59. Lee, J.; Lee, J.; Farquhar, K.S.; Yun, J.; Frankenberger, C.A.; Bevilacqua, E.; Yeung, K.; Kim, E.-J.; Balázsi, G.; Rosner, M.R. Network of mutually repressive metastasis regulators can promote cell heterogeneity and metastatic transitions. Proc. Natl. Acad. Sci. USA 2014, 111, 364–373. [Google Scholar] [CrossRef] [Green Version]
  60. Xin, Y.; Cummins, B.; Gedeon, T. Multistability in the epithelial-mesenchymal transition network. BMC Bioinform. 2020, 21, 71. [Google Scholar] [CrossRef] [Green Version]
  61. Huang, B.; Jolly, M.K.; Lu, M.; Tsarfaty, I.; Onuchic, J.N.; Ben-Jacob, E. Modeling the transitions between collective and solitary migration phenotypes in cancer metastasis. Sci. Rep. 2015, 5, 17379. [Google Scholar] [CrossRef] [PubMed]
  62. Balaban, N.Q.; Merrin, J.; Chait, R.; Kowalik, L.; Leibler, S. Bacterial Persistence as a Phenotypic Switch. Science 2004, 305, 1622–1625. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  63. Padmanabhan, P.; Garaigorta, U.; Dixit, N.M. Emergent properties of the interferon-signalling network may underlie the success of hepatitis C treatment. Nat. Commun. 2014, 5, 3872. [Google Scholar] [CrossRef] [Green Version]
  64. Elowitz, M.B.; Levine, A.J.; Siggia, E.D.; Swain, P.S. Stochastic gene expression in a single cell. Science 2002, 297, 1183–1186. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  65. Mooney, S.M.; Jolly, M.K.; Levine, H.; Kulkarni, P. Phenotypic plasticity in prostate cancer: Role of intrinsically disordered proteins. Asian J. Androl. 2016, 18, 704–710. [Google Scholar] [CrossRef]
  66. Brock, A.; Chang, H.; Huang, S. Non-genetic heterogeneity—A mutation-independent driving force for the somatic evolution of tumours. Nat. Rev. Genet. 2009, 10, 336–342. [Google Scholar] [CrossRef]
  67. Iliopoulos, D.; Hirsch, H.A.; Wang, G.; Struhl, K. Inducible formation of breast cancer stem cells and their dynamic equilibrium with non-stem cancer cells via IL6 secretion. Proc. Natl. Acad. Sci. USA 2011, 108, 1397–1402. [Google Scholar] [CrossRef] [Green Version]
  68. Leclerc, G.M.; Leclerc, G.J.; Fu, G.; Barredo, J.C. AMPK-induced activation of Akt by AICAR is mediated by IGF-1R dependent and independent mechanisms in acute lymphoblastic leukemia. J. Mol. Signal. 2010, 5, 15. [Google Scholar] [CrossRef] [Green Version]
  69. Han, F.; Li, C.F.; Cai, Z.; Zhang, X.; Jin, G.; Zhang, W.N.; Xu, C.; Wang, C.Y.; Morrow, J.; Zhang, S.; et al. The critical role of AMPK in driving Akt activation under stress, tumorigenesis and drug resistance. Nat. Commun. 2018, 9, 4728. [Google Scholar] [CrossRef]
  70. Hung, Y.P.; Teragawa, C.; Kosaisawe, N.; Gillies, T.E.; Pargett, M.; Minguet, M.; Distor, K.; Rocha-Gregg, B.L.; Coloff, J.L.; Keibler, M.A.; et al. Akt regulation of glycolysis mediates bioenergetic stability in epithelial cells. eLife 2017, 6, 27293. [Google Scholar] [CrossRef]
  71. Holczer, M.; Hajdú, B.; Lőrincz, T.; Szarka, A.; Bánhegyi, G.; Kapuy, O. A double negative feedback loop between MTORC1 and AMPK kinases guarantees precise autophagy induction upon cellular stress. Int. J. Mol. Sci. 2019, 20, 5543. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  72. Szymanska, P.; Martin, K.R.; MacKeigan, J.P.; Hlavacek, W.S.; Lipniacki, T. Computational analysis of an autophagy/translation switch based on mutual inhibition of MTORC1 and ULK1. PLoS ONE 2015, 10, 0116550. [Google Scholar] [CrossRef] [PubMed]
  73. Duddu, A.S.; Sahoo, S.; Hati, S.; Jhunjhunwala, S.; Jolly, M.K. Multi-stability in cellular differentiation enabled by a network of three mutually repressing master regulators. J. R. Soc. Interface 2020, 17, 20200631. [Google Scholar] [CrossRef] [PubMed]
  74. Ozbudak, E.M.; Thattai, M.; Lim, H.N.; Shraiman, B.I.; van Oudenaarden, A. Multistability in the lactose utilization network of Escherichia coli. Nature 2004, 427, 737–740. [Google Scholar] [CrossRef] [PubMed]
  75. Kumar, S.; Park, S.H.; Cieply, B.; Schupp, J.; Killiam, E.; Zhang, F.; Rimm, D.L.; Frisch, S.M. A Pathway for the Control of Anoikis Sensitivity by E-Cadherin and Epithelial-to-Mesenchymal Transition. Mol. Cell. Biol. 2011, 31, 4036–4051. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  76. Huang, R.Y.-J.; Wong, M.K.; Tan, T.Z.; Kuay, K.T.; Ng, A.H.C.; Chung, V.Y.; Chu, Y.-S.; Matsumura, N.; Lai, H.-C.; Lee, Y.F.; et al. An EMT spectrum defines an anoikis-resistant and spheroidogenic intermediate mesenchymal state that is sensitive to e-cadherin restoration by a src-kinase inhibitor, saracatinib (AZD0530). Cell Death Dis. 2013, 4, 915. [Google Scholar] [CrossRef] [Green Version]
  77. Jolly, M.K.; Ware, K.E.; Xu, S.; Gilja, S.; Shetler, S.; Yang, Y.; Wang, X.; Austin, R.G.; Runyambo, D.; Hish, A.J.; et al. E-Cadherin Represses Anchorage-Independent Growth in Sarcomas through Both Signaling and Mechanical Mechanisms. Mol. Cancer Res. 2019, 17, 1391–1402. [Google Scholar] [CrossRef] [Green Version]
  78. Kim, S.Y.; Hong, S.H.; Basse, P.H.; Wu, C.; Bartlett, D.L.; Kwon, Y.T.; Lee, Y.J. Cancer Stem Cells Protect Non-Stem Cells From Anoikis: Bystander Effects. J. Cell. Biochem. 2016, 2301, 2289–2301. [Google Scholar] [CrossRef]
  79. Zhang, P.; Song, Y.; Sun, Y.; Li, X.; Chen, L.; Yang, L. AMPK/GSK3b/b-catenin cascade—Triggered overexpression of CEMIP promotes migration and invasion in anoikis-resistant prostate cancer cells by enhancing metabolic reprogramming. FASEB J. 2018, 32, 3924–3935. [Google Scholar] [CrossRef] [Green Version]
  80. Kamarajugadda, S.; Stemboroski, L.; Cai, Q.; Simpson, N.E.; Nayak, S.; Tan, M.; Lu, J. Glucose Oxidation Modulates Anoikis and Tumor Metastasis. Mol. Cell. Biol. 2012, 32, 1893–1907. [Google Scholar] [CrossRef] [Green Version]
  81. Caneba, C.A.; Bellance, N.; Yang, L.; Pabst, L.; Nagrath, D. Pyruvate uptake is increased in highly invasive ovarian cancer cells under anoikis conditions for anaplerosis, mitochondrial function, and migration. Am. J. Physiol. Endocrinol. Metab. 2012, 303, 1036–1052. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Figure 1. AMPK-Akt feedback loop. Regulatory network between AMPK and Akt. Green arrows and associated green ‘P’ circles denote activation by phosphorylation. Red arrows indicate deactivation by phosphorylation, and red hammerheads show deactivation by dephosphorylation. Solid lines show known molecular mechanisms, and dashed lines show scenarios where such information is not available.
Figure 1. AMPK-Akt feedback loop. Regulatory network between AMPK and Akt. Green arrows and associated green ‘P’ circles denote activation by phosphorylation. Red arrows indicate deactivation by phosphorylation, and red hammerheads show deactivation by dephosphorylation. Solid lines show known molecular mechanisms, and dashed lines show scenarios where such information is not available.
Jcm 10 00472 g001
Figure 2. Dynamics of AMPK-Akt feedback loop. (A) Reduced network considered in the model containing a total of eight species – each of the four molecules (AMPK, Akt, PHLPP2, PP2Cα) has an active (rectangle) and an inactive (oval) form. Black arrows show the conversion of species from active to inactive and vice-versa. Red hammerheads represent inhibition, green arrows represent activation. Unless stated otherwise, further analysis shows the levels of only active species of the four molecules. (B) Heatmap of steady states attained by 10,000 random parameter sets generated from 1000 random initial conditions of the active levels of AMPK, Akt, PHLPP2, and PP2Cα. Color is based on the z-score calculated for the whole set of simulations across all parameter sets, orange represents positive z-score (high) and purple represents negative z-score (low). LL, HL, LH and HH denote the four states when considering the steady state levels of pAMPK, pAKT - pAMPKlow/ pAktlow, pAMPKhigh/pAktlow, pAMPKlow/ pAkthigh and pAMPKhigh/pAkthigh. (C) Scatter plot of z-scores of AMPK and Akt steady state values represented in heatmap showing the state distribution. Pearson correlation coefficient, p-value are reported.
Figure 2. Dynamics of AMPK-Akt feedback loop. (A) Reduced network considered in the model containing a total of eight species – each of the four molecules (AMPK, Akt, PHLPP2, PP2Cα) has an active (rectangle) and an inactive (oval) form. Black arrows show the conversion of species from active to inactive and vice-versa. Red hammerheads represent inhibition, green arrows represent activation. Unless stated otherwise, further analysis shows the levels of only active species of the four molecules. (B) Heatmap of steady states attained by 10,000 random parameter sets generated from 1000 random initial conditions of the active levels of AMPK, Akt, PHLPP2, and PP2Cα. Color is based on the z-score calculated for the whole set of simulations across all parameter sets, orange represents positive z-score (high) and purple represents negative z-score (low). LL, HL, LH and HH denote the four states when considering the steady state levels of pAMPK, pAKT - pAMPKlow/ pAktlow, pAMPKhigh/pAktlow, pAMPKlow/ pAkthigh and pAMPKhigh/pAkthigh. (C) Scatter plot of z-scores of AMPK and Akt steady state values represented in heatmap showing the state distribution. Pearson correlation coefficient, p-value are reported.
Jcm 10 00472 g002
Figure 3. Nullcline, bifurcation and sensitivity analysis. (A) Nullclines for three representative bistable parameter sets (rows 1–3 in Table S3). Green curve is AMPK Nullcline (d/dt (AMPK)=0). Red curve is Akt Nullcline (d/dt (AKT)=0). Blue circles represent stable states, white circles represent unstable steady state. (B) Bifurcation of Akt steady state levels with respect to the activation rate of AMPK (k_ac AMPK) for three representative bistable parameter sets (rows 1–3 in Table S3). Green curves denote stable states and red curves denote unstable states. Region bound by the blue dashed line represents the bistable region, where both states can co-exist. (C) (I) Representative bifurcation of Akt levels with respect to activation rate of AMPK (k_ac AMPK), drawn for three different levels (control, ±10% change in deactivation rate of AMPK (k_dac AMPK). Green curve shows the control case, blue curve for +10% and magenta Curve for −10% of k_dac AMPK. (II) Sensitivity analysis of the width of the bistable region (length of the segment of the x axis between the black dotted lines in bifurcation diagram) for a parameter set (row #1 in Table S3) with changes in individual parameter by ±10% from the original value. Blue dotted line represents the ±10% change. Arrows show possible transitions between the different states.
Figure 3. Nullcline, bifurcation and sensitivity analysis. (A) Nullclines for three representative bistable parameter sets (rows 1–3 in Table S3). Green curve is AMPK Nullcline (d/dt (AMPK)=0). Red curve is Akt Nullcline (d/dt (AKT)=0). Blue circles represent stable states, white circles represent unstable steady state. (B) Bifurcation of Akt steady state levels with respect to the activation rate of AMPK (k_ac AMPK) for three representative bistable parameter sets (rows 1–3 in Table S3). Green curves denote stable states and red curves denote unstable states. Region bound by the blue dashed line represents the bistable region, where both states can co-exist. (C) (I) Representative bifurcation of Akt levels with respect to activation rate of AMPK (k_ac AMPK), drawn for three different levels (control, ±10% change in deactivation rate of AMPK (k_dac AMPK). Green curve shows the control case, blue curve for +10% and magenta Curve for −10% of k_dac AMPK. (II) Sensitivity analysis of the width of the bistable region (length of the segment of the x axis between the black dotted lines in bifurcation diagram) for a parameter set (row #1 in Table S3) with changes in individual parameter by ±10% from the original value. Blue dotted line represents the ±10% change. Arrows show possible transitions between the different states.
Jcm 10 00472 g003
Figure 4. Phase plots and stochastic simulations. (A) (I) Phase diagram for two parameters –activation rates of AMPK and Akt – showing monostable and bistable regions. All other parameter values correspond to those given in row #1 in Table S3. (II) Bifurcation diagram of AMPK levels with respect to k_ac AMPK for a constant value of k_ac Akt = 0.15. (III) Same as (II) but with respect to k_ac Akt for a constant value of k_ac AMPK=0.15. Green curve shows stable states, red curve shows unstable states. Blue dotted lines show region of bistability. (B) Stochastic simulations showing trajectories of AMPK, Akt values under the influence of noise for three representative parameter sets (rows 1–3 in Table S3). Noise parameter value η = 20.
Figure 4. Phase plots and stochastic simulations. (A) (I) Phase diagram for two parameters –activation rates of AMPK and Akt – showing monostable and bistable regions. All other parameter values correspond to those given in row #1 in Table S3. (II) Bifurcation diagram of AMPK levels with respect to k_ac AMPK for a constant value of k_ac Akt = 0.15. (III) Same as (II) but with respect to k_ac Akt for a constant value of k_ac AMPK=0.15. Green curve shows stable states, red curve shows unstable states. Blue dotted lines show region of bistability. (B) Stochastic simulations showing trajectories of AMPK, Akt values under the influence of noise for three representative parameter sets (rows 1–3 in Table S3). Noise parameter value η = 20.
Jcm 10 00472 g004
Figure 5. Stochastic state transitions. (A) Experimental validation of AMPK-Akt feedback loop using MDA-MB-231 EGR1-Turbo RFP cell lines sorted for high and low RFP expressing population using FACS. The RFP high population (red) corresponds to low pAMPK and high pAkt, and RFPlow (green) population corresponds to low EGR activity, high pAMPK, and low pAkt. Grey population correspond to cells that have lost the vector. Histograms show the population composition after 1, 7 and 9 days when started with distinct RFPhigh (red, top panel) and RFPlow (green, bottom panel) populations. (B) State transition graph: each node is a cell phenotype colored the same as FACS data for representative purposes, and each edge represents transition between the corresponding phenotypes. Transition rates (per day) calculated from the Markov model are shown on the corresponding arrows. (C) Predicted evolution of population heterogeneity when started from the initial population fractions, equivalent to distinct populations in the FACS data, using the transition probabilities calculated here. Error bars represent mean ± standard deviation from 1000 such simulations.
Figure 5. Stochastic state transitions. (A) Experimental validation of AMPK-Akt feedback loop using MDA-MB-231 EGR1-Turbo RFP cell lines sorted for high and low RFP expressing population using FACS. The RFP high population (red) corresponds to low pAMPK and high pAkt, and RFPlow (green) population corresponds to low EGR activity, high pAMPK, and low pAkt. Grey population correspond to cells that have lost the vector. Histograms show the population composition after 1, 7 and 9 days when started with distinct RFPhigh (red, top panel) and RFPlow (green, bottom panel) populations. (B) State transition graph: each node is a cell phenotype colored the same as FACS data for representative purposes, and each edge represents transition between the corresponding phenotypes. Transition rates (per day) calculated from the Markov model are shown on the corresponding arrows. (C) Predicted evolution of population heterogeneity when started from the initial population fractions, equivalent to distinct populations in the FACS data, using the transition probabilities calculated here. Error bars represent mean ± standard deviation from 1000 such simulations.
Jcm 10 00472 g005
Figure 6. Clinical validation of AMPK-Akt double negative feedback loop. Scatter plot of active levels of AMPK (AMPK pT172) and Akt (pT308 and pS473) in (A) Breast cancer cohort (TCGA-BRCA) (n = 901), (B) Pan Cancer cohort of 32 cancer types (n = 7694) and (C) Sarcoma cohort (TCGA-SARC) (n = 221). Each dot represents one patient and coordinates correspond to the protein expression levels capture using RPPA (reverse phase protein array) with Akt and AMPK active forms (AMPK-pT172, Akt-pT308 and Akt-pS473). Pearson’s correlation coefficients and p-value are reported.
Figure 6. Clinical validation of AMPK-Akt double negative feedback loop. Scatter plot of active levels of AMPK (AMPK pT172) and Akt (pT308 and pS473) in (A) Breast cancer cohort (TCGA-BRCA) (n = 901), (B) Pan Cancer cohort of 32 cancer types (n = 7694) and (C) Sarcoma cohort (TCGA-SARC) (n = 221). Each dot represents one patient and coordinates correspond to the protein expression levels capture using RPPA (reverse phase protein array) with Akt and AMPK active forms (AMPK-pT172, Akt-pT308 and Akt-pS473). Pearson’s correlation coefficients and p-value are reported.
Jcm 10 00472 g006
Table 1. Parameter values, their description and ranges used for random circuit perturbation simulation.
Table 1. Parameter values, their description and ranges used for random circuit perturbation simulation.
ParameterDescriptionValue Range
totalAMPKTotal level of AMPK100
totalAKTTotal level of AKT100
totalPHLPP2Total level of PHLPP2100
totalPP2CαTotal level of PP2Cα100
kac_AMPKActivation rate of AMPK(0.02–0.2)
kac_AKTActivation rate of AKT(0.02–0.2)
kac_PHLPP2Activation rate of PHLPP2(0.02–0.2)
kac_PP2CαActivation rate of PP2Cα(0.02–0.2)
kdac_AMPKDeactivation rate of AMPK(0.02–0.2)
kdac_AKTDeactivation rate of AKT(0.02–0.2)
kdac_PHLPP2Deactivation rate of PHLPP2(0.02–0.2)
kdac_PP2CαDeactivation rate of PP2Cα(0.02–0.2)
λPP2cαEffect of PP2Cα on AMPK(5–10)
λPHLPP2Effect of PHLPP2 on AKT(5–10)
λAKTEffect of AKT on PP2Cα(5–10)
λAMPKEffect of AMPK on PHLPP2(5–10)
nPP2CαHill coefficient of PP2Cα for deactivation of AMPK4, 5, 6
nPHLPP2Hill coefficient of PHLPP2 for deactivation of AKT4, 5, 6
nAKTHill coefficient of AKT for activation of PP2Cα4, 5, 6
nAMPKHill coefficient of AMPK for activation of PHLPP24, 5, 6
PP2Cα0Threshold value of PP2Cα for deactivation of AMPK(0.25–0.75) × totalPP2Cα
PHLPP20Threshold value of PHKLPP2 for deactivation of AKT(0.25–0.75) × totalPHLPP2
AMPK0Threshold value of AMPK for activation of PHLPP2(0.25–0.75) × totalAMPK
AKT0Threshold value of AKT for activation of PP2Cα(0.25–0.75) × totalAKT
AMPK: AMP-activated Protein Kinase, PHLPP2: Plectistin Homology Domain And Leucine Rich Repeat Protein Phosphatase 2, PP2Cα: Protein Phosphatase 2C alpha.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Chedere, A.; Hari, K.; Kumar, S.; Rangarajan, A.; Jolly, M.K. Multi-Stability and Consequent Phenotypic Plasticity in AMPK-Akt Double Negative Feedback Loop in Cancer Cells. J. Clin. Med. 2021, 10, 472. https://doi.org/10.3390/jcm10030472

AMA Style

Chedere A, Hari K, Kumar S, Rangarajan A, Jolly MK. Multi-Stability and Consequent Phenotypic Plasticity in AMPK-Akt Double Negative Feedback Loop in Cancer Cells. Journal of Clinical Medicine. 2021; 10(3):472. https://doi.org/10.3390/jcm10030472

Chicago/Turabian Style

Chedere, Adithya, Kishore Hari, Saurav Kumar, Annapoorni Rangarajan, and Mohit Kumar Jolly. 2021. "Multi-Stability and Consequent Phenotypic Plasticity in AMPK-Akt Double Negative Feedback Loop in Cancer Cells" Journal of Clinical Medicine 10, no. 3: 472. https://doi.org/10.3390/jcm10030472

APA Style

Chedere, A., Hari, K., Kumar, S., Rangarajan, A., & Jolly, M. K. (2021). Multi-Stability and Consequent Phenotypic Plasticity in AMPK-Akt Double Negative Feedback Loop in Cancer Cells. Journal of Clinical Medicine, 10(3), 472. https://doi.org/10.3390/jcm10030472

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop