 Research article
 Open Access
 Published:
Timedependent antagonistagonist switching in receptor tyrosine kinasemediated signaling
BMC Bioinformaticsvolume 20, Article number: 242 (2019)
Abstract
Background
ErbB4/HER4 is a unique member of the ErbB family of receptor tyrosine kinases concerning its activation of antiproliferative JAK2STAT5 pathway when stimulated by ligand Neuregulin (NRG). Activation of this pathway leads to expression of genes like βcasein which promote cell differentiation. Recent experimental studies on mouse HC11 mammary epithelial cells stimulated by ligand Neuregulin (NRG) showed a timedependent switching behavior in the βcasein expression. This behavior cannot be explained using currently available mechanistic models of the JAKSTAT pathway. We constructed an improved mechanistic model which introduces two crucial modifications to the canonical HER4JAK2STAT5 pathway based on literature findings. These modifications include competitive HER4 heterodimerization with other members of the ErbB family and a slower JAK2 independent activation STAT5 through HER4. We also performed global sensitivity analysis on the model to test the robustness of the predictions and parameter combinations that are sensitive to the outcome.
Results
Our model was able to reproduce the timedependent switching behavior of βcasein and also establish that the modifications mentioned above to the canonical JAKSTAT pathway are necessary to reproduce this behavior. The sensitivity studies show that the competitive HER4 heterodimerization reactions have a profound impact on the sensitivity of the pathway to NRG stimulation, while the slower JAK2independent pathway is necessary for the late stage promotion of βcasein mRNA transcription. The difference in the time scales of the JAKdependent and JAKindependent pathways was found to be the main contributing factor to the timedependent switch. The transport rates controlling activated STAT5 dimer nuclear import and βcasein mRNA export to cytoplasm affected the time delay between NRG stimulation and peak βcasein mRNA activity.
Conclusion
This study highlights the effect of competitive and parallel reaction pathways on both short and longterm dynamics of receptormediated signaling. It provides robust and testable predictions of the dynamical behavior of the HER4 mediated JAKSTAT pathway which could be useful in designing treatments for various cancers where this pathway is activated/altered.
Background
ErbB4/HER4 belongs to the ErbB receptor tyrosine kinase family (consisting of EGFR, ErbB2, ErbB3, and ErbB4). These receptors regulate several critical signaling pathways that are frequently altered in cancers of lung, breast, prostate, etc. When activated by specific growth factors, they initiate multiple signaling cascades leading to the transcription of genes responsible for the determination of cell fate such as proliferation, differentiation or apoptosis [1]. Overexpression of these receptors or development of domainspecific mutations that allow them to be constitutively activated, causes them to promote various important pathways that drive the cell towards a program of proliferation and suppresses those that lead to apoptosis (cell death) or growth arrest through cell senescence. Compared to the other members of the ErbB family, the role of HER4 is incompletely understood [2]. In part, this has to do with the fact that this receptor has several unique features compared to other members of the family. It has four structurally and functionally different isoforms generated by mRNA splicing. Some of these isoforms undergo shedding of their ectodomain by proteolytic cleavage reaction mediated by Tumor Necrosis Factoralpha converting enzyme (TACE). The remaining 80 kDa fragment is further cleaved by enzyme γsecretase and is transported to the nucleus and other parts of the cell thereby taking part in critical cellular reactions that affect cell fate [3]. Unlike the other members of the ErbB family, HER4 has an antiproliferative role through its activation of the JAK2STAT5 pathway. Ligand Neuregulin (NRG), which is one of the ten natural ErbB family ligands and is often expressed in human tumors, can bind to and activate HER4. Upon activation and binding with JAK2, activated HER4 s80 domain can activate Signal Transducer and Activator of Transcription 5 (STAT5) causing it to translocate to the nucleus and initiate transcription of genes that promote differentiation [2, 3]. Hence it is essential to gain a mechanistic understanding of this pathway, particularly its dynamical behavior in response to NRG. Constructing and solving systems of nonlinear ordinary differential equations is a typical approach of modeling these pathways. The models are constructed based on experimental measurements of protein activities and rate constants and are used to predict the timedependent dynamics of target gene/protein expression for different ligand stimulations. The predictions are validated through experimental measurements of the timedependent protein activity through Western Blots or Immunoprecipitation studies [4]. Such an approach has been used for HER4 mediated JAKSTAT pathway before [5, 6]. However, these models failed to reproduce some interesting recent observations from experiments conducted in HC11 mouse mammary epithelial cell lines which prompted this computational study. In the experiments, it was observed that HER4 mediated transcription of antiproliferative genes like βcasein milk genes follows a timedependent switching behavior which cannot be explained by the previously published models. Here we propose an improved model of the HER4JAKSTAT pathway that incorporates additional interactions which have been previously reported in the literature but were not part of the original models. Some of these reactions include competitive heterodimerization of HER4 receptors with other members of ErbB family and JAK2 independent activation of HER4. The latter phenomenon of signaling through parallel pathways has some similarities with signaling pathways driven by constitutively activated mutant forms of EGFR [7]. Identification of such pathways can be critical to gain a better understanding of the HER4 induced JAK/STAT pathway and leverage its anticancer role by designing an appropriate treatment strategy to treat cancers in cell lines where this receptor is significantly expressed.
Results
HER4 signaling through JAKdependent pathway
The computational study was motivated by some intriguing results from experiments conducted on HC11 mouse mammary epithelial cells [8] stimulated with different ligands including NRG. In these studies, HC11 cells were grown and maintained at 5% CO2 in complete medium (RPMI 1640, 10% FBS, 2 mM Lglutamine, 100 U/mL penicillin, 100 μg/mL streptomycin, 1 μg/mL hydrocortisone, 10 ng/mL murine EGF and 5 μg/mL insulin). Cells were seeded in 6well plates at a density of 2 × 104 cells/cm^{2} and allowed to grow to 100% confluence to induce differentiation. The cells were maintained at confluence for 2 days in serumfree/EGFfree medium to induce competence. The competent cells were then stimulated with NRG to induce differentiation. Three stages of the differentiation process were identified: Stage 1 refers to growing the cells to confluence, Stage 2 refers to maintaining the cells at confluence for 48 hours to induce a state of competence, and Stage 3 refers to cell stimulation with HER4 ligand. RTPCR and ELISA based Transcription Factor activation assay were performed in differentiating HC11 cells to characterize the signaling dynamics in the HER4mediated JAK/STAT pathway. Specifically, the following readouts of HC11 mammary differentiation were assayed:

