Sampling-based Bayesian approaches reveal the importance of quasi-bistable behavior in cellular decision processes on the example of the MAPK signaling pathway in PC-12 cell lines

Abstract

Background: Positive and negative feedback loops are ubiquitous motifs in biochemical signaling pathways. The mitogen-activated protein kinase (MAPK) pathway module is part of many distinct signaling networks and comprises several of these motifs, whose functioning depends on the cell line at hand and on the particular context. The maintainance of specificity of the response of the MAPK module to distinct stimuli has become a key paradigm especially in PC-12 cells, where the same module leads to different cell fates, depending on the stimulating growth factor. This cell fate is regulated by differences in the ERK (MAPK) activation profile, which shows a transient response upon stimulation with EGF, while the response is sustained in case of NGF. This behavior was explained by different effective network topologies. It is widely believed that this sustained response requires a bistable system.

Results: In this study we present a sampling-based Bayesian model analysis on a dataset, in which PC-12 cells have been stimulated with different growth factors. This is combined with novel analysis methods to investigate the role of feedback interconnections to shape ERK response. Results strongly suggest that, besides bistability, an additional effect called quasi-bistability can contribute to explain the observed responses of the system to different stimuli. Quasi-bistability is the ability of a monostable system to maintain two distinct states over a long time period upon a transient signal, which is also related to positive feedback, but cannot be detected by standard steady state analysis methods.

Conclusions: Although applied on a specific example, our framework is generic enough to be also relevant for other regulatory network modeling studies that comprise positive feedback to explain cellular decision making processes. Overall, this study advices to focus not only on steady states, but also to take transient behavior into account in the analysis.

Keywords: Cellular decision making; MAPK signaling pathway; Quasi-bistability.

Fig. 1

Activities of Raf, MEK and ERK after stimulation. Scaled activities of Raf, MEK and ERK measured by polychromatic flow cytometry (by visual inspection from Fig. S1b in [12])

Fig. 2

Data from modular response analysis. Table. Means and standard deviations of the global response coefficients extracted from the silencing experiments via modular response analysis. These were calculated from replicates in Table S1d in [12]. Figure. Illustration of respective changes in protein concentrations in response to silencing relative to the control experiments (without silencing)

Fig. 3

Model structure of the MAPK module. a Reaction scheme of the MAPK module. Upon addition of growth factors, Raf, MEK and ERK are successively activated in a phosphorylation cascade. Different feedback topologies are assumed to shape context dependent ERK response: Effective negative feedback from ERK to Raf upon EGF stimulation (dotted line from ppERK to dephosphorylation of pRaf), and positive feedback in case of NGF stimulation (dashed line from ppERK to phosphorylation of Raf). b Differential equation model of the MAPK cascade. Bold parameters are the unknown constants, collected in the parameter vector θ, while gray parameters define the specific experimental condition for the simulation

Fig. 4

Calibrated model using a Bayesian approach. Dynamic responses of pRaf, ppMEK and ppERK after stimulation with EGF a and NGF b. Shown are the values acquired from flow cytometry experiments (Fig. 1) in comparison to the respective PPDs predicted by the model. Data have been normalized to t=5 min. c Comparison of data for the global response coefficients (GRC) extracted from the siRNA perturbation experiments (Fig. 2), and respective simulated distributions, here for clarity represented with the first and second moment. Data are taken from [12]

Fig. 5

Model validation. a Dose-response profiles of ERK activation were mimicked by simulating the model with increasing input strength parameter k _u for stimulation with EGF (left) and NGF (right). The system shows a unimodal and ultrasensitively increasing ppERK concentration after stimulation with EGF (t=5 min after stimulation) and a bimodal distribution when stimulated with NGF exceeding a threshold concentration (t=60 min after stimulation) (compare data in [12], Subfigs 2c and d). b Inhibition of MEK (left) results in the loss of sustained Raf activation upon stimulation with NGF (gray dashed PPDs) compared to the control case (blue continuous PPDs). Inhibition of PKC via Gö7874 (right) causes the loss of sustained ERK activation upon NGF stimulation (data from [12], Fig. 4a). This was simulated by switching off the feedback connection. c Irreveversibility in MAPK activation upon NGF stimulation was investigated via mimicking treatment of the cell culture with neutralizing antibodies (left) and TrkA inhibitors (right) (compare data in [12], Subfigs 3a and c)

Fig. 6

Steady state analysis using the circuit-breaking algorithm. The CBA is used for an efficient calculation of the steady states of the system for the MCMC parameter samples and subsequent automatic classification into mono- and bistable systems (a-c). d Result of this classification analysis. Depicted is also the distribution of the second stable steady state ${\overline{z}}_3\ne 0$ in case of a bistable system, which corresponds to the concentration of active ERK normalized to t=5 min (see Additional file 2)

Fig. 7

Simulation-based analysis of the long-term behavior of the cell population. a Trajectories are automatically classified into bistable, quasi-bistable and monostable systems, as described in the text. For stimulation with NGF in the control experiment 90% show a quasi-bistable behavior, while only 10% are really bistable. b Distribution of the second stable steady state distinct from zero of the bistable trajectories. c For the choice of threshold parameters that were used for the classification scheme we evaluated the 0% and 100% percentiles of ERK trajectories in the simulation of the NGF control experiment

Fig. 8

Quasi-bistability phenomenon. The CBA is used for the investigation of the quasi-bistability phenomenon, in which the system, despite being monostable, shows a very prolonged sustained response. The first column shows the time course of normalized ppERK for a representative parameter sample from class 2 with switching time at t _switch≈440 min. Columns 2,3 and 4 show the circuit-characteristic c(κ,u(t)), along with the actual normalized state ppERK(t) for 12 different time points. After a fast transient dynamic (a1) the circuit-characteristic has three zeros (A2-B1), which disappear at a later time point, here t=102 min (b2), via a saddle-node bifurcation. After 60 min the input is almost zero and the vector field and therefore the circuit-characteristic changes only slowly. The system state has almost approached the higher fixed point. b1-c3 are eyeglass views on the dynamics near this second fixed point. These plots show that, even if the fixed point has disappeared, the system trajectory moves very slowly through the state space for a rather long time, since $\dot{x}$ is still small. Only after about 440 min the system has overcome this slow region of the state space, and from here on rapidly moves towards its globally asymptotically stable steady state $\overline{x}=0$

Fig. 9

Combination of two delay mechanisms in quasi-bistable systems. a Scheme of a bifurcation diagram for a quasi-bistable system. The system is monostable for u=0 and has a saddle-node bifurcation u ^SNB close to u=0, where it becomes bistable. A sufficiently strong transient signal u(t) pushes the system state into the basin of attraction of the higher stable steady state (1). As long as the change in u(t) is not slow compared to the dynamics of the system, the system cannot be considered in quasi-steady state, and we observe a transient dynamics (2). When u(t) is almost back to 0, two delay effects lead to quasi-bistable behavior (3). First, the system remains in the upper stable steady state as long as u(t) is still above the saddle-node bifurcation. Second, for u(t)<u ^SNB the acceleration remains very small in this region of the state space. b Absolute value of the vector field along the model trajectory for the same model parameters that have been used in Fig. 8. c Two respresentative bifurcation diagrams for a quasi-bistable system belonging to class 2 of the classification scheme, and a bistable system belonging to class 1