 Research article
 Open Access
 Published:
In silico modeling of phosphorylation dependent and independent cMyc degradation
BMC Bioinformatics volume 20, Article number: 230 (2019)
Abstract
Background
cMyc plays an important role in cell proliferation, cell growth and in differentiation, making it a key regulator for carcinogenesis and pluripotency. Tight control of cmyc turnover is required by ubiquitinmediated degradation. This is achieved in the system by two Fbox proteins Skp2 and FBXW7.
Results
Dynamic modelling technique was used to build two exclusive models for phosphorylation dependent degradation of Myc by FBXW7 (Model 1) and phosphorylation independent degradation by Skp2 (Model 2). Sensitivity analysis performed on these two models revealed that these models were corroborating experimental studies. It was also seen that Model 1 was more robust and perhaps more efficient in degrading cMyc. These results questioned the existence of the two models in the system and to answer the question a combined model was hypothesised which had a decision making switch. The combined model had both Skp2 and FBXW7 mediated degradation where again the latter played a more important role. This model was able to achieve the lowest levels of ubiquitylated Myc and therefore functioned most efficiently in degradation of Myc.
Conclusion
In this report, cMyc degradation by two Fbox proteins was mathematically evaluated based on the importance of cMyc turnover. The study was performed in a homeostatic system and therefore, prompts the exploration of cMyc degradation in cancer state and in pluripotent state.
Background
cMyc protein is a shortlived transcription factor that plays an important role in cell proliferation, apoptosis, cell growth, angiogenesis and pluripotency [1]. The structure of cMyc protein is divided into N terminal transcription activation domain (TAD) and Cterminal basichelixloophelixleucinezipper (bHLHLZip) domains. The bHLHLZip region is responsible for dimerization with Max and binding to Eboxes of target gene promoters. Whereas Max is expressed constitutively, Myc is transient. The halflife of cmyc is short (~30mins) in proliferating cells [2] and needs to be tightly regulated as overexpression leads to tumour formation. It is known that cMyc undergoes ubiquitination followed by degradation in the proteasome [3] and the regions responsible for this signalling is believed to be Myc Box 1 (MBI) and Myc Box 2 (MBII) in the TAD [4]. In fact it is seen that these regions, especially MBI are the hotspots for mutations present in cancer [5].
Ubiquitin mediated degradation in the proteasomes is enabled by three enzymes; E1 for ubiquitin activation, E2 for ubiquitin conjugation and E3 for ubiquitin ligation which confers substrate specificity [6]. Two types of E3 ubiquitin ligases are there in the system, namely SCF complex and anaphasepromoting complex or cyclosome (APC /C). The SCF complex itself has four components: SKP1, Cul1, Rbx1 and a variable Fbox protein which determines target specificity [7]. Among the many Fbox proteins identified, Skp2 and FBXW7 have been well characterized in their role of degradation of p27^{Kip2} [8] and Cyclin E [9] respectively. Evidence for binding of both Skp2 and FBXW7 to cMyc protein is present in literature [10]. Even though both are Fbox proteins responsible for degradation, a few key differences may be noted between them. Skp2 contains leucinerich repeats along with an FBox domain and binds to cMyc via MBII and HLHLZip domains. It degrades cMyc independent of phosphorylation [11]. FBXW7 on the other hand, contains WD40 repeats along with FBox and binds to cMyc via MBI. It requires phosphorylation at S62 and T58 positions for ubiquitination to take place [11] (Fig. 1). Not only do the 2 Fbox proteins have different structures and modes of interaction with cMyc, another important difference lies in the fact that Skp2 is an oncoprotein whereas FBXW7 is a tumour suppressor [12].
Although it is known that both FBXW7 and Skp2 regulate cMyc by degradation, information regarding the distinctions of the tightly regulated degradation process remains to be explored. The turnover of cMyc plays an important role in the cells decision to proliferate or differentiate in pluripotent cells [13]. Overexpression of Myc is also linked with cancer development [14]. Hence, it is of utmost importance to understand the degradation pathway of cmyc so that any aberrant overexpression can be therapeutically targeted. One of the ways in which control on Myc level is achieved is at the posttranslational level via protein stability modulation [15]. Therefore, in this study, role of posttranslational modifications in cMyc function is evaluated giving prime focus on phosphorylation and ubiquitination.
It is already known that sequential phosphorylation of cMyc takes place by Erk and GSK3β, which is followed by ubiquitination by FBXW7 and degradation in proteasome [16]. Degradation of Myc may also take place independent of the phosphorylation steps by Skp2 protein [17]. Therefore, for our study we have selected these two pathways that degrade cMyc protein. In addition, questions remain as to how the system chooses which pathway to take for the turnover of Myc. Since experimental studies are time consuming and perturbation of cellular systems require resources, we have aimed to do a preliminary insilico study to establish the critical steps required for cMyc degradation. Despite many studies with cMyc, which resulted in thousands of publication in the past three decades, how Myc expression is still able to get past this regulation and cause cancer is a question that remains unanswered. The results would give us a better insight into the biology of ubiquitination and degradation of cMyc.
Mathematical modelling has been valuable in assessing the viability of potential therapeutic strategies by identifying critical steps required for regulation. Therefore, to understand the relation between the two types of regulation of cMyc we have constructed two dynamical models responsible for cMyc degradation. In the first model, which is a complete and improved version of a model made by Lee et al., we have considered the cMyc phosphorylation at two residues and then FBXW7 dependent ubiquitination followed by degradation [18]. In the second model, we have dealt with Skp2 mediated ubiquitination of cMyc independent of phosphorylation. The synchronized consequences of the three signals make the system control the diverse cell fates [19]. We have performed a parameter sensitivity analysis on the two models within the framework of various correlation coefficients to identify the contribution of the modular structures in signal propagation for both the models independently. In the view of signalling and cMyc degradation, the Fbxw7 mediated pathway shows more robustness over the Skp2 pathway for a large range of alteration in the system components as well as signal components. Questions may arise though as to how Model 1, which is energetically more expensive than Model 2, is the more robust mechanism for cMyc degradation. To address this query, we have constructed a full model by combining the two independent models in order to see how they perform together. Since it did not fit well with the experimental observations and was not able to efficiently degrade Myc, we incorporated an exclusive ‘onoff switch’ between Model 1 and 2 and tried to find out how and under what circumstances will the system flip from one model to the other.
Results
Model building
In the past few decades, bioinformatics tools have been extensively used in protein degradation related studies. They have also been instrumental in construction of many large scale databases. Nevertheless, these largescale efforts are mostly restricted in providing only a static picture of a network. We have hypothesized two separate dynamic models (Fig. 2) and explored the possibility of their coexistence for cMyc degradation taking into consideration normal homeostasis state. This may allow us to predict the changes that may occur in disease state In this study, we have used dynamical model based approach, which is fast evolving into the most promising tool for such purpose.
In Model 1, the levels of the signals ERK and GSK were evaluated using the model from Lee et al. and deriving the GSK signal from it [18] (Fig. 3a). As is reported, Myc is highly unstable when synthesized and is stabilized when it undergoes phosphorylation at Serine 62 by Erk [20]. Phosphorylation at Threonine 58 by Gsk3β starts to destabilize Myc and primes it for ubiquitination [21]. It is also seen that in some cell lines Erk has an early transient pulse when stimulated by growth factors (GFs) [22]. On the other hand, PI3K, an antagonist of Gsk3β, has two peaks in some cell lines when stimulated by fetal bovine serum (FBS) [23]. This leads us to derive a signal for Gsk3β having the kind of curve depicted in Fig. 3a [18] (See Methods). It has been experimentally shown that Myc phosphorylations are essential for its binding to Fbxw7 and that Gsk3β inhibitor reduced the interaction between Myc and Fbxw7 [11]. Therefore, we can say that the activation feedback on Fbxw7 comes from x_{3}, which is Myc phosphorylated at T58 by GSK, and that Fbxw7 triggers the ubiquitination. The existence of this delayed activation is clearly depicted in the Fig. 3a. Fbxw7 pulse is also generated through assumptions from previous models and experimental results [11] in form of an ODE (Eq) described in the methods section. This part of the Model 1 has been added by our group to the derived version from Lee et al. and extended to phosphorylation dependent ubiquitination models from Nguyen et al. [24]. In this consequence, the onoff timing of the Erk and GSK3β is very crucial. The activation of Myc also integrates and correlates with the input signals. Until Erk stays on (~ 2 h) the pool of cMyc (x_{1–2}) starts growing in a steep fashion, though x_{3} and x_{4} do appear after the Erk signal turns off (Fig. 3b). The most complicated pulsed signal of GSK3β helps the system by introducing delay in the uiquitination process. In presence of these variations in signals, the time profiles of four states of Myc proteins (x_{1}–_{4}) are generated by solving the set of ODEs (Eq.in Methods) (Fig. 3b).
In the Model 2, the levels of Skp2 signal is regulated in feedback by cMyc, and encourages the direct ubiquitination process. Experimental evidence suggests that overexpression of Skp2 in Rat1 cells caused enhanced degradation of Myc although in control cells Myc was stabilized [10]. We have derived Model 2 from a generic phosphorylation independent ubiquitination model explained by Nguyen et al [24]. and used experimental data to make approximations of the parameters involved. The kinetics of the two states of cMyc and Skp2 are evaluated by solving the corresponding dynamical equations (Eq.) and represented in (Fig. 3c). The parameter sensitivity analysis along with the robustness estimation we have performed on the two models gives an idea about the comparative importance of the defined model parameters.
Sensitivity analysis
For sensitivity analysis, a dimensionless quantity γ (ratio of ubiquitylated Myc with total Myc in all four populations) was calculated for Model 1 and γ’ (ratio of ubiquitylated Myc with total Myc in two populations) for Model 2 (see Methods section). Three types of correlation values were calculated to get an idea of which parameter is most sensitive (Tables 1 and 2). Additionally, scatter plots for the correlation values clearly show the positive and negative correlations (see Methods section for rate constants and other parameters) (Fig. 4a and b). Here, we have computed the correlation indices between the individual model input parameters and the model output. Although we have sampled the model parameters from a normal distribution, it does not make sure that the model output is also normally distributed. In general, for such network of interacting variables, if there exists any nonlinear relation in the model equation, the output variable does not mimic the same distribution as the inputs. If they do, as our results show, then the CC and RCC values will show similar values. In the present study, we are able to check the presence of any nonlinearity in the model by calculating both the CC and RCC. PRCC, on the other hand, can magnify the onetoone correlation between two variables after removing all the controlling or confounding effects of other model variables. For the case of low strength confounders, the cofactors of each elements would show similar values as the corresponding elements leading to PRCC values same as the RCC values as seen in our results.
In Model 1, results show that rate constant of GSK3β mediated phosphorylation of x_{2} (k_{5}) and GSK3β levels have high positive correlation with γ whereas degradation at Myc T58 phosphorylated x_{3} state (k_{6}) as well as degradation of ubiquitylated state x_{4} (k_{11}) have negative correlation (Table 1). Experimental evidence from Yada et al. shows that using GSK3 inhibitor reduced the association of Fbxw7 with cMyc and delayed Myc turnover [11]. This explains the high positive correlation we have found in our model. The correlation coefficients when intensity of perturbation is changed in the range of 5–25% also observes the same trend (data not shown). In Model 2, we see positive correlation of Myc interacting with Skp2 (k_{12}) and Skp2 ubiquitylating Myc (k_{13}) with γ’, whereas degradation of Skp2 (k_{14}) shows a negative correlation (Table 2). This result is reflected in experimental studies by von del Lehr et al. as they show that Skp2 is interacting with Myc in a positive feedback loop where Skp2 deletion mutation led to accumulation of cMyc [17].
We observe that k_{1} and k_{2} have higher correlation values with γ’ in Model 2 as compared to correlation values with γ in Model 1. It is however noteworthy that degradation of ubiquitylated Myc (k_{11}) has higher correlation values (negative correlation) in Model 1 as compared to Model 2. This could indicate that Model 1 is more efficient in degrading Myc. This is also reflected in studies by von der Lehr group which indicates that Skp2 binds with Myc during S phase entry and promotes transcription targets of Myc and Skp2. As they hypothesised that this makes Skp2 mediated degradation a necessary step for activation of transcription [17]. Additionally, Yada et al. have also suggested that inhibition of degradation would be seen only in the initial phase if Skp^{−}/^{−} mutation is present in MEFs, as Fbxw7 or other degradation factors will take over the responsibility of degrading Myc later [11]. These evidences may seem to establish Fbxw7 mediated degradation pathway as the more efficient degradation system.
Although ttest was performed for all the correlation coefficients but since the system is linear, performing this analysis did not add any significance to the data (Additional file 1: Table S3 and S4). From these results, we can say that the two models have parity with invivo and invitro conditions of Myc degradation and correspond to results generated by experimental studies [11, 17].
Robustness
In the parameter sensitivity analysis for the two models, we have explored the weight of contribution of each of the model parameters on the model outputs. To check the efficiency of stability of the models, we calculate the robustness in terms of model output with respect to the intensity of perturbation. With the rise of the perturbation intensity on the model parameters, the models generally to lose their stability and deviate far from the stable steady state. From the spread of the scatter plots and the deviations in the moving averages of the model outputs, we can compare the two models quantitatively (Fig. 5a and b). The moving average value calculated for Model 1 stays approximately invariant in the range of the total parameter variation while the same for Model 2 start fluctuating at higher parameter variations. These comparative results simply can infer that Model 1 is more robust than Model 2.
It is known that Skp2 is expressed in Sphase while Fbxw7 is expressed throughout the cell cycle [10]. It is also seen that, overexpression of Fbxw7 reduced transactivation by Myc as compared to overexpression of Skp2, which contributes to growth promoting effects [12]. Finally, we know that mutations in Thr58 and Ser62 contributes to cancer [25] therefore giving an importance to these residues for degradation step. From these evidences listed above, we can explain why Model 1 should in fact be more robust than Model 2.
Decisionmaking
In a combined model, if we incorporated every aspect from the two models in logical succession and evaluated the time profiles for the different states of cMyc it was seen that low concentration of x_{4} could not be reached (Additional file 1: Figure S3). Therefore, based on the information reported till date, it may be concluded that the two models perform exclusively perhaps in different phases of the cell cycle [26, 27]. Infact, as already mentioned before, Skp2 expression is specific to Sphase of the cell cycle. The cell cycle stage or other spatiotemporal factors may therefore lead to the degradation flipping from one pathway to the other. To mimic their exclusive contribution in homeostasis, i.e., the degradation of ubiquitinated cMyc, we have designed an onoff parameter α that switches between the activation of Fbxw7 and Skp2 (See methods). As model 1 is more dominant in the sense of robustness, we consider it as the on state of α and the model 2 as the off state of α. On sampling the switching function α with different durations and onoff times iteratively, we found two states of α, one ‘on’ states and one ‘off’ state. α stays off up to 3 h from the initial time point and gets on until 30 h, the time course of the observation (Fig. 6a). As an effect of this exclusive activation of the two models, the ubiquitinated cMyc population is attenuated (Fig. 6b). Our prediction is in line with experimental evidence which suggests that in Fbw7^{−/−} ES cells degradation of Myc was impaired for up to 60 min whereas in Skp2^{−/−} MEFs Myc degradation was impaired for 20 min [11]. As described before, in homeostasis, it can be concluded that FBXW7 mediated degradation plays the dominant role in degradation of Myc.
Discussion
The oncoprotein cMyc, a basic helixloophelix/leucine zipper (bHLH/Zip)  type transcription factor, is a master regulator of cell proliferation [28]. The expression of this protein is transient and is responsible for cell proliferation. Its amplification or mutation is present in most types of cancers and therefore, its turnover is a critical determinant of carcinogenesis [29]. Myc is ubiquitinated and marked for degradation in the proteasome [15]. Although we know that regions of Myc like MBI and MBII are not responsible for degradation directly, but they play an important role in binding of E3 ligases [4]. Two such ligases that have been indicated in this paper are FBXW7 and Skp2. Based on experimental evidences and models made by other groups we have discussed two separate models for cMyc degradation in this paper. Model 1 is a phosphorylation dependent model via FBXW7 whereas Model 2 is a phosphorylation independent model via Skp2.
Dynamic modelling method was used for designing the two models and rate parameters were estimated from experiments as well as models designed by other groups. In this study, for the sake of reducing complexity, we have ignored the reverse reactions in all cases. Other assumptions made include the fact that we have considered the system to be at homeostasis and how the model will function in the diseased state has not been explored. Other limitations include the fact that Skp2 and FBXW7 mediated degradations are not the only two pathways that are responsible for degrading Myc. Other pathways have been ignored for making the model simple.
From the results obtained in our models and comparing them with experimental data cited we can conclude that the mechanisms of degradation hypothesised by us is in line with experimentally proved in the biological system. The results also indicate that Model 1 is more robust and degradation rate constant k_{11} has a higher correlation with output signal γ. This happens because efficient degradation of cMyc requires a protein, which is a tumour suppressor, like FBXW7. Skp2 on the other hand is an oncoprotein, which along with ubiquitylating Myc also helps in its transactivation. The role of Skp2 remains unclear, but it is evident that Model 1 is a more reliable method of regulating cMyc. This is also evident as mutations in S62 and T58, residues involved in phosphorylation according to model 1, are responsible for cancer phenotype [25].
To find out why two different mechanisms exist in vivo and to understand its role on degradation of Myc, we have combined the two models. Results indicate that when these two models work together simultaneously, they are not efficient to reduce the ubiquitylated Myc population (Additional file 1: Figure S3). This may imply that even in the biological system, the two methods for degradation work exclusively. Therefore, a switch was incorporated to decide which E3 ligase would work at which point of time. In the cell, it is perhaps the cell cycle state, which decides when the two E3 ligases operate. It can be concluded that by combining the two models in a decision dependent manner lower levels of ubiquitylated Myc could be achieved. It should also be noted that in the combined model, FBXW7 plays a more predominant role than Skp2 in cMyc degradation.
We can find some clues from this combined model as to which steps are critical in causing Myc overexpression leading to cancer. However, the entire dynamic modelling was done for the system in homeostasis. Myc is known to play an important role in cancer state as well as in pluripotency. It is an important player, which modulates expression of multiple other genes to switch from normal state to cancer or from nondividing state to proliferative state [30]. This area remains to be explored further and the use of modelbased methods to decipher therapeutic strategies remains to be the target for researchers. It is also expected that modelling and modelbased approaches in integration with experimentation will become a major tool for data interpretation and hypothesis generation in all fields of biology.
Conclusion
cMyc is a transcription factor responsible for cell fate decisions. The presence of many Fbox proteins functioning as E3 ubiquitin ligase for degradation of cMyc protein has been reported. In this study two such proteins Skp2 and FBXW7 have been explored for their role in degradation. It has been shown that phosphorylation dependent degradation via FBXW7 is a more robust mechanism for degradation in spite of which phosphorylation independent degradation via Skp2 also takes place. Therefore a putative mechanism of degradation that takes place in the system has been hypothesised by combining the two models in a decision dependent manner. This gives way to understanding how cMyc may be regulated in the cell and question how in diseasestate as well as during pluripotency this process is altered.
Methods
Phosphorylation and ubiquitination often have opposing effects on target proteins and require the interaction of different partners. This creates complications and makes the model quite large, which renders model construction and parameter estimation quite challenging. In this study, not only have we combined phosphorylation with ubiquitination, but we have also given two alternative pathways by which Myc can be ubiquitinated. Finally, we have combined the two models to give a better idea as to what happens in the system.
Phosphorylation dependent cMyc degradation: model 1
The role of posttranslational modifications (PTMs) in cMyc was investigated in a kinetic modeling framework. Myc is stimulated by external growth factors (GF) and further gets a sequential phosphorylation by Erk (E) and GSK3β (G) at Serine 62 (S62) and Threonine 58 (T58) respectively. The T58 phosphorylation occurs with the dephosphorylation at S62. The T58 phosphorylation leads to activation of the Fbox protein FBXW7, which eventually triggers the ubiquitination of Myc and hence degradation (Fig. 2). To make the model simplified and tractable, we have not considered the phosphatases and the DUBs in the present network. The four states of cMyc are depicted by x_{1}, x_{2}, x_{3}, and x_{4} for GF stimulated cMyc, cMyc phosphorylated at S62, cMyc phosphorylated at T58 and ubiquitinated, respectively. The kinetic reactions of the cascade have been depicted in Additional file 1: Table S1 along with reaction rates. The rate constants were derived from previous models and were subjected to estimation and assumptions [18, 24, 31, 32].
The Erk and GSK3β signals were introduced in pulses of different structures, while the GF induces the cMyc activation constantly. The Pulse amplitudes were kept at 1.0 but the onoff switch durations for Erk and GSK3β are different. Erk was on up to 1.0 h and got turned off as soon as GSK3β signal turned on. The nature of the GSK3β signal was complicated; it turned on at 1.0 h, was turned off at 2.0 h and again turned on at 7.0 h and stayed on until the end of the time course. All rate parameters can be found in the Additional file 1: Table S1.The previously mentioned degradation cascade can be presented in a set of coupled ordinary differential equations as
Where, the pulse of Erk and GSK3β are constructed by the combination of several Heaviside step functions as
The system of ODEs has been solved for 30 h, as described in a previous model [18]. A dimensionless fraction of ubiquitinated cMyc (x_{4}/(x_{1} + x_{2} + x_{3} + x_{4}) = γ) has been considered as the model output. At steady state, we performed an input parameter sensitivity analysis and checked the robustness of the model with respect to this model output. The methodology of the correlation coefficient based sensitivity analysis and robustness calculation will be described later.
Concentration of growth factor: GF = 1.0;
Details of the Erk pulse: E_{Max} = 0.9; E_{R} = 0.1; Dur_{E} = 1.0;
Details of the Gsk pulse: G_{Max} = 0.9; G_{R} = 0.1; Dur_{G1} = 2; Dur_{G2} = 7; Width_{G} = 0.5;
Details of the Fbxw7 pulse: F_{T} = 5.0.
Phosphorylation independent cMyc degradation: model 2
Apart from the Fbxw7 dependent pathway of cMyc degradation, a Skp2 mediated pathway exists and it is independent of phosphorylation. The Skp2 pathway of cMyc degradation is triggered by the activation of Skp2 by cMyc. Skp2 starts ubiquitination of cMyc for degradation in the proteasome. We have framed the model in a set of ordinary differential equations as
where, x_{1}, x_{4} are the states of cMyc before and after ubiquitination and S stands for Skp2. The detail of the reactions is tabulated in Additional file 1: Table S2. Similarly, a dimensionless fraction of ubiquitinated cMyc (x_{4}/(x_{1} + x_{4}) = γ’) has been considered as the model output. Some parameters were taken from the previous model.
Parameter sensitivity
To quantify the contribution of each model parameter simultaneously, we prefer to utilize the correlation coefficient based parameter sensitivity analysis. We calculate three kinds of correlations, the Pearson’s (CC), Spearman’s rank (RCC), and the partial rank (PRCC) correlation coefficients as reported in a previous article [33]. All the correlation coefficients were measured at steady state only with increasing the strength of perturbation, but it does not show any notable change until 20%. Within the period of 30 h, the Erk signal goes ‘off’ from ‘on’ state once and the Gsk signal flips thrice. We try to find the variation in the parameter sensitivity profiles near these flip points for Model 1, but we cannot observe any difference. p –values was calculated for both models using ttest and results were tabulated in Additional file 1: Tables S3 and S4.
Robustness
The efficiency of a biochemical network to maintain its functionality in the very altering environment can be quantitated. Measure of the model robustness is one of the simple and reliable approaches. In the current paper, robustness of the two models are analysed in terms of ensembles of perturbed systems. The rate parameters and system inputs like the duration and width of the Erk and Gsk signals were sampled randomly from Gaussian distribution about their native values and corresponding spread of the network output (γ and γ’) is measured. The robustness was characterized in terms of the summation of normalized alteration in the parameter set for an ensemble of 5000 models [34]. We increase the perturbation intensity and note the ensemble of the variation in the γ and γ ‘value from the original.
Decision process
The previously mentioned two models of cMyc degradation exist in the system but perform exclusively. We predict there must be some hidden agents, who decide the switch between the two models. To monitor the decision making process, we assemble the two models in to a single one and consider the activation of Fbxw7 and Skp2 as exclusive events.
where α is the switching parameter, flips between 0 and 1. When Fbxw7 is activated, the model 1 will be on board whether activation of Skp2 triggers the model 2 in action. To explore the switching time and preference of the two models we have defined a switching parameter α that triggers activation of the models exclusively at preferred time points. The combined model was framed as an optimization problem with an objective to reduce the population of x_{4} in the system over the time course. To perform the optimization, we have constructed an objective function as
to be minimized. To optimize the temporal evolution of α throughout the course of the full signaling network, we need to construct a generalized step switch with flexible ‘onoff’ at any time. A generalized Pi function has been considered with duration and ‘onoff’ time sampled from a uniform distribution over the full timescale. On running 10,000 independent iterations, we have found the optimal structure of the switching function. For the optimization runs and series of sampled structures of the switch please find the Additional file 1: Figure S1S2.
Abbreviations
 APC /C:

Anaphasepromoting complex or cyclosome
 bHLHLZip:

Basichelixloophelixleucinezipper
 CC:

Pearson’s correlation coefficient
 Erk:

Extracellular signal–regulated kinases
 ES cells:

Embryonic stem cells
 FBS:

Fetal bovine serum
 FBXW7:

FBox And WD Repeat Domain Containing 7
 GF:

Growth factors
 GSK3β:

Glycogen synthase kinase3 beta
 MBI:

Myc Box 1
 MBII:

Myc Box 2
 MEFs:

Mouse embryonic fibroblasts
 ODEs:

Ordinary differential equations
 PI3K:

Phosphatidylinositol3kinases
 PRCC:

Partial rank correlation coefficient
 Rbx1:

RINGbox protein 1
 RCC:

Spearman’s rank correlation coefficient
 S62:

Serine at position 62
 SCF:

Skp, Cullin, Fbox containing complex
 Skp1:

Sphase kinaseassociated protein 1
 Skp2:

Sphase kinaseassociated protein 2
 T58:

Threonine at position 58
 TAD:

Transcription activation domain
References
 1.
Dang CV, Lewis BC. Role of oncogenic transcription factor cMyc in cell cycle regulation, apoptosis and metabolism. J Biomed Sci. 1997;4(6):269–78.
 2.