Levels of activated STAT5A and glucocorticoid receptor (GR) in the nucleus

Expression of the βcasein milk gene mRNA.
STAT5 and GR are the two key transcription factors which are both activated during HC11 cell differentiation and synergize on the βcasein gene promoter [9]. In the experimental studies, the steroid hydrocortisone (HC), which signals through GR was included in the medium. It was shown that in cells stimulated with NRG in the absence of HC in the medium, no βcasein was expressed. However, upon addition of HC, HER4 induced robust expression of βcasein. HC was required in Stage 2 of the differentiation process as well as Stage 3. Hence the βcasein expressions obtained by varying NRG levels were normalized by the control case of only HC present in the medium.
In tracking the response of the HC11 cells to HER4 stimulation by NRG in the presence of HC at various time points, it was observed (Fig. 1) that at early time points poststimulation, NRG does not exhibit any effect on βcasein expression. However, at 12h poststimulation, it was found that increased stimulation with NRG decreased levels of βcasein, relative to the control (hydrocortisone). At 24–48 h poststimulation, NRG began to increase βcasein expression levels at sufficiently high ligand concentrations. This result has not been discovered previously, as the few studies examining the HER4STAT5A pathway focused mostly on gene expression at earlier time points post stimulation [9, 10].
As a first step towards understanding the mechanism of this timedependent switching behavior observed in the experiments, we used mathematical models of JAKSTAT pathways reported before in literature [5, 11]. In these models, HER4 is activated in a JAK2dependent fashion, and the HER4JAK2 dimer activates cytoplasmic STAT5 which dimerizes and translocates to the nucleus. This model was run with increasing levels of ligand (NRG) stimulation in the same range as used in the experiments [8]. Although these models were developed for shorter time ranges of 0–12 h we wanted to see to what extent this model can capture the early time behavior observed in experiments. In Fig. 2a and b, we plot the timeintegrated βcasein gene expression as a function of ligand NRG stimulation and instantaneous time profile of βcasein gene. To show how model parameter uncertainties influence its predictions we sampled the parameter space using Latin Hypercube Sampling (see Methods) after setting biologically proper bounds on each parameter. This sampling was used to generate an ensemble of model parameter values which were used to run the model and obtain a set of output values. The average value of this set is reported, and the average dispersion is shown in the form of error bars. All the timeintegrated plots in Figs. 2, 3, 4 and 5 are generated using this procedure.
The plots corresponding to the literature model show little sensitivity towards NRG in terms of βcasein expression at all three time intervals. The literature model has several negative regulators including two cytoplasmic phosphatases for JAK2 and STAT5, one nuclear phosphatase for STAT5 and negative feedback through SOCS protein deactivating HER4JAK complex. SOCS is another product of STAT5 mediated transcription apart from βcasein. The latter is expected to limit the response and prevent any increase in transcription due to higher ligand stimulation. However, the insensitivity of βcasein towards NRG is inconsistent with the observations from the experiments. In the 0–12 h interval, the experiments show a decrease in the βcasein gene expression with increasing ligand stimulation. This observation suggested that there are additional mechanisms which are not considered in the canonical JAKSTAT pathways published earlier in the literature.
Competitive heterodimerization involving HER4 upon ligand binding
It has been previously reported in the literature that competitive heterodimerization of HER4 with other members of ErbB family can result in a decrease in HER4 activity for a particular reaction pathway with increasing ligand stimulation (socalled partial agonists) [12,13,14,15]. To explore this effect, we introduced additional heterodimerization reactions of HER4 with ErbB3. Introduction of such reactions resulted in a decrease in the βcasein activity with increasing NRG stimulation consistent with experimental observations. The reduction is evident from Fig. 2c and d where we have again plotted both integrated and instantaneous βcasein mRNA concentrations. Since there are several competing reactions possible for a fixed number of free HER4 receptors at the cell surface, increasing ligand stimulation promotes the more favorable (heterodimerization) reaction at the expense of the less favorable (HER4JAK dimerization). Addition of competitive heterodimerization, however, does not increase βcasein activity in the 24–48h interval as observed in the experiments. This observation suggests additional reactions are at play during this late period contributing to an increase in the βcasein gene expression.
Combined model with JAKdependent and JAKindependent HER4 signaling
It has been reported before in the literature [16, 17] that HER4 can activate STAT5 independent of JAK2 although with lower rates compared to the canonical JAKdependent pathway. We hypothesized that such a slower JAKindependent pathway might likely become active in the later stages (24–48 h) where it can contribute to an increase in gene expression with increasing ligand stimulation after the faster and more favorable heterodimerization reactions have equilibrated. Also, this pathway is not affected by negative feedback through SOCS, which can only act on the JAKbound complex. We incorporated such a pathway in our model to explore late time behavior (see the method section for more details). The combined pathway is shown in the schematic diagram in Fig. 3 where the JAKdependent and JAKindependent pathways are highlighted.
When this combined pathway was stimulated with increasing NRG stimulations (plotted in Fig. 4), we were able to reproduce the timedependent switching behavior seen in the experiments: that is, in the 0–24 h period the βcasein activity decreases with ligand stimulation, and it increases in the later stages (24–48 h). From the βcasein mRNA time profiles, we see that after the 24h mark there is a transition from negative to positive dependence of βcasein mRNA towards ligand stimulation. In this regime, the βcasein levels continue to increase.
Time delay in βcasein mRNA transcription
In the HC11 cell line experiments mentioned above additional ELISAbased transcription factor (TF) activation assays were also performed to measure the activity of STAT5A and Glucocorticoid receptor (GR) which are the two transcription factors necessary for transcription of βcasein mRNA. It was observed that even though STAT5A and GR activity was highest in the nucleus 15–30 min post ligand stimulation, a significant activity of these transcription factors persisted even 24 h post stimulation. These findings were consistent with previous ChIPSeq studies [18, 19] which assayed for binding of various TFs (including STAT5A and GR) to the βcasein gene promoter at different time points following stimulation with prolactin (PRL). These studies showed that although several of the TFs assemble on the promoter at early time points, the RNA polymerase does not bind and commence transcription until 24 h poststimulation.
To incorporate this delay in the transcription and detection of βcasein mRNA, we introduced two significant transport rates in our model — nuclear import of activated STAT5 dimer and nuclear export of βcasein mRNA. The plots of instantaneous βcasein mRNA and activated STAT5 dimer at the nucleus between 0 and 12 h time points for two different mRNA transport rates show that with low rates of mRNA export we get a significant delay in the mRNA peak activity. These results are consistent with previous experimental and modeling studies [5, 17] which identified these transport rates as essential determinants of the signaling activity of this pathway. Also, such delay in transcription did not affect the timedependent switching behavior (Fig. 5e and f) in the βcasein expression which occurs over a much longer time scale.
Additional file 1: Figures S5S14 explore the effect of the reaction parameters such as initial receptor numbers, phosphatase activity, transcription rates, HER4 homo, and heterodimerization rates and the rates of JAKindependent activation of STAT5 by HER4 on the βcasein expression in the 0–12, 12–24 and 24–48 h time intervals. These results complement the findings from the Global Sensitivity Analysis (next section).
Global sensitivity analysis
To systematically explore the effect of model parameters on various output quantities of interest (such as βcasein expression, timedependent switch in βcasein expression, βcasein transcription time delay, etc.), we performed sensitivity analysis. The HER4JAKSTAT signaling network is a highly nonlinear system, and most of the parameters can vary over a wide range in the corresponding biological system. Hence, a Global Sensitivity Analysis (GSA) which considers the combined effect of the model parameters rather than one at a time is a more appropriate method here [20]. We use a particular type of global sensitivity analysis called Sobol Sensitivity Analysis [21] method for this study (details in the Methods section). The analysis procedure involves:

