|Home | About | Journals | Submit | Contact Us | Français|
Expression of the two transcription factors microphthalmia-associated transcription factor (MITF) and signal transducer and activator of transcription 3 (STAT3) are tightly connected to cell proliferation and survival, and are important for melanocyte development. The co-regulation of MITF and STAT3 via their binding to a common inhibitor Protein Inhibitor of Activated STAT3 (PIAS3) is intriguing. A better quantitative understanding of this regulation is likely to be important for elucidation of the melanocyte biology.
We present a mathematical model describing the MITF-PIAS3-STAT3 signalling network. A default parameter set was developed, partly informed by the literature and partly by constraining the model to mimic reported behavioural features of the system. In addition, a set of experiment-specific parameters was derived for each of 28 experiments reported in the literature. The model seems capable of accounting for most of these experiments in terms of observed temporal development of protein amounts and phosphorylation states. Further, the results also suggest that this system possesses some regulatory features yet to be elucidated.
We find that the experimentally observed crosstalk between MITF and STAT3 via PIAS3 in melanocytes is faithfully reproduced in our model, offering mechanistic explanations for this behaviour, as well as providing a scaffold for further studies of MITF signalling in melanoma.
The melanocytes are skin cells of neural crest origin that constitute 5% - 20% of the basal layer of human epidermis [1-6]. The cell type is responsible for the melanin pigment production and thus the colour patterning of skin and hair in mammals. Melanoma, a cancer originating in melanocytes, is in its later stages notoriously resistant to treatment, and although good prognostic markers exist, the understanding of the underlying biology is only slowly forthcoming . While knowledge about each single protein and gene involved in melanocyte development and regulation of homeostasis is important, developing an understanding of the signalling networks connecting the receptors on the surface to the regulating effect on gene transcription in the nucleus appears crucial in implementing efficient molecular treatment strategies in the dawning era of personalized cancer therapy.
Expression of microphthalmia-associated transcription factor (MITF), the signal transducer and activator of transcription 3 (STAT3), and their co-regulation via protein inhibitor of activated STAT3 (PIAS3), are all tightly connected to cell differentiation, proliferation and survival. MITF is considered to be a master regulatory gene for melanocytes, and has been shown to play important roles in the regulation of genes involved in cell cycle progression, including Bcl-2 and CDK2 [8-10]. MITF is also of clinical significance, as MITF mutations in humans cause Waardenburg syndrome type II , and a significant number of malignant melanomas harbour MITF amplifications. MITF has also been proposed to be important for both differentiation of melanocytes and for tumour transformation . MITF has two phosphorylation sites influencing the PIAS3 binding: S73 and S409. These sites are phosphorylated by different kinases in the MAPK pathway, the ERK and RSK, respectively [13,14].
STAT3 is a transcription factor involved in signal transduction pathways that are activated by several extracellular stimuli, including the IL-6 family of cytokines. It is tyrosine phosphorylated by the Janus kinase (JAK) or SRC. The resulting signal mediates cell growth, differentiation, and survival [15-17]. The underlying molecular details have only partly been elucidated .
PIAS3 has been identified as an inhibitor of both activated STAT3 and MITF [19-23]. PIAS3 can bind activated STAT3, as well as non-activated MITF in one of its two inactive complexes. The phosphorylation of MITF at S409 results in MITF dissociation from the complex, and more PIAS3 is thereby made available. As a result, more STAT3 is bound in complex with PIAS3 and is thus prevented from binding DNA and activating target genes [22,24,25]. Similarly, expression of constitutively active STAT3 will complex with unbound PIAS3, resulting in less PIAS3 being available for binding to MITF. Consequently, more active MITF is observed .
The connection between MITF, STAT3 and PIAS3 (Figure (Figure1)1) has several interesting features: (1) MITF and STAT3 interacts through binding and sequestration of their common inhibitor PIAS3 [19-22], (2) PIAS3 binds to phosphorylated (activated) STAT3, but disassociates from activated MITF  which introduces an asymmetry to the network, (3) MITF has two phosphorylation sites interfering with PIAS3 binding, and all four resulting phosphorylation states have different binding affinities to PIAS3 . We have developed a mathematical model to incorporate quantitative aspects of the system in order to both test if the current conceptions of the system can account for observed results, and to serve as a framework for further studies of this module's interaction with other pathways. See Figure Figure22 for a graphical representation of the model. This dynamic model of the MITF-PIAS3-STAT3 system was designed to be simple, while still being capable of reproducing the available results. The inputs to the model are the activation events of the two transcription factors MITF and PIAS3. When MITF becomes phosphorylated at S73 and/or S409, the affinity to PIAS3, the degradation rate and the transcriptional activity are altered. On the other hand, STAT3 has only one relevant phosphorylation site, which is phosphorylated by JAK. When STAT3 becomes phosphorylated, the affinity to PIAS3 and the transcriptional activity is altered. As output, the activation levels of the two transcription factors are monitored: For STAT3, this implies the amount of phosphorylated protein, while for MITF it represents a complex feature of the distribution of MITF among the four phosphorylation states.
Default values for the set of core parameters common to all experimental settings were obtained partly directly from the literature and partly by fitting model behaviour to available data. These core parameters describe degradation rate, production rate, binding affinity, and phosphorylation and de-phosphorylation rate (Table (Table1).1). For further details on the parameter value deduction, see Methods. We next singled out a representative selection consisting of 28 well-reported experiments addressing various relevant aspects of the MITF-STAT3-PIAS3 system (Table (Table2),2), and tested to what degree the model was able to account for the available experimental data. For each experiment, a set of experiment-specific parameters, such as the effect of transfection on production rates or the activation level in growing cells, was specifically incorporated to reflect the experimental setup and conditions. These parameters were specifically set for each experiment. In those cases where the experiments were replicates, the experiment-specific parameters were given identical values. While keeping the core parameters at default values, we manually searched for experiment specific parameter values that minimized the distance between the simulation output and the reported results of the focal lab-experiment. The experiment-specific perturbations of the model are provided in the Methods section. We found that we were able to account reasonably well for 27 of the 28 experiments. Only experiment #26 seems to be beyond explanatory reach by our model framework.
We next investigated how sensitive the obtained fit between the experimental data and the model results was with respect to variation in the core parameter values, in the 27 successful cases examined. While keeping the experiment-specific parameters unchanged, 106 sets of the core parameters were sampled uniformly on a logarithmic scale from a hypercube centred at the default parameter values and within the range 0.5-2.0× of each default value. For all sampled parameter sets, we simulated the 27 successful experiments in Table Table22 and recorded success or failure accordingly. By analysing the 30 × 106 parameter matrix and the corresponding 28 × 106 result matrix we studied the sensitivity of the success rate of each experiment to variation in a given parameter by computing success rates in bins of sorted parameter sets, and used the sum of absolute deviations from the overall success rate to range the parameters according to sensitivity (Figure (Figure3).3). Each experiment tested different features of the system, which is reflected in the parameter sensitivity. Some experiments are sensitive only to one parameter. For example experiment #1, which investigated the temporal development of ERK and RSK1 kinase activities, and experiment #14, which investigated the transcriptional activity of the MITF-S73D mutant in response to transfection of PIAS3 were both sensitive to a single parameter, the RSK1 phosphorylation rate constant kRp+ and the PIAS3 production rate PPLAS3, respectively. On the other hand, experiment #17, which investigated the MITF-PIAS3 association in response to activation, and experiment #18, which investigated S409A mutated MITF-PIAS3 association in response to activation were both sensitive to more than half the parameters. Further, it was observed that the experiments that dealt only with MITF and the MITF-PIAS3 connection, and not with STAT3 (experiments #2 - #16), were only affected by the parameters regarding MITF. On the other hand, experiments #19 to #24 were sensitive only to perturbations of the STAT3 phosphorylation and de-phosphorylation constants, as well as the level of phosphorylated JAK. Notably, three parameters did not affect the results of any experiment: kMp409ass, kMp409diss, and PMITF. In addition, three parameters had only weak effects: PSTAT3, γSTAT3and ku.
Most of the experiment-specific simulations suggest that the main underlying biological mechanisms seem to be captured by the model. As an example, the result of experiment #25 is plotted alongside the original figure from  in Figure Figure4.4. The temporal development for all variables for all experiment simulations can be found in Additional file 1 - simulationFigures.zip. In the following, we consider only those experiments in need of further interpretation and qualification.
In experiment #7, the MITF-PIAS3 association after activation was investigated. BL6-B16 melanoma cells were co-transfected with MITF and PIAS3. After 48 hours, the cells were treated with tetradecanoyl phorbol acetate (TPA) for 30 minutes. The amount of MITF-PIAS3 complex was measured before, and after 10 and 30 minutes of TPA treatment (Figure (Figure5A).5A). Two different interpretations of the two bands representing MITF have been proposed. In all cases, the lower band is considered as representing un-phosphorylated MITF, while the upper band may represent all phosphorylated states, or alternatively only S73 phosphorylated MITF [14,26]. If we assume that the upper band corresponds to all phosphorylated MITF-PIAS3 complexes (which is in accordance with what is assumed about the bands representing MITF in ), we can predict the temporal development of the distribution among the different phosphorylation states. In our simulations (Figure (Figure5B),5B), the amount of un-phosphorylated MITF-PIAS3 complex decreases rapidly, which is in accordance with the lower band in Figure Figure5A.5A. The sum of all the phosphorylated states is high after 10 minutes, and is falling again after 30 min, which is in accordance with the higher band in Figure Figure5I.5I. The distribution among the phosphorylated states can thus be viewed as a prediction. However, it has been suggested that the upper band represents only S73 phosphorylated MITF , which is also supported by our model. The S73 phosphorylated complex is also high after 10 minutes, and falls again after 30 minutes. If this is the case, the temporal development of the other phosphorylated complexes (MITF-PIAS3p409, MITF-PIAS3pp) can be viewed as predictions of undetected MITF-PIAS3 phosphorylation states.
In experiments #19 - 24, the transcriptional activity of STAT3 was measured in NIH 3T3 fibroblasts that were co-transfected with STAT3, PIAS3, and with different doses (0.05 μg, 0.1 μg) of MITF (wild type [WT] or S409A [MUT]). After 48 hours of incubation, the cells were activated with IL6/IL6R for 6 h. Transcriptional activity of STAT3 was measured, and the mean ± the SEM of three experiments is shown in Figure Figure6A.6A. The simulations of these experiments, which were performed with equal experiment-specific parameters, captured the reported results well (Figure (Figure6B).6B). However, there is a striking, but un-noticed result both in the lab experiment and the simulation in that there is no effect of adding the smaller amount of wtMITF, while there is a rather big effect of adding the greater amount of wtMITF. Because of this similarity in the nonlinearity between the model simulation and the data, an inspection of the model and the simulation results may provide insight to the underlying mechanism. During the incubation time, when the transfected genes are being transcribed and translated to their respective proteins, complexes will be formed continuously. When MITF is added, the PIAS3-MITF complex will form. The sequestering of PIAS3 is not only making it unavailable for STAT3 inhibition, but also rescues PIAS3 from degradation. In the case where the smaller amount of MITF was added, these two effects would have cancelled each other out, while in the case where the greater amount of MITF was added, all PIAS3 would be sequestered, resulting in a saturated rescuing effect.
In experiment #26, a comparison was performed between MITF wild type and a mutant that was transcriptionally active, but unable to bind PIAS3. Transcriptional activity was measured by rtPCR on the mRNA of MITF and STAT3 target genes at 30 minutes and at 4 hours after activation. When the cells were transfected with wild type MITF, elevated transcriptional activity was observed for both transcription factors after activation. However, when the cells were transfected with the MITF mutant, neither MITF nor STAT3 displayed elevated transcriptional activity (Figure (Figure7A).7A). The model was able to mimic the change in transcriptional activity of MITF in reaction to activation for both wild type and mutated MITF (Figure (Figure7B),7B), while the transcriptional activity of STAT3 was not in accordance with the experimental data. In the simulation of this experiment, the transcriptional activity of STAT3 increased in response to activation, regardless of MITF mutation status (Figure (Figure7C).7C). In the lab experiment, this increased transcriptional activity was only seen in the cells with wild type MITF and not in the cells with mutated MITF. The authors of the original paper conclude that PIAS3, being unable to bind the mutated MITF, is accessible for binding and thus inhibition of STAT3. In the model simulations, there was indeed more PIAS3 available in the case with mutated MITF, which causes the lower peak in the level of phosphorylated STAT3 in this case (Figure (Figure7C).7C). However, due to MITF degradation, similar levels of PIAS3 becomes available also in the case of wild type MITF, leading to no difference in level of free PIAS3 to inhibit phosphorylated STAT3 activity. Thus, according to our results, the proposed mechanism is not a plausible explanation of the observed data.
The focal melanocyte regulatory module may reflect the need for co-regulation of two transcription factors. It seems to provide one type of solution for obtaining different co-regulation effects on different time scales, from the immediate effects on the activity of one transcription factor upon activation of the other, to the long term effect of sequestering and later release. Because this module has components working on many different time scales, the history leading up to an event may be determinative. This may be of importance when interpreting results from experiments that operate on different time scales, such as transfection of constitutively active kinases compared to activation of pre-existing kinases.
The asymmetry of the network, introduced by PIAS3 binding to activated STAT3, while most efficiently binding to non-activated MITF, yields a particular regulatory effect. We observe that in the model, STAT3 is more affected by the activation of MITF than MITF is affected by the activation of STAT3. This phenomenon, which is in accordance with earlier observations [19,20], can be elucidated from some features of the model. In the resting cell, PIAS3 is to a larger extent bound to un-phosphorylated MITF and MITFp73, which are the major phosphorylation states in the resting cell. On the other hand, STAT3 only binds PIAS3 when activated. Thus, even if both transcription factors are able to alter the amount of free PIAS3, only STAT3 is affected by it: Activated MITF does not interact with PIAS3, and is therefore unaffected by the amount of free PIAS3. STAT3, on the other hand, binds PIAS3 only in activated state, and the strength of the signal is therefore highly affected by the level of free PIAS3. This is a mechanistic explanation for the observation that STAT3 does not associate with PIAS3 that is already associated with MITF .
Our results indicate that it is not possible to mimic experiment 26 with this model. This finding contradicts the mechanism suggested in  as an explanation for the observed data. In our opinion, these observations must either be caused by un-expected properties in the MITF mutant, and/or by un-modelled connections in the network. The latter is supported by the fact that the module modelled here is an abstraction of a small part of a vastly complex regulatory network. Mechanisms not accounted for in the present model include regulation of MITF [10,27,28] and upstream co-regulators in the JAK-STAT  and the SRC-STAT pathways.
In the sensitivity analysis the model parameters were perturbed in a range from half to double their default value. Some of the parameters exposed negligible effect on the simulation results under these perturbations. For instance there are no experiments that show high sensitivity to the association and dissociation constants kMp409ass, kMp409diss. A possible explanation for this is that the affinity of S409 phosphorylated MITF to PIAS3 is so low that PIAS3-MITFp409 is rare compared to the other PIAS3-MITF complexes, and modest perturbations of the association and dissociation constants does not change this situation. The MITF and STAT3 production rates (PMITF andPSTAT3) are affecting the total amount of protein, but still none of the experiments are particularly sensitive to these parameters. The reason for this seems to be that the relative response measured in the experiments is not affected much by this total amount.
In the skin, the major contributors to melanocyte signalling are dermal fibroblasts and epidermal keratinocytes [30,31]. For MITF, the major regulatory pathways are thought to be the TGF-β and the POMC derivatives ACTH and α-MSH signalling pathways . The TGF-β pathway, which is often dysregulated in melanoma [33,34], down-regulates MITF, while the α-MSH pathway up-regulates MITF expression [28,35,36]. The up-regulation occurs via α-MSH binding and activation of MC1R, which leads to cAMP/PKA activation and activation of MAPK [36,37]. Thus, the effect of the MAPK pathway on MITF is complex, including up-regulation, activation and tagging for degradation. For STAT3, in addition to the IL-6 cytokine family, growth factors such as EGF, PDGF and KITL have been reported to activate STAT3 [38-40]. Here, KITL is of special interest in melanocyte biology. In addition to being essential for melanocyte development, it has recently been reported that the hereto unidentified coupling of MC1R activation to MAPK proceeds via MC1R-coupled KIT activation of MAPK . With a focus on developing potential therapies toward melanoma, both MITF and STAT3 have been implicated in neoplastic progression. MITF has been classified as a bona fide context dependent oncogene . STAT3 has been implicated in neoplastic progression with an increased contribution especially in metastatic melanoma . Strikingly, no spontaneous or inherited mutations have been isolated . Targeting STAT3 is, however, feasible from both upstream and downstream elements. In addition to being a melanoma oncogene, depletion of MITF has been reported to induce arrest, senescence and cell death in melanoma identifying it as a potential target for therapy [43,44]. Indeed, a PIAS3 based peptide has been developed to target both MITF and STAT3 via a PIAS3 23-aa mimic [21,23].
With regards to MITF and STAT3 and their roles in melanocyte biology, it is intriguing that α-MSH signalling is somewhat equivalent to KITL signalling with the additional activation of the cAMP/PKA axis. This leads to an increase in MITF production and dual phosphorylation of both MITF and STAT3. Likewise, α-MSH signalling in melanoma predisposing MC1R RHC (red hair colour) mutant backgrounds, which does not elevate cAMP, would alone be equivalent to KITL signalling in a wild type background with regards to MITF and STAT3, the outcome being STAT3 phosphorylation and MITF activation and subsequent depletion. With this in mind, the emergence of MITF's role in melanoma and its high frequency of duplication, it seems contradictory that MC1R mutants predispose individuals to this disease. It will be interesting to see if compensatory mutations in cAMP elevating factors such as G protein-coupled receptors (e.g. the histamine receptor H2) are found to be necessary for progression and if so to model the outcomes [45-47].
Special conditions like JAK-STAT3 activation duration due to SOCS feedback inhibition versus SRC-STAT3 activation, which is not influenced by SOCS , may be simulated through the model. In addition, other relevant factors such as NF-κ-B, which is also regulated by PIAS3 binding , may be added by incorporating an NF-κ-B module. Another interesting module would, for example, be GM-CSF/KITL or GM-CSF/α-MSH signalling. GM-CSF activates MAPK, but has been reported to inhibit KIT signalling via direct binding of CSF2RA to KIT .
With this model at hand, more detailed issues in this network can be addressed. The previously posed question on which phosphorylation states of MITF are represented by the upper band interpreted to be phosphorylated MITF in Western blot analysis [13,14], is such an issue. With the ability to test different hypothesis in simulations, both the experiment design and the interpretation of the experimental results could benefit.
Further work could include refinement of various parts of the model by integrated wet-lab modelling efforts. The three parts of the model that could benefit the most from structural enhancement, are the following: (i) The mapping function from the MITF phosphorylation states to the MITF transcriptional activity could be represented on a higher resolution level. This is needed for the model to address quantitative data with absolute values. (ii) The production regime of MITF, PIAS3, STAT3 and RSK1 could be represented by both transcription and translation. (iii) The enzymatic equations could be represented by Michaelis-Menten kinetics. All these structural changes will introduce more state variables and model parameters, thus nothing would be gained by these efforts without the generation of accurate quantitative data to pin down both the model structure and the parameter values.
In this work, we have provided a mathematical model of the MITF-PIAS3-STAT3 network and have mimicked a representative selection of lab experiments that explore the features of this network. The analyses of this model have revealed explanations to the observed phenomena, as well as recommended reconsideration of previously proposed explanations. This model provides a framework for further investigation of this interesting crosstalk, and can be used as a tool for experimental design and as starting point for further modelling efforts.
We have developed an ordinary differential equation (ODE) model of the MITF-PIAS3-STAT3 system. The graphical representation of the model given in Figure Figure22 is represented in Systems Biology Graphical Notation as implemented in CellDesigner . The ODE-model was implemented in MATLAB (Additional file 2 - MATLABcode.zip). In the following the model is presented in terms of the chemical reaction equations and the ordinary differential equations. All protein amounts represented by the state variables in the model are presented in an arbitrary unit (au).
The amount of phosphorylated ERK is represented by a constant (a MAPK signal is represented by an altering of this constant). This kinase is involved in two processes; phosphorylation of RSK1 (1) and phosphorylation of MITF at S73 (2).
In all the equations the subscript p after the protein denotes phosphorylated form. The phosphorylation and de-phosphorylation rate constants are represented by kRp+ and kRp- respectively.
All the equations in (2) represent ERKp catalyzed phosphorylation of MITF at S73. Subscripts after MITF denote phosphorylation state (un-phosphorylated, S73 phosphorylated, S409 phosphorylated and phosphorylated at both sites). Complexes are represented by underscore separated protein abbreviations. Note that the same phosphorylation (kMp73E+) and de-phosphorylation (kMp73-) rate constants are used in all four equations. In addition MITF is auto-phosphorylated at S73 governed by a separate rate constant, kMp73a+:
Further, the phosphorylation of MITF at S409 is catalyzed by the phosphorylated RSK1:
The four different phosphorylation states have different association and dissociation rate constants:
STAT3 is also binding PIAS3:
And finally the phosphorylation of STAT3 is catalyzed by phosphorylated JAK (which is represented by a constant):
A signal propagating through the JAK-STAT3-pathway is represented in the model simulations by an elevation of the constant representing the amount of phosphorylated JAK.
These chemical reactions were translated into an ODE model by anticipating simple mass action reaction kinetics. See Table Table22 for a comprehensive list of parameters. Equation (11), describing the dynamics of unbound and un-phosphorylated MITF, is provided in detail below. One differential equation for each phosphorylation and complex state of each protein is given. These equations are in the same way determined from the chemical equations and are therefore only commented when features differ.
In this model, MITF can be present in 8 different states and is therefore represented by 8 differential equations, each containing parts reflecting the chemical reactions involving that particular state.
In differential equation (11), where the rate of change of un-phosphorylated and un-bound MITF is described, PMITF is the production rate, kMp73E+[MITF][ERKp] is the rate of MITF leaving this state because it gets phosphorylated at S73 (reaction equation (2)), kMp73-[MITFp73] is the rate of MITF entering this state from the S73 single phosphorylated state due to de-phosphorylation (reaction equation (2)), kMp73a+[MITF] is the rate of MITF leaving this state because of auto-phosphorylation at S73 (reaction equation (3)), kMp409+[MITF][RSK1p] is the rate of MITF leaving this state because of phosphorylation at S409 (reaction equation (4)), kMp409-[MITFp409] is the rate of MITF entering this state from the S409 single phosphorylated MITF state because of de-phosphorylation (reaction equation (4)), kMass[MITF][PIAS3] is the rate of MITF leaving this state because of MITF-PIAS3 association (reaction equation (5)) and kMdiss[PIAS3_MITF] is the rate of MITF entering this state due to MITF-PIAS3 dissociation (reaction equation (5)). The last term γMITF[MITF], represents the MITF degradation. The total amount of RSK1 is kept constant, while its phosphorylation state is dynamically determined by equations (12) and (13).
The dynamics of S73 phosphorylated MITF (MITFp73) is given in equation (14)
Where γMITFp73 is the degradation rate that applies to the ubiquitinated ratio R of the amount of MITFp73. In the dynamics of S409 phosphorylated MITF,
simple linear degradation is applied. When MITF is phosphorylated on both sites, the sum of the two degradation rates applies:
Further, these four phosphorylation states of MITF bind PIAS3, each with their own association and dissociation rate constants, yielding differential equation (17) for PIAS3, (18) for the MITF-PIAS3-complex, (19) for the S73 phosphorylated MITF-PIAS3-complex, (20) for the S409 phosphorylated MITF-PIAS3-complex and (21) for the double phosphorylated MITF-PIAS3-complex.
PIAS3 has a linear degradation.
STAT3 is binding PIAS3 while phosphorylated. The dynamics of the PIAS3-STAT3-complex is given in equation (22) and un-phosphorylated and phosphorylated STAT3 in equation (23) and (24) respectively.
To represent higher stability by proteins in complexes, the degradation rate of the complex is set to 20% of the mean of the degradation rates assigned to the constituting proteins in un-bound form.
S73 phosphorylated MITF get ubiquitin tagged for degradation by phosphorylated ERK. Introducing another modification site would cause another doubling of the number of state variables representing MITF. To reduce complexity, we chose to represent the ubiquitination with one differential equation (25), describing the ratio R of the total amount of S73 phosphorylated MITF that also is ubiquitinated.
Where A is defined as
The three parts of equation (25) correspond to the three processes altering the ratio R: (i) The ubiquitination of un-ubiquitinated MITF, (ii) the ubiquitination dependent degradation of ubiquitinated MITF and (iii) the production of new un-ubiquitinated MITF. In (i), the ubiquitination rate constant ku, denotes the rate of change of R and thus have the unit au-1 min-1. To derive part (ii), let R' be the ratio R after one minute of ubiquitin mediated degradation. Then the rate of change due to this effect is -(R - R') per minute, where
and B = [MITFp73] + [MITFpp]. The term a = 1 and has the unit minute and does not affect the calculation, but is needed for the correct unit representation. Part (iii) is derived in a similar way: Let R'' be the ratio R after one minute of MITF production. Then the rate of change due to this effect is -(R - R'') per minute, where
The readouts of the model are the amounts of the different phosphorylation and complex states of the proteins involved. To be able to view this model in the light of available experimental results, these levels have to be interpreted in a way comparable to the lab-generated data. The transcriptional activity of MITF and STAT3 is measured either as luciferase activity of a transfected reporter gene or as mRNA-levels of downstream genes.
We thus devised a function f([MITF],[MITFp73],[MITFp409],[MITFpp]) to map the amounts of the different phosphorylation states of MITF to a downstream transcriptional activity. The production of any gene product can be modelled by a function f(S1,S2,....,SN), where Sj is the concentration of the N species affecting the transcription and translation of the particular gene product. This function can be described at various levels of precision from the basic principles of chemical kinetics, via Michaelis-Menten enzymatic kinetics to low order polynomials that function as an ansatz for the actual molecular mechanisms. The function that maps from the amounts of the MITF phosphorylation states to production rate of downstream genes will probably differ between target genes on a higher resolution level. However, for the resolution level of the current effort, one mapping function is sufficient. We have approximated the MITF activity by a sum of first order polynomials, which in effect is a weighted sum:
Values for the transcriptional activity of different MITF mutants are given in . From these we have calculated the factors to be AM = -0.11, AM73 = 0.44, AM409 = 0.11 and AMpp = 0.56. This simple representation of MITF transcriptional activity is able to take negative values, which does not make sense biologically; therefore, A0 is set to the relative high arbitrary value 10, which is sufficiently high to avoid negative MITF activity in this study.
While this model contains a large number of parameters, we can make at least a reasonable order of magnitude estimate for most of them from empirical biological data. In the experiments used for model development and parameter fitting, neither implicit nor explicit absolute protein amounts have been provided. The cells are either transfected, activated or both, and thus the input protein amounts and perturbations will by this fact be qualitative. The results are presented as qualitative protein or complex amounts (Western blot) or transcriptional activity for MITF and PIAS3 (measured by luciferase assays or PCR-measurements of target genes mRNA levels). These results are quantified relative to a control, but not in absolute terms.
An estimate for the degradation rates of PIAS3 and STAT3 is calculated directly from the protein stability index (PSI) as presented in  (γSTAT3 = 0.002 and γPIAS3 = 0.008). The PSI derived degradation rate for MITF is 0.008, but MITF is known to be relatively stable in un-phosphorylated form while ERK mediated phosphorylation on S73 also tags MITF for proteasome mediated degradation [13,14]. To match degradation time series data [13,14] MITF was assigned a low degradation rate in the un-phosphorylated state, γMITF = 0.0012, and relatively high degradation rate when phosphorylated on S409 (regardless of phosphorylation status on S73), γMITFp73 = 0.01. The proportion of MITF phosphorylated on S73 (regardless of phosphorylation status on S409) is degraded only when ubiquitinated (and then with a rather high degradation rate, γMITFp73 = 0.02). The ubiquitination rate constant is adjusted to ku = 0.0001 to yield degradation in accordance with data [13,14]. When the proteins are in complex, they are thought to be more stable than as free molecules. As an approximation the degradation rates of all complexes were set to 20% of the average of the degradation rate of the two components.
The production rates for MITF, PIAS3 and STAT3 are complex functions of the concentration of a large number of substances, of which none are involved in this model. Therefore these production rates are represented here by constants whose values are relative to the corresponding degradation rates determined by the steady state levels of each of the three proteins. It is important for the proper function of this module that these three proteins are present at approximately equimolar levels. The values are therefore set to pMITF = 1, pPLAS3 = 1.262 and pSTAT3 = 0.211 to meet this criterion. This yields a steady state amount of approximately 100 (au). The level of phosphorylated ERK in the resting cell is set to ERKp = 10 and is elevated to approximately 1000 whenever a MAPK signal is simulated. The total amount of RSK1 is set to 500 of which the phosphorylation is determined dynamically (equation (12) and (13)). The phosphorylation speed for RSK1 is found in , and the phosphorylation and de-phosphorylation rate constants calculated to kRp+ = 0.0004 and kRp- = 0.04, respectively. The Western blot analysis of MITF presented in  and  has one band interpreted as un-phosphorylated MITF and one band interpreted as phosphorylated MITF. No interpretation is provided on which phosphorylation states this latter band may represent. By interpreting this band as the sum of all phosphorylated MITF, we have adjusted the phosphorylation rate constants of MITF to make the model behaviour fit the experimental data. The rate constants are set to kMp73E+ = 0.00015, kMp73- = 0.03, kMp73a+ = 0.025, kMp409+ = 0.0001 and kMp409- = 0.04. With no direct measurements of the STAT3 phosphorylation rate, we assigned the same values to the STAT3 phosphorylation/de-phosphorylation rate constants as we did for RSK1: kSp+ = 0.0002 and kSp- = 0.04. The constant representing phosphorylated JAK is set to 10 in the un-activated state and approximately 1000 in the activated state.
With no direct measurements of the MITF-PIAS3 and STAT3-PIAS3 association and dissociation rate constants, they were inferred from qualitative results and assumptions. We anticipated that this process is faster than the phosphorylation process, which means that both association and dissociation rate constants should be relatively high. MITF has lower affinity to PIAS3 when phosphorylated at S409 and higher affinity when phosphorylated at S73, compared to un-phosphorylated MITF . The double phosphorylated MITF exhibits an intermediate affinity comparable to that of un-phosphorylated MITF . When the MITF-PIAS3 association/dissociation rate constants were adjusted so that the model emulated the results from , they were assigned the following values: kMiss = 0.01, kMdiss = 1, kMp73ass = 0.03, kMp73diss = 0.5, kMp409ass = 0.0001, kMp409diss = 1, kMppass = 0.01 and kMppdss = 1. The STAT3-PIAS3 association/dissociation rate constants were adjusted in the same way to data from  and are assigned the values kSpass = 0.005 and kSpass = 0.2.
We extracted a set of representative experiments from four relevant publications [14,19,20,22] (Table (Table2).2). The numbering in the table provided a unique ID for each experiment, and was used throughout the article. The perturbations of the cells that were reported in the published experiments, like transfection of mutated genes or receptor activation, were assigned digital counterparts in the model through alterations of the input values, starting values, or some of the parameters. In the following, these details, as well as the criteria used for the sensitivity analysis, are given for each of the experiments.
Experiment # 1 is a simulation of the experiment presented in  Figure Figure5B.5B. Here, the authors use kinase assays to investigate the temporal development of ERK and RSK1 kinase activity after stimulation of melanoma cells. The stimulation was simulated by elevation of the amount of phosphorylated ERK to 1000 from the default value of 10. The success criterion used in the sensitivity analysis was that the phosphorylated RSK1 level should reach 90% of its maximum between 3 and 10 minutes after activation.
Experiment #2 simulates the pre-stimulation distribution among the MITF phosphorylation states in the experiments presented in , Figure Figure11 and and2.2. To emulate the growing cell culture the level of phosphorylated ERK was elevated from the default value of 10 to the intermediate level of 80. Steady state values were found for the state variables and the total amount of phosphorylated and un-phosphorylated MITF was compared. The success criterion used in the sensitivity analysis was that both categories should be between 25% and 75% of the total amount.
Experiment #3 is also based on , Figure Figure11 and and2.2. Here, the phosphorylation state of MITF after 30 minutes of activation is considered. The cell stimulation is simulated by elevation of the levels of phosphorylated ERK and JAK from the default value of 10 to 1000. The success criterion used in the sensitivity analysis was that more than 80% of the MITF is phosphorylated after 30 minutes.
Experiment #4 is also based on , Figure Figure11 and and2.2. Here, the degradation profile is considered. The cell stimulation is simulated by elevation of the levels of phosphorylated ERK and JAK from the default value of 10 to 950. The success criterion used in the sensitivity analysis is that less than 50% of the MITF is degraded after 1 h and at least two thirds are degraded after 5 h.
Experiment #5 is a simulation of the experiment presented in , Figure Figure6.6. Here, the authors have transfected MITF, a MITF reporter gene, and different amounts of PIAS3 into NIH 3T3 fibroblasts (endogenously not expressing MITF). To mimic MITF transfection, the MITF production rate is increased from 1 to 18. To mimic the four different amounts of PIAS3 transfection, the PIAS3 production rate was increased from 0.211 to 0.5, 1, 2 and 4. The model was simulated for 2880 minutes (2 days) for each perturbation and the resulting MITF activity was compared for the four PIAS3 levels. The success criterion used in the sensitivity analysis was that MITF activity should be monotonically decreasing with increasing PIAS3 and that the MITF activity at the highest PIAS3 amount should be less than the half that of the lowest PIAS3 amount.
In experiments #6 and #7, we have mimicked the experiments presented in , Figure Figure2.2. Here the authors have investigated the association between MITF and PIAS3 in BL6-B16 cells before and after 10 and 30 minutes of activation with c-kit ligand. In experiment #6 only GFP-PIAS3 was transfected, while in experiment #7, the cells were transfected with both MITF and PIAS3. In experiment #6, activation is performed by increase of the levels of phosphorylated ERK and JAK from 10 to 500 and 50, respectively. In experiment #7, transfection was simulated by elevation of the MITF production rate from 1 to 10 and the PIAS3 production rate from 1.262 to 10 and simulated for 2880 minutes (2 days). Thereafter, activation was simulated by increase of the amount of phosphorylated ERK and JAK from 10 to 1200 and 250, respectively. The success criterion used for the sensitivity analysis in both experiments was that the total amount of MITF-PIAS3-complex should be higher after 10 minutes of stimulation compared to before stimulation and that the amount of complex should be lower after 30 minutes of stimulation than after 10 minutes.
Experiment #8 simulated the experiment from , Figure Figure3A,3A, left. The authors co-transfected NIH 3T3 cells using MITF, with a constitutively active RSK1 plasmid and with a reporter gene to read MITF activity. To simulate transfection of MITF, the MITF production rate was elevated from 1 to 7. The amount of RSK1 is not dynamically determined in the model, and thus the transfection of this product is approximated with an elevation of the static amount from 500 to 5000. The manipulated model was simulated for 2880 minutes (2 days), and the MITF activities with and without RSK1 were compared. The success criteria used in the sensitivity analysis was that the two cases should not differ more than twofold.
Experiments #9 to #12 represent simulations of the experiment presented in  Figure Figure3B.3B. Here, the authors transfected NIH 3T3 cells with MITF or the MITF-S409A mutant, with constitutively active RSK1 and with a reporter gene to monitor MITF transcriptional activity. To simulate the transfection of MITF, the MITF production rate was elevated from 1 to 7. To simulate the transfection of PIAS3, the PIAS3 production rate was elevated from 1.262 to 4. The amount of RSK1 is not dynamically determined in the model, and thus the transfection of this product is approximated with an elevation of the static amount from 500 to 5000. The activation level was elevated from 10 to 20 to emulate the activity in growing cells. After simulation for 2880 minutes (2 days), the MITF activities, with and without PIAS3 were compared for each of the four experiments. For each experiment the manipulations were as follows: without RSK1, with wild-type MITF (experiment #9), with RSK1 and wild-type MITF (experiment #10), without RSK1, with mutant MITF (experiment #10) and with RSK1, with mutant MITF (experiment #11). The success criterion used in the sensitivity analysis was for each experiment: Experiment #9: MITF activity in the case with PIAS3 is less than 42% and greater than 2% of that without PIAS3. Experiment #10: MITF activity in the case of PIAS3 is less than 84% and greater than 44% of that without PIAS3. Experiment #11: MITF activity in the case of PIAS3 is less than 55% and greater than 15% of that without PIAS3. Experiment #12: MITF activity in the case of PIAS3 is less than 50% and greater than 10% of that without PIAS3.
Experiment #13 to #16 are simulations of the experiment presented in  Figure Figure6.6. In this experiment the authors investigate the inhibition of transcriptional activity of S/D mutants of MITF by PIAS3. This is achieved by transfection of wild-type MITF (experiment #13), MITF-S73D (experiment #14), MITF-S409D (experiment # 15) or MITF-S73/409D (experiment #16) and PIAS3 and a reporter gene into NIH 3T3 cells. We simulated the MITF transfection by elevation of the MITF production rate from 1 to 1.5, and the PIAS3 transfection was simulated by elevation of the PIAS3 production rate from 1.262 to 4.5. Here we have also set all MITF levels to zero before starting the simulation, since NIH 3T3 fibroblasts do not express endogenous MITF. To simulate MITF-S73D mutation we have set the MITF S73 de-phosphorylation rate constant to zero while the MITF S73 auto-phosphorylation rate constant is elevated from 0.025 to 5. To simulate MITF-S409D mutation, the RSK1 phosphorylation rate constant is elevated form 0.0004 to 0.04, the RSK1 de-phosphorylation rate constant is decreased from 0.04 to 0.004, the MITF S409 phosphorylation rate constant is elevated from 0.0001 to 0.01 and the MITF S409 de-phosphorylation rate constant is set to zero. The model was simulated for 2880 minutes (2 days) and the MITF transcriptional activity with PIAS3 was compared with the case without PIAS3 for wild-type MITF, and for each mutation form. The success criterion used for the sensitivity analysis was for each experiment: Experiment #13: MITF activity in the case of PIAS3 is less than 70% and greater than 30% of that without PIAS3. Experiment #14: MITF activity in the case of PIAS3 is less than 50% and greater than 10% of that without PIAS3. Experiment #15: MITF activity in the case of PIAS3 is less than 115% and greater than 75% of that without PIAS3. Experiment #16: MITF activity in the case of PIAS3 is less than 85% and greater than 45% of that without PIAS3.
Experiments #17 and #18 are simulations of the experiment presented in  Figure Figure3A.3A. Here, the authors investigate PIAS3-STAT3 association in response to stimulation of NIH 3T3 cells transfected with STAT3, PIAS3 and MITF or the MITF S409A-mutant. The transfection is simulated by elevation of the MITF production rate from 1 to 10, the PIAS3 production rate from 1.262 to 10 and the STAT3 production rate from 0.211 to 10. The MITF S409A mutation was simulated by setting the MITF S409 phosphorylation rate constant to zero. The activation was simulated by elevation of the amount of phosphorylated ERK and JAK from 10 to 1000. In experiment #17, wild-type MITF, PIAS3 and STAT3 were all transfected and simulated for 2880 minutes (2 days), thereafter simulated for 15 more minutes with and without stimulation. The amount of PIAS3-STAT3-complex was compared. The success criterion used in the sensitivity analysis was that the amount of PIAS3-STAT3-complex should be at least 25% higher with stimulation compared to the case without stimulation. In experiment #18, S409A MITF, PIAS3 and STAT3 were transfected and simulated for 2880 minutes (2 days), and thereafter simulated for 15 more minutes with or without stimulation. The amount of PIAS3-STAT3-complex was compared with the stimulated case in experiment #17. The success criterion used in the sensitivity analysis was that the amount of PIAS3-STAT3-complex in the stimulated case in experiment #17 should be at least 25% higher than the amount of PIAS3-STAT3-complex in any of the two cases in #18.
Experiments #19 to #24 are simulations of the experiment presented in  Figure Figure3B.3B. Here, the authors have investigated STAT3 transcriptional activity as a response to activation, transfection of PIAS3, and transfection of various amounts of wild-type and S409A mutated MITF in NIH T3T cells. The activation was simulated by elevation of the amounts of phosphorylated ERK and JAK from 10 to 1000. The transfections of PIAS3, STAT3 and the two different MITF amounts were simulated by elevation of the PIAS3 production rate from 1.262 to 7, the STAT3 production rate from 0.262 to 5, and the MITF production rate from 1 to 10 or 50, respectively. The S409A mutation was simulated by setting the MITF phosphorylation rate constant to zero. The model was run for 2880 minutes to simulate the incubation and then for 360 minutes to simulate the activation. The background luciferase activity of STAT3 was determined by elevation of the STAT3 production rate and running a model simulation for 3240 (2880+360) minutes. The amount of phosphorylated STAT3 in the end was interpreted as an estimate for the luciferase activity of STAT3 without activation and without PIAS3 or MITF. For experiments #19 to #24, the amount of phosphorylated STAT3 was compared to this background level. In experiment #19, the cells were transfected with STAT3 and activated. The success criterion used in the sensitivity analysis was the level of phosphorylated STAT3 between 10 and 20 times the background level. In experiment #20, the cells were transfected with STAT3 and PIAS3 and activated. The success criterion used in the sensitivity analysis was the level of phosphorylated STAT3 between 2.67 and 5.33 times the background level. In experiment #21 the cells were transfected with STAT3, PIAS3 and the smaller amount of MITF and activated. The success criterion used in the sensitivity analysis was the level of phosphorylated STAT3 between 2.67 and 5.33 times the background level. In experiment #22 the cells were transfected with STAT3, PIAS3 and the larger amount of MITF and activated. The success criterion used in the sensitivity analysis was the level of phosphorylated STAT3 between 6.67 and 13.33 times the background level. In experiment #23 the cells were transfected with STAT3, PIAS3 and the smaller amount of S409A mutated MITF and activated. The success criterion used in the sensitivity analysis was the level of phosphorylated STAT3 between 5.67 and 11.33 times the background level. In experiment #24 the cells were transfected with STAT3, PIAS3 and the larger amount of S409A mutated MITF. The success criterion used in the sensitivity analysis was the level of phosphorylated STAT3 between 10 and 20 times the background level.
Experiment #25 simulated the experiment presented in  Figure Figure4B.4B. Here, the authors have investigated MITF transcriptional activity in response to transfection of PIAS3, STAT3-Y705F mutant and a constitutively active STAT3 mutant (STAT3-C) in NIH 3T3 cells. The transfections were simulated by elevation of the production rate of MITF from 1 to 2, PIAS3 from 1.262 to 2 and STAT3 from 0.211 to 20. The NIH 3T3 cells were simulated by decreasing the starting values of all MITF and PIAS3 states to 1 and in simulations where MITF or PIAS3 were not transfected; their production rate were set to zero. The STAT3-Y705F mutant was simulated by setting the STAT3 phosphorylation rate constant to zero. STAT3-C was simulated by increase the STAT3 phosphorylation rate constant from 0.0002 to 0.02 and setting the de-phosphorylation rate constant to zero. Six different simulations were performed and the MITF transcriptional activity was compared: (i) Cells transfected with STAT3-Y705F, (ii) cells transfected with STAT3-C, (iii) cells transfected with MITF and STAT3-Y705F, (iv) cells transfected with MITF and STAT3-C, (v) cells transfected with MITF, PIAS3 and STAT3-Y705F and (vi) cells transfected with MITF, PIAS3 and STAT3-Y705F. The success criterion used in the sensitivity analysis was that the MITF transcriptional activity from simulation (i) should be between 50% and 150% of the MITF activity form (ii) and the MITF activity from (iii) should be between 80% and 120% of the MITF activity from (iv) and the MITF activity from (v) should be between 12.5% and 100% of the MITF activity from (vi) and the MITF activity from (i) should be less than the MITF activity from (iii) and the MITF activity from (iii) should be greater than 200% of the MITF activity from (v).
Experiment #26 simulates the experiment presented in  Figure Figure6D6D and and6E.6E. Here, the authors use mRNA levels of MITF and STAT3 target genes to investigate MITF and STAT3 transcriptional activity after stimulation of mast cells derived from wild type or mutant mice. The MITF protein from the mutated mice lacks its ability to bind PIAS3. The activation is simulated by elevation of the amounts of phosphorylated ERK and JAK from 10 to 1000. The mutation is simulated by setting all four MITF-PIAS3 association rate constants to zero. The success criterion used in the sensitivity analysis was for cells from wild type mice: MITF activity after 30 minutes should be more than twice the MITF activity without stimulation, STAT3 activity at 4 hours should be more than twice the STAT3 activity without stimulation and for cells from mutated mice: MITF activity at 30 minutes and at 4 hours should be less than twice the MITF activity without stimulation and STAT3 activity at 30 minutes and 4 hours should be less than twice the STAT3 activity without stimulation.
Experiments #27 and #28 are simulations of the experiments presented in  Figure Figure7C7C and and7D.7D. In experiment #27 the authors have investigated the STAT3 transcriptional activity in response to transfection of wild type MITF and a MITF mutant emulating MITF phosphorylated at S409. In experiment #28 the authors have investigated MITF transcriptional activity in response to transfection of constitutively active STAT3 (STAT3-C). The transfections were simulated by elevation of the STAT3 production rate from 0.211 to 5 and the MITF production rate from 1 to 5. The MITF mutant was simulated by setting the MITF-S409 phosphorylation rate constant to 5 and the MITF-S409 de-phosphorylation rate constant to zero. The STAT3-C mutant was simulated by increasing the STAT3 phosphorylation rate by 25 times and setting the de-phosphorylation rate to zero. The model was simulated for 2880 minutes without and with each of the three transfections. The success criterion used in the sensitivity analysis was for experiment #27 that the STAT3 activity had increased when the cells were transfected with MITF or the MITF mutant compared to without transfection and that the increase was higher when cells were transfected with wild type MITF compared to mutated MITF. The success criterion used in the sensitivity analysis for experiment #28 was MITF activity increased by at least 10%.
The ability of the model to mimic the experiments described above was tested for values of the core parameters in the multidimensional neighbourhood of the default values. Each parameter was sampled from a distribution being uniform on the logarithmic scale, with a sampling area reaching from one half to double the default value. For each set of values of the core parameters, the experiments were simulated and the result (success or failure) was recorded for each experiment. This procedure was repeated 106 times resulting in 106 parameter sets, each with a corresponding result vector recording success or failure for each experiment.
In order to study the sensitivity of each experiment i's success rate to variation in a parameter j, we sorted all parameter sets according to parameter j and divided them into 100 bins of 10000 parameter sets each. We used fi , the o verall success rate of experiment i (see Table Table2,2, column 6) and , the success rate in bin k to computed the sensitivity measure . In order to remove clearly insensitive parameters, we set to zero all si, j that did not exceed the maximum value observed under 10000 permutations of the parameter values.
JT, SWO and EH conceived the paper. JT designed the model, generated the parameter set, performed the programming and analysis of the model, and drafted the manuscript. TJL, ABG and SWO participated in model definition and revisions. TJL and EH contributed with knowledge on the biological scope, and with parameter estimation. ABG and JT performed the sensitivity analysis. All authors revised the manuscript. All authors read and approved the final manuscript.
A zip-file containing temporal plots for all variables for all simulated experiments with default parameters. The figures are in jpg-format.
A zip-file containing an MATLAB implementation of the model and scripts used to carry out this study. In addition to the MATLAB scripts, a readme file is provided with instructions on how to get the model running.
This study was supported by The National Programme for Research in Functional Genomics (FUGE) in the Research Council of Norway