Hann SR, Eisenman RN. Proteins encoded by the human cmyc oncogene: differential expression in neoplastic cells. Mol Cell Biol. 1984;4(11):2486–97.
 3.
Ciechanover A, DiGiuseppe JA, Schwartz AL, Brodeur GM. Degradation of MYCN oncoprotein by the ubiquitin system. Prog Clin Biol Res. 1991;366:37–43.
 4.
Grandori C, Cowley SM, James LP, Eisenman RN. The Myc/max/mad network and the transcriptional control of cell behavior. Annu Rev Cell Dev Biol. 2000;16:653–99.
 5.
Bahram F, von der Lehr N, Cetinkaya C, Larsson LG. CMyc hot spot mutations in lymphomas result in inefficient ubiquitination and decreased proteasomemediated turnover. Blood. 2000;95(6):2104–10.
 6.
Ciechanover A. Proteolysis: from the lysosome to ubiquitin and the proteasome. Nat Rev Mol Cell Biol. 2005;6(1):79–87.
 7.
Elledge SJ, Harper JW. The role of protein stability in the cell cycle and cancer. Biochim Biophys Acta. 1998;1377(2):M61–70.
 8.
Carrano AC, Eytan E, Hershko A, Pagano M. SKP2 is required for ubiquitinmediated degradation of the CDK inhibitor p27. Nat Cell Biol. 1999;1(4):193–9.
 9.