Sampling over the space of model parameters to create an ensemble of models.

Running the simulations for each member of the ensemble.

Calculating the output quantity of interest from the results for each member of the ensemble.

Determining the effect of the model parameter variation on the variation in the output quantity using Sobol Sensitivity Analysis.
The input parameters here are the initial concentrations of all the species and the kinetic rate constants of the reactions. Three main output quantities were calculated from the results obtained by running the ensemble of models, and the Sobol sensitivity coefficients were calculated for each of these three quantities. These are:
Integrated βcasein levels
βcasein mRNA levels integrated over the whole time interval (0–48 h) is a direct output of the model which we have used to present the results in the previous section. This quantity is an automatic choice for the sensitivity analysis.
Timedependent switch in the βcasein expression
we consider the relative contributions to the timeintegrated βcasein levels of the JAKdependent and the JAKindependent parts of the pathway. Since we know from the above results that the JAKdependent pathway is operational in the earlier part (0–12 h) while the JAKindependent pathway in the later part (24–48 h) we can calculate the integrated mRNA levels separately for these time intervals and calculate their ratio. When this ratio is near unity, it indicates equal contributions of both these pathways to the net βcasein expression. On the other hand, a significant value of this ratio will indicate that the main contribution comes from JAKdependent pathway and the model does not show a timedependent switch. A similar interpretation can be made when the ratio is small.
Transcription time delay
We get the delay between the times when activated STAT5 dimer reaches the nucleus and when βcasein mRNA transcription starts by getting the difference between the peak times of nuclear STAT5 dimer and cytoplasmic βcasein mRNA. The previous results also showed that the dynamical behavior of the combined pathway is much dependent on the level of ligand NRG stimulation. So, we decided to conduct the global sensitivity analysis at three different NRG stimulation – low (10 nM), medium (20 nM) and high (50 nM).
Timeintegrated βcasein levels
We performed the global sensitivity analysis concerning the total integrated βcasein mRNA levels in the cytoplasm. The results are presented separately for the initial amounts of species and the kinetic rate constants (Fig. 6a and b). The plots show both the first order (S1) and the total effect sensitivities (ST) of the top 10 most sensitive species and parameters of the model. Table 1 below (left column) shows the names of the parameters and the reactions they belong to. From the results, we see the βcasein mRNA expression is sensitive to initial numbers of Her4 and Her2 receptors. It is also sensitive concerning the phosphatases SHP and PPX and JAK which are expected. As for the reactions the sensitivity is high for the transcription parameter (including the Hill coefficient), suggesting that the βcasein expression is sensitive to the mechanism of the transcription reaction. We want to find if the same holds when we consider our main output quantities of interests – timedependent switch and transcription delay. Among other reaction parameters, we see that the rates of JAKindependent activation of STAT also features among the top sensitive parameters. The sensitivity plots also show that the total effect sensitivity S_{T} ≈ S_{1} here for all the sensitive parameters which suggests that there is not much higher order effect of parameters on the absolute βcasein levels.
Timedependent switch in βcasein expression
The next output parameter we considered is the ratio of βcasein mRNA expression during the JAKdependent pathway (0–12 h) to that of JAKindependent pathway (12–48 h). As explained previously, this quantity can be an indicator of whether and to what extent there is a timedependent switch in βcasein mRNA expression. This output quantity was found (Fig. 6c and d) to be particularly sensitive to HER4 heterodimerization and homodimerization reactions (both liganddependent and constitutive) as well as the JAKindependent HER4 activation. The results confirm our hypothesis that both these additions to the literature JAKdependent pathway are sufficient for producing a timedependent switch. Interestingly the transcription reaction parameters which influenced the absolute time integrated mRNA levels to a high degree (from the previous section) did not appear in the list of top parameters for the timedependent switch. These insights highlight the importance of selecting a proper output quantity of interest in Global Sensitivity Analysis as the results can vary depending on the choice of the output parameter. Among the species we see HER4, and HER2 numbers are again significant, along with JAK and cytoplasmic JAK phosphatase SHP2. Here also we see that the total effect sensitivity S_{T} ≈ S_{1} here for all the sensitive parameters which suggests that there is not much higher order effect of parameters on the timedependent switch.
Delay in STAT nuclear translocation and βcasein transcription
The next output parameter of interest is the time delay between STAT5 nuclear translocation and βcasein transcription and transport outside the nucleus. Here we make several interesting observations (Fig. 6e and f). The various transport rates in the model especially the rate of activated STAT5 dimer nuclear import and rate of βcasein mRNA nuclear export were prominent in the list of sensitive parameters. The effect of these transport rates on the pathway activity has been reported before in literature [5, 22] and is expected from the model. The SOCS mediated negative feedback rate also features on the sensitive parameters list. The negative feedback does not directly influence the βcasein mRNA levels but has a more indirect effect. The result shows that in a nonlinear system like HER4JAKSTAT pathway parallel/side reactions can influence the dynamic behavior and can be exploited as therapeutic targets. In contrast to the absolute βcasein sensitivity we see here that S_{T} ≫ S_{1} which suggests that there is a high degree of nonlinear effect of the sensitive parameters to the overall sensitivity of transcription time delay.
In the of the Additional file 1: Figures S1S2 we explore the effect of high and low NRG stimulation on the global sensitivity coefficients. We see some NRG induced changes in the order of parameters and the degree of nonlinearity, but there were no significant changes in the list of most sensitive parameters.
Discussion
Epidermal growth factor family of receptors plays a crucial role in different cancers by activating critical signaling pathways that control cell fate decisions. HER4, a member of this family has received attention due to its several distinctive properties. Some isoforms of this receptor can undergo enzymemediated cleavage reactions resulting in shedding of its ectodomain and cytoplasmic domain leaving an 80 kDa transmembrane domain which can translocate to the nucleus and promote transcription of various genes. On activation by ligand NRG, it can activate the antiproliferative JAK2STAT5A pathway which results in activation, dimerization and nuclear localization of the transcription factor STAT5A which promotes transcription of genes mediating differentiation in particular breast cancer cell lines. This study was motivated by experimental work conducted on mouse HC11 mammary epithelial cell lines stimulated by various ligands including NRG using RTPCR to measure the expression of the important differentiation marker gene βcasein. One of the intriguing results from the study was that the response of these cells to the NRG stimulation (agonistic versus antagonistic) is timedependent. It was observed that whereas NRG suppressed the transcription of βcasein at early time points (012 h), at later time points (2448 h), it promoted transcription. Using ELISA based TF activation assays of STAT5A and GR  two transcription factors required for transcription of βcasein, the same study also found that the activity of these transcription factors persisted even 24 h after NRG stimulation. These observations mirrored similar findings in the literature from ChIPseq studies that showed there is a significant time delay between the STAT5A entry to the nucleus and transcription of the βcasein gene. To obtain a mechanistic understanding of these observations we turned to mathematical modeling. Although the JAKSTAT pathway has been extensively modeled before, the literature models failed to reproduce these experimental findings particularly the timedependent switch in βcasein gene transcription. The failure of the existing models led us to develop a new model for the HER4JAK2STAT5 system that retains the core reactions of the literature models but adds two essential components which have been reported before in the literature but have not been included in the models. These are:

Competitive binding and heterodimerization of HER4 receptor with other members of the ErbB family like HER3.

A HER4 mediated JAKindependent activation of STAT5 which proceeds at a lower rate than the canonical JAKdependent activation.
Including these reactions in our model, we were able to reproduce the experimental findings. By systematically turning these reactions on and off we showed how the signaling dynamics shifted to the pattern seen in the experiments at late time points. We performed extensive parameter sweep studies to understand how different individual components of the model influenced the signaling dynamics. We also conducted global sensitivity analysis to test the robustness of model predictions and to obtain the most sensitive parameters of the model for different output dynamical quantities. These studies show that the competitive HER4 heterodimerization reactions have a profound impact on the sensitivity of the pathway to NRG stimulation at earlier time points. This reaction is a necessary condition for the observed suppression of transcription by NRG.
Along with this competitive heterodimerization, the addition of a slower JAKindependent mechanism of activation of HER4 was sufficient to reproduce the timedependent switch in the transcription of the βcasein gene observed in the experiments. We also found that the various transport rates in the model such as STAT5 dimer nuclear import and βcasein mRNA export influences the time delays associated with transcription. The Global sensitivity analysis results confirm these model findings by showing that the parameters of the above reactions are most sensitive to the corresponding model output such as delay. This study highlights the effect of competitive and parallel reaction pathways on both short and longterm dynamics of receptormediated signaling. Such timedependent alterations in the signaling behavior of these pathways are highly consequential in cancer and are one of the main ways tumor cells develop resistance to targeted inhibitors. Identification of such timedependent changes in signaling dynamics can also help in designing optimal treatment strategies with different dose interval and durations [23, 24]. By obtaining a deeper understanding of the dynamics of such pathways, we will be able to design more efficient drug dosing regimens that can target and exploit the differential dynamics. However, because of the uncertainties inherent in these mechanistic models, there may be alternative mechanisms that might explain the observed data. In those cases, a datadriven clustering based approach [25, 26] which has received attention recently can be used to perform system identification and map dynamical behavior the observed data.
Conclusion
Using a mathematical model of HER4 mediated JAKSTAT pathway that incorporates recently identified reactions of the HER4 receptor, we were able to reproduce the timedependent switching behavior of the differentiation marker βcasein milk gene observed in HC11 mouse mammary epithelial cells. The model demonstrated the necessity of competitive heterodimerization of HER4 receptor and a slower JAKindependent activation of transcription factor STAT5 to reproduce the timedependent switch. Through global sensitivity analysis, we also identified the critical parameters such as HER4 competitive heterodimerization rates and rates of transport of STAT5 transcription factor and βcasein mRNA that profoundly influence the dynamical behavior. This study illustrates the importance of parallel reaction pathways and differences in reaction time scales to produce essential changes in the dynamics of cellular pathways. Such models are valuable means of exploring and designing effective treatment strategies in cancers that can exploit the dynamical behavior of the pathways.
Materials and methods
Pathway description
The HER4JAKSTAT system was modeled as a deterministic reaction network using mass action kinetic equations modeled by ordinary differential equations. The system was modeled as a twocompartment system (for cytoplasm and nucleus) with NRG in the extracellular medium. The species were assumed to be in sufficiently large amounts so that the deterministic approximation applies. The model builds on the canonical JAKSTAT pathway model in the literature [5, 6]. In this pathway, the receptorJAK2 complex is activated by the ligand. This activated receptor complex, in turn, activates cytoplasmic STAT5 which then dimerizes and translocates to the nucleus. The STAT5 dimer in the nucleus initiates transcription of various genes of which the SOCS mRNA which translates to SOCS protein and exerts negative feedback on the circuit by deactivating the activated receptorJAK2 complex is of interest. The other gene of interest for the current HC11 system is, of course, βcasein, which on transcription is transported outside the nucleus. This cytoplasmic βcasein is reported in the experiments.
To explain the experimental findings which showed that the canonical pathway is inadequate to explain things like a timedependent switch in βcasein mRNA levels we incorporated two main modifications of the canonical pathway based on recent literature findings