Koepp DM, Schaefer LK, Ye X, Keyomarsi K, Chu C, Harper JW, et al. Phosphorylationdependent ubiquitination of cyclin E by the SCFFbw7 ubiquitin ligase. Science. 2001;294(5540):173–7.
 10.
Kim SY, Herbst A, Tworkowski KA, Salghetti SE, Tansey WP. Skp2 regulates Myc protein stability and activity. Mol Cell. 2003;11(5):1177–88.
 11.
Yada M, Hatakeyama S, Kamura T, Nishiyama M, Tsunematsu R, Imaki H, et al. Phosphorylationdependent degradation of cMyc is mediated by the Fbox protein Fbw7. EMBO J. 2004;23(10):2116–25.
 12.
Bashir T, Pagano M. Aberrant ubiquitinmediated proteolysis of cell cycle regulatory proteins and oncogenesis. Adv Cancer Res. 2003;88:101–44.
 13.
Varlakhanova NV, Cotterman RF, deVries WN, Morgan J, Donahue LR, Murray S, et al. Myc maintains embryonic stem cell pluripotency and selfrenewal. Differentiation. 2010;80(1):9–19.
 14.
Dang CV. MYC on the path to Cancer. Cell. 2012;149(1):22–35.
 15.
Farrell AS, Sears RC. MYC Degradation. Cold Spring Harb Perspect Med. 2014;4(3):a014365.
 16.
Lutterbach B, Hann SR. Hierarchical phosphorylation at Nterminal transformationsensitive sites in CMyc protein is regulated by mitogens and in mitosis. Mol Cell Biol. 1994;14(8):5510–22.
 17.
von der Lehr N, Johansson S, Wu SQ, Bahram F, Castell A, Cetinkaya C, et al. The Fbox protein Skp2 participates in cMyc protelosornal degradation and acts as a cofactor for cMycregulated transcription. Mol Cell. 2003;11(5):1189–200.
 18.
Lee T, Yao G, Nevins J, You L. Sensing and integration of Erk and PI3K signals by Myc. PLoS Comput Biol. 2008;4(2):e1000013.
 19.
Novak B, Tyson JJ. A model for restriction point control of the mammalian cell cycle. J Theor Biol. 2004;230(4):563–79.
 20.
Gregory MA, Hann SR. CMyc proteolysis by the ubiquitinproteasome pathway: stabilization of cMyc in Burkitt's lymphoma cells. Mol Cell Biol. 2000;20(7):2423–35.
 21.
Yeh E, Cunningham M, Arnold H, Chasse D, Monteith T, Ivaldi G, et al. A signalling pathway controlling cMyc degradation that impacts oncogenic transformation of human cells. Nat Cell Biol. 2004;6(4):308–18.
 22.