a)
HER4 can form both homo and heterodimers with other members of the ErbB family on activation by Neuregulin. These reactions compete with the dimerization and activation reactions of HER4JAK2 complex. Such competition for ligand has been shown to induce an inverse ligand dependence (signaling activity decreases with increasing ligand stimulation) in these pathways which is also observed in the experiments. These additional homodimerization and heterodimerization steps were modeled using three reactions.

b)
The requirement of JAK2 for tyrosine phosphorylation of STAT5 was only demonstrated in HBEGF and Prolactin stimulated pathways [3]. There are other shreds of evidence in the literature that suggest STAT5A can directly interact with HER4 s80 when stimulated by HRG/NRG. Hence in the present work, we incorporated this possibility by introducing a JAK2 independent pathway which allows a direct interaction and activation of the s80/4ICD domain of HER4 with STAT5a. The most important effect of adding this pathway is that being JAK independent, this complex will not be negatively regulated by SOCS to the extent when JAK is present [17] which can increase the overall βcasein gene expression in a timedependent manner. The rate constants of the activation reactions for this JAK2 independent pathway were taken as smaller than that of the JAK2dependent pathway since we expect the active kinase domain of JAK2 to catalyze the phosphorylation of STAT5.
The rate constants were taken from published models [5, 6, 27] for the reactions which were common to the current model. For the new reactions, we estimated the rate constants by starting with similar values as the related reactions for which the parameters are known, and then doing sensitivity analysis and comparing with experimental results. Bounds for the rate constants were also confirmed independently using values/estimates of the diffusion coefficient. The model was solved using COPASI [28]. A full description of the model along with the initial expressions and rate constants are provided in the SBML (http://sbml.org/Documents/Specifications) format in Additional file 2. We used Latin Hypercube sampling [29] to sample the parameter space and to create an ensemble of models. Then for each member of the ensemble, we calculate and plot the timeintegrated value of βcasein mRNA and average them over the ensemble. This task ensures that the uncertainties in the model parameters are reflected in the predictions to a certain degree and that the model is robust to small perturbations in the model parameters.
Global sensitivity analysis
Real biological pathways like HER4JAK2STAT5 are nonlinear and various parameters such as initial species expression, kinetic rate constants, etc. can undergo large deviations. Sensitivity analysis techniques are useful to understand how perturbations in these parameters influence the model outputs. Because these systems have a high degree of nonlinearity, simple local sensitivity analysis where each parameter is altered one at a time keeping the others fixed is not sufficient [20]. Hence, we use a global sensitivity analysis for this model – more specifically the Sobol Sensitivity Analysis which is based on analysis of variances in the model parameters and output.
Sobol sensitivity analysis
Sobol sensitivity analysis is a variancebased technique and is uniquely well suited for complex nonlinear systems of moderate size. Here we give a very brief description of the method that is relevant to the sensitivity analysis done in this paper. A detailed description can be found in references [21, 30]. We use many of the notations from [30] below.
Suppose we have a model with k parameters X_{1}, X_{2}, … , X_{k} which are assumed to be independent random variables. The model output Y is related to these parameters through the relation Y = f(X_{1}, X_{2}, … , X_{k}) where f is a general nonlinear function. In a variancebased sensitivity analysis, we want to understand how variances in the individual and combination of parameters X_{1}, X_{2}, … , X_{k} factor into the variances in Y. To determine this, we can first fix a parameter X_{i} to a value (say V_{i}) and then determine the model output averaged over all remaining parameters X_{j}, i ≠ j (which is denoted with a condensed notation \( {\overline{X}}_{\sim i} \)). This average is \( {E}_{{\overline{X}}_{\sim i}}\left(Y{X}_i={V}_i\right) \) which will be different for different V_{j}. The variance of this average over all possible V_{i} which is \( {V}_{X_i}\left({E}_{{\overline{X}}_{\sim i}}\left(Y{X}_i={V}_i\right)\right) \) will give us the net first order effect of variation in X_{i} on the variation in Y. The first order sensitivity S_{i} associated with parameter X_{i} is defined as:
In the above V(Y) is the overall variance in Y.
The other sensitivity parameter of interest is the total effect sensitivity \( {S}_{T_i} \) which represents the first and all higher order effects of the parameter X_{i} on the model output. To determine this, we can start with determining the first order effect of all parameters except X_{i} which is denoted by \( {\overline{X}}_{\sim i} \). Again, we first find the value of Y averaged over all X_{i} keeping all other parameters \( {\overline{X}}_{\sim i} \) fixed which is \( {E}_{X_i}\left(Y{\overline{X}}_{\sim i}={\overline{V}}_{\sim i}\right) \). We then find the variance of this quantity over all possible \( {\overline{X}}_{\sim i} \) or \( {V}_{{\overline{X}}_{\sim i}}\left(\ {E}_{X_i}\left(Y{\overline{X}}_{\sim i}={\overline{V}}_{\sim i}\right)\right) \). Then \( V(Y){V}_{{\overline{X}}_{\sim i}}\left(\ {E}_{X_i}\left(Y{\overline{X}}_{\sim i}={\overline{V}}_{\sim i}\right)\right) \) must represent the contribution of all terms where X_{i} appears. Dividing this by V(Y) we get the total effect sensitivity:
This variancebased sensitivity analysis framework is based on a functional decomposition scheme where a square integrable function Y = f(X_{1}, X_{2}, … , X_{k}) defined over Ω, the kdimensional unit hypercube, can be expressed as follows:
The normalization condition here is
In the above w = i_{1}, … , i_{s}. Using this the various terms are calculated as
Taking the variances of both sides of these equations give us
All these variances are linked by
Dividing the above by V(Y), we get
Since the sensitivity coefficients above are multidimensional integrals, the standard way of computing them is by using Monte Carlo type sampling. In Monte Carlo sampling an integral I[f] = ∫ f(x)dx is computed by generating a sequence of uniformly distributed random numbers and computing their expectation \( {I}_N\left[f\right]=\frac{1}{N}{\sum}_{i=1}^Nf\left({x}_i\right) \). In this case, we have a kdimensional function Y = f(X_{1}, X_{2}, … , X_{k}). Hence, we need to sample N times for each of these k parameters. These samples can be represented by a N × k matrix. For the calculation of the above sensitivities, the standard procedure is to start with two independent N × k sampling matrices \( \overline{\overline{A}} \) and \( \overline{\overline{B}} \). We can compute matrices \( {\overline{\overline{A}}}_B^i \) which is obtained by taking \( \overline{\overline{A}} \) and replacing the ith column (for parameter X_{i}) by the corresponding column from \( \overline{\overline{B}} \). Similarly, we can define \( {\overline{\overline{B}}}_A^i \). It can be shown that [30] the variances in the equations for S_{i} and \( {S}_{T_i} \) can be estimated using
The above equations form the basis of the computation of the sensitivity coefficients for the model.
Quasirandom sequences
Monte Carlo method using a sequence of random or pseudorandom numbers is extensively used to compute multidimensional integrals like above. The Central Limit Theorem of probability shows that the error in the Monte Carlo estimate is equal to the product of the standard deviation of the function and \( {N}^{\frac{1}{2}} \) where N is the number of samples. Hence the convergence of this method is \( O\left({N}^{\frac{1}{2}}\right) \) which can be very slow [31]. This method of sampling using pseudorandom numbers also suffers from a related problem of clumping where the sample points often tend to clump together and leaves empty spaces in between which is magnified in higher dimensions. One alternative to obtaining a more uniform distribution of points is by using a stratified sampling method like Latin Hypercube Sampling which divides the intervals into equally spaced points. However, this only works when the dimensionality is low. For integrations in higher dimensions, an alternative sampling technique is applied called quasirandom sampling. A quantitative measure of uniformity of a sequence is a factor termed “discrepancy”. Lower the discrepancy, more uniform is the sequence. Suppose we have a sequence of N points {x_{N}} in the k dimensional unit cube I^{k}. For any subset J ∈ I^{k} one can define the error in Monte Carlo estimate of the volume of J as [31]:
The discrepancy is then defined as some norm of R_{N}(J). Formally, if J is restricted to rectangular set and E is all possible such sets then the discrepancy D_{N} is defined as
The KoksmaHlawka inequality provides an upper bound for a Monte Carlo integration error as a product of the variance of the function and the discrepancy of the sequence. Hence, if one can generate a sequence of points in such a way as to minimize the discrepancy, then it can be used to obtain an improved estimate of the integral. Such sequences are called quasirandom sequences. These are not random at all but are generated deterministically to minimize discrepancy. Since they are not random, quasirandom sequences are more limited in scope than pseudorandom numbers. However, it can be shown that for integration the convergence rate of quasirandom sequences is O(N^{−1}(logN)^{k}) where k is any number which is considerably faster than the \( O\left({N}^{\frac{1}{2}}\right) \) convergence of the standard MonteCarlo method using pseudorandom sequences. There are various techniques for determining quasi random sequences. We use Sobol sequences [32] using a method suggested by Saltelli [30]. The software package SALib [33] was used for the computation of the Sobol coefficients along with custom python scripts and matplotlib [34] for plotting.
References
 1.
Citri A, Yarden Y. EGFERBB signalling: towards the systems level. Nat Rev Mol Cell Biol. 2006;7(7):505–16.
 2.
Thor AD, Edgerton SM, Jones FE. Subcellular localization of the HER4 intracellular domain, 4ICD, identifies distinct prognostic outcomes for breast cancer patients. Am J Pathol. 2009;175(5):1802–9.
 3.
Williams CC, et al. The ERBB4/HER4 receptor tyrosine kinase regulates gene expression by functioning as a STAT5A nuclear chaperone. J Cell Biol. 2004;167(3):469–78.
 4.
Chen WW, et al. Inputoutput behavior of ErbB signaling pathways as revealed by a mass action model trained against dynamic data. Mol Syst Biol. 2009;5:239.
 5.
Swameye I, et al. Identification of nucleocytoplasmic cycling as a remote sensor in cellular signaling by databased modeling. Proc Natl Acad Sci. 2003;100(3):1028–33.
 6.
Yamada S, et al. Control mechanism of JAK/STAT signal transduction pathway. FEBS Lett. 2003;534(1–3):190–6.
 7.
Greulich H, et al. Oncogenic transformation by inhibitorsensitive andresistant EGFR mutants. PLoS Med. 2005;2(11):1167.
 8.
Telesco ST. Multiscale modeling of the ErbB receptor tyrosine kinase signaling network through theory and experiment: PhD Thesis, University of Pennsylvania; Publicly Accessible Penn Dissertations. 971;2011. https://repository.upenn.edu/edissertations/971.
 9.
MuraokaCook RS, et al. The intracellular domain of ErbB4 induces differentiation of mammary epithelial cells. Mol Biol Cell. 2006;17(9):4118–29.
 10.
MuraokaCook RS, et al. Prolactin and ErbB4/HER4 signaling interact via Janus kinase 2 to induce mammary epithelial cell gene expression differentiation. Mol Endocrinol. 2008;22(10):2307–21.
 11.
Koibuchi H, et al. Grand canonical Monte Carlo simulations of elastic membranes with fluidity. Phys Lett A. 2003;319:44.
 12.
Riese DJ. Ligandbased receptor tyrosine kinase partial agonists: new paradigm for cancer drug discovery? Expert Opin Drug Discovery. 2011;6(2):185–93.
 13.
Gilmore JL, Riese DJ 2nd. secErbB426/549 antagonizes ligandinduced ErbB4 tyrosine phosphorylation. Oncol Res. 2004;14(11–12):589–602.
 14.
Gilmore JL, Scott JA, Bouizar Z, Robling A, Pitfield SE, Riese DJ, Foley J. AmphiregulinEGFR signaling regulates PTHrP gene expression in breast cancer cells. Breast Cancer Res Treat. 2008;110(3):493–505.
 15.
Willmarth NE, Baillo A, Dziubinski ML, Wilson K, Riese DJ 2nd, Ethier SP. Altered EGFR localization and degradation in human breast cancer cells with an amphiregulin/EGFR autocrine loop. Cell Signal. 2009;21(2):212–9.
 16.
Clark DE, et al. ERBB4/HER4 potentiates STAT5A transcriptional activity by regulating novel STAT5A serine phosphorylation events. J Biol Chem. 2005;280(25):24175–80.
 17.
Hennighausen L, Robinson GW. Interpretation of cytokine signaling through the transcription factors STAT5A and STAT5B. Genes Dev. 2008;22(6):711–21.
 18.
Kabotyanski EB, et al. Integration of prolactin and glucocorticoid signaling at the betacasein promoter and enhancer by ordered recruitment of specific transcription factors and chromatin modifiers. Mol Endocrinol. 2006;20(10):2355–68.
 19.
Kabotyanski EB, et al. Lactogenic hormonal induction of long distance interactions between betacasein gene regulatory elements. J Biol Chem. 2009;284(34):22815–24.
 20.
Zi Z. Sensitivity analysis approaches applied to systems biology models. IET Syst Biol. 2011;5(6):336–46.
 21.
Sobol IM. Sensitivity estimates for nonlinear mathematical models. Math Modeling Comput Experiment. 1993;1(4):407–14.
 22.
Gambin A, et al. Computational models of the JAK1/2STAT1 signaling. Jakstat. 2013;2(3):e24672.
 23.
Kolch W, Halasz M, Granovskaya M, Kholodenko BN. The dynamic control of signal transduction networks in cancer cells. Nat Rev Cancer. 2015;15(9):515.
 24.
Lee MJ, Ye AS, Gardino AK, Heijink AM, Sorger PK, MacBeath G, Yaffe MB. Sequential application of anticancer drugs enhances cell death by rewiring apoptotic signaling networks. Cell. 2012;149(4):780–94.
 25.
Kwon J, Jayaraman A, Lee D. Identification of a timevarying intracellular signaling model through data clustering and parameter selection: application to NFkB signaling pathway induced by LPS in the presence of BFA. IET Syst Biol. 2019; in press. https://doi.org/10.1049/ietsyb.2018.5079.
 26.
Mangan NM, Askham T, Brunton SL, Kutz JN, Proctor JL. Model selection for hybrid dynamical systems via sparse regression. Proc R Soc A. 2019;475(2223):20180534.
 27.
Schoeberl B, et al. Therapeutically targeting ErbB3: a key node in ligandinduced activation of the ErbB receptorPI3K axis. Sci Signal. 2009;2(77):ra31.
 28.
Hoops S, et al. COPASI—a complex pathway simulator. Bioinformatics. 2006;22(24):3067–74.
 29.
McKay MD, Beckman RJ, Conover WJ. Comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics. 1979;21(2):239–45.
 30.
Saltelli A, et al. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. Comput Phys Commun. 2010;181(2):259–70.
 31.
Caflisch RE. Monte carlo and quasimonte carlo methods. Acta Numerica. 1998;7:1–49.
 32.
Sobol IM. Uniformly distributed sequences with an additional uniform property. USSR Comput Math Math Phys. 1976;16(5):236–42.
 33.
Herman J, Usher W. SALib: an opensource Python library for sensitivity analysis. J Open Source Softw. 2017;2(9):90–95.
 34.
Hunter JD. Matplotlib: a 2D graphics environment. Comput Sci Eng. 2007;9(3):90–5.
Acknowledgments
The authors thank the members of the CHIC consortium and PSOC members at Penn for helpful discussions.
Funding
The authors acknowledge funding from the EU through grant FP7ICT20119600841 and from National Institutes of Health through grants U54 CA193417, and U01 CA227550. These funding bodies played no roles in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Availability of data and materials
All datasets and models are after available after publication upon request to the corresponding author.
Author information
Affiliations
Contributions
RR and AG developed the molecular model. All authors participated in data analysis and in writing the manuscript. All authors have read and approved the final manuscript, and ensure that this is the case.
Corresponding author
Correspondence to Ravi Radhakrishnan.
Ethics declarations
Ethics approval and consent to participate
None required.
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:
Collection of all the supplementary figures showing the results of sensitivity analysis and parameter sweep studies (PDF 2789 kb)
Additional file 2:
The detailed HER4JAKSTAT model with the reactions, initial expression of the proteins and the reaction rate constants in Systems Biology Markup Language (SBML) format. (XML 58 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
Received
Accepted
Published
DOI
Keywords
 RTK signaling
 JAKSTAT pathway
 Timedependent switch
 Global sensitivity