Sasagawa S, Ozaki Y, Fujita K, Kuroda S. Prediction and validation of the distinct dynamics of transient and sustained ERK activation. Nat Cell Biol. 2005;7(4):365–73.
 23.
Kumar A, Marques M, Carrera AC. Phosphoinositide 3kinase activation in late G1 is required for cMyc stabilization and S phase entry. Mol Cell Biol. 2006;26(23):9116–25.
 24.
Nguyen LK, Kolch W, Kholodenko BN. When ubiquitination meets phosphorylation: a systems biology perspective of EGFR/MAPK signalling. Cell Commun Signal. 2013;11:52.
 25.
Chakravorty D, Jana T, Das Mandal S, Seth A, Bhattacharya A, Saha S. MYCbase: a database of functional sites and biochemical properties of Myc in both normal and cancer cells. BMC bioinformatics. 2017;18(1):224.
 26.
Hara T, Kamura T, Nakayama K, Oshikawa K, Hatakeyama S, Nakayama K. Degradation of p27(Kip1) at the G(0)G(1) transition mediated by a Skp2independent ubiquitination pathway. J Biol Chem. 2001;276(52):48937–43.
 27.
Spruck CH, Strohmaier H, Sangfelt O, Muller HM, Hubalek M, MullerHolzner E, et al. hCDC4 gene mutations in endometrial cancer. Cancer Res. 2002;62(16):4535–9.
 28.
Oster SK, Ho CS, Soucie EL, Penn LZ. The myc oncogene: MarvelouslY complex. Adv Cancer Res. 2002;84:81–154.
 29.
DallaFavera R, Bregni M, Erikson J, Patterson D, Gallo RC, Croce CM. Human cmyc onc gene is located on the region of chromosome 8 that is translocated in Burkitt lymphoma cells. Proc Natl Acad Sci U S A. 1982;79(24):7824–7.
 30.
Rangarajan N, Fox Z, Singh A, Kulkarni P, Rangarajan G. Disorder, oscillatory dynamics and state switching: the role of cMyc. J Theor Biol. 2015;386:105–14.
 31.
Nguyen LK. Dynamics of ubiquitinmediated signalling: insights from mathematical modelling and experimental studies. Brief Bioinform. 2016;17(3):479–93.
 32.
Xu L, Qu Z. Roles of protein ubiquitination and degradation kinetics in biological oscillations. PLoS One. 2012;7(4):e34616.
 33.
Mapder T, Talukder S, Chattopadhyay S, Banik SK. Deciphering parameter sensitivity in the BvgAS signal transduction. PLoS One. 2016;11(1):e0147281.
 34.
Barkai N, Leibler S. Robustness in simple biochemical networks. Nature. 1997;387(6636):913–7.
Acknowledgements
The authors thank Bioinformatics Centre of Bose Institute for providing the infrastructure for the studies carried out for this paper. We thank Prof. Indrani Bose an Emeritus Scientist, Department of Physics, Bose Institute for her valuable feedback. We would also like to thank ICMR for funding the project number BIC/12(30) /2012.
Funding
DC is supported by Bose Institute Senior Research Fellowship. Bose Institute is not involved in the design or conclusion of the study. KM is supported by ICMR project number BIC/12(30) /2012. ICMR did not have a role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Availability of data and materials
Parameters considered in the models and correlation coefficients derived during this study are included in this published article [and its Additional file 1]. Software specific codes are available as pdf files zipped and saved under Additional file 2.
Author information
Affiliations
Contributions
SS and TM conceived the experimental design. DC and TM performed the experimental studies. SS, DC, KB and TM analysed the data. SS, DC and TM wrote the manuscript. All authors have read and approved the manuscript.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional files
Additional file 1:
Figure S1. Few representative replicas from the sample set of the switching function gives better understanding of the different switching time and states. Figure S2. The profile of the objective function (J) over 100 samples. We have minimized J over 10,000 independent sample runs. Figure S3. Combined Model of Model 1 and Model 2. A: All input signals of the two models are combined and shown vs time; B: The concentration of cmyc populations with respect to time is shown. Colour codes are given for both graphs. Additionally we can observe that × 4 value is not as low as expected (Compare with Fig. 6). Table S1. Reaction scheme of the phosphorylation dependent degradation of cMyc, Model 1. Table S2. Reaction scheme of the phosphorylation independent degradation of cMyc, Model 2. Table S3. pvalues for the correlation coefficients in Model 1. Table S4. pvalues for the correlation coefficients in Model 2. (DOCX 560 kb) (DOCX 560 kb)
Additional file 2:
This file contains software specific codes for all the analysis performed in this study including the codes for the two combined models in a pfd format. (RAR 10429 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Chakravorty, D., Banerjee, K., Mapder, T. et al. In silico modeling of phosphorylation dependent and independent cMyc degradation. BMC Bioinformatics 20, 230 (2019). https://doi.org/10.1186/s128590192846x
Received:
Accepted:
Published:
Keywords
 cMyc
 Degradation
 Ubiquitination
 Fbox proteins