Research ArticleCOMPUTATIONAL BIOLOGY

Automatic Generation of Predictive Dynamic Models Reveals Nuclear Phosphorylation as the Key Msn2 Control Mechanism

See allHide authors and affiliations

Science Signaling  28 May 2013:
Vol. 6, Issue 277, pp. ra41
DOI: 10.1126/scisignal.2003621

Abstract

Predictive dynamical models are critical for the analysis of complex biological systems. However, methods to systematically develop and discriminate among systems biology models are still lacking. We describe a computational method that incorporates all hypothetical mechanisms about the architecture of a biological system into a single model and automatically generates a set of simpler models compatible with observational data. As a proof of principle, we analyzed the dynamic control of the transcription factor Msn2 in Saccharomyces cerevisiae, specifically the short-term mechanisms mediating the cells’ recovery after release from starvation stress. Our method determined that 12 of 192 possible models were compatible with available Msn2 localization data. Iterations between model predictions and rationally designed phosphoproteomics and imaging experiments identified a single-circuit topology with a relative probability of 99% among the 192 models. Model analysis revealed that the coupling of dynamic phenomena in Msn2 phosphorylation and transport could lead to efficient stress response signaling by establishing a rate-of-change sensor. Similar principles could apply to mammalian stress response pathways. Systematic construction of dynamic models may yield detailed insight into nonobvious molecular mechanisms.

Introduction

The complexity of biological networks and processes makes mathematical models useful as compact representations of the available data, and as highly structured, mechanistic representations of biological hypotheses that can be used to generate quantitative testable predictions of the network behavior (1, 2). When the dynamics of a system is of interest, such as in cell signaling, and when stochastic noise due to low molecular copy numbers can be neglected, these models typically consist of systems of ordinary differential equations (ODEs) to capture the kinetics of proteins, mRNAs, and small molecules. In these models, there is a distinction between state variables, such as time-varying concentrations of components, and parameters, which usually represent kinetics constants of biochemical reactions; the evolution of state variables in time is dictated by both the pathway topology and the parameter values. Modelers often face the harsh reality of conflicting biochemical hypotheses, or even a complete lack of hypotheses, about multiple mechanistic details. Such knowledge gaps result from current—and often fundamental—limitations in generating comprehensive quantitative experimental data. They lead to two classes of ambiguities of mathematical models in systems biology, namely, ambiguities in model topologies (for example, pathway structures and kinetic rate laws) and uncertainties in parameter values (for example, kinetic constants for processes such as association and dissociation of molecules).

Uncertainties in model structure make it necessary to evaluate more than one candidate model topology, where each topology represents a different biological hypothesis. Model discrimination then consists of evaluating alternative topologies—or candidate models—for their ability to describe experimental observations, with the ultimate aim of identifying the “real” system structure (35). The set of candidate models may, for example, take the form of a core model, which incorporates the consensus understanding of a pathway, and a set of additional hypotheses to extend this core to an ensemble of hypothetical models (3). Such approaches have been successful in identifying detailed mechanisms of epidermal growth factor (EGF) (6) and mammalian target of rapamycin (mTOR) (7) signaling. These studies, however, compared small numbers (four and three, respectively) of hypotheses. In general, the ensemble size grows exponentially with the number of hypotheses (for example, 20 hypotheses yield 220 ≈ 106 candidate models), such that one cannot simply enumerate and evaluate all alternatives individually. Systematic methods for model generation that deal with this complexity of biological hypothesis spaces are currently lacking.

The computational problems are particularly challenging because of our currently limited quantitative knowledge about many model parameters. Even if a model is structurally identifiable, meaning that a unique assignment of parameter values is possible in principle (8), it may be unidentifiable with the available data. This can translate into ambiguous predictions of a system’s quantitative or even qualitative behavior. Therefore, a single set of parameter values (referred to as a parameter point) may be insufficient to predict the behavior of nonmeasured molecular species and to select between model topologies. Consequently, the evaluation of a candidate model also requires an exploration of the model’s parameter space, and the associated computational effort increases with the dimensionality of this space. With efficient exploration methods, however, one can use Bayesian approaches (6, 9, 10) to take into account that different model parameters and topologies can be compatible with the data. With the Bayesian approach, the probability that various hypotheses are correct is updated by sequentially incorporating experimental observations. This helps in ranking alternative topologies (4, 11) and in setting bounds for the predictions of individual models (12, 13).

Here, we propose a method for automatic model generation termed “topological filtering” that only needs to evaluate a subset of all possible topologies in the Bayesian framework but still provides insights about the mechanisms of the underlying biochemical system. After introducing the method with a small example network, we applied it to large previously published models of EGF signaling (6), reproducing and extending the published results. Finally, we combined iterations of experimental and theoretical approaches to study the short-term dynamics of the transcription factor Msn2, which is a key regulator of the general stress response in Saccharomyces cerevisiae (1416). Nuclear localization, and hence Msn2 activity, is under the tight control of environmental conditions, and it is triggered by several stresses, including osmotic stress, heat shock, and carbon or nitrogen limitation. Msn2 localization is regulated by direct phosphorylation through cyclic adenosine monophosphate (cAMP)–dependent protein kinase A (PKA), which targets phosphorylation sites within both the nuclear localization signal (NLS) and the nuclear export signal (NES). Under favorable conditions in the presence of glucose, the second messenger cAMP is produced from adenosine triphosphate, and it activates the catalytic subunits of PKA (encoded by Tpk1 to Tpk3) by binding to PKA’s inhibitory subunit Bcy1. Addition of glucose to starved cells triggers a rapid, but transient, accumulation of cAMP because of strong PKA-dependent feedback regulation of cAMP turnover, leading to translocation of Msn2 from the nucleus to the cytoplasm. Altered dynamics of Msn2 translocation induced by different stresses substantially affect Msn2-dependent gene expression (17), and it is still unclear how Msn2 can serve as an integration point for such a wide range of stresses (16). Our analysis suggests rapid, predominantly nuclear phosphorylation as the most important mechanism for glucose-induced deactivation of Msn2, and that the network controlling the translocation of Msn2 senses the rate of change in the cAMP input, and not other features, such as the total input. The kinetic principles and sensor function may apply to the control of other eukaryotic transcription factors or to any signaling molecule that exists in at least four different states with respect to activity, modification, or localization.

Results

The method of topological filtering

Topological filtering starts from the most complex model that includes all of the known or hypothesized interactions. This defines a single model topology, which we term the original model (Fig. 1A). It then explores (samples) regions of parameter space where the model is consistent with an experimentally observed behavior, given the uncertainty in the experimental data. The method uses the resulting set of parameters to evaluate the model topology with a Bayesian approach and to identify parameters (that is, kinetic constants) as candidates to be eliminated from the model. The elimination of a parameter that is proportional to the reaction rate will inactivate the corresponding reaction. This yields one or several reduced models that are then subject to further cycles of exploration and reduction, thus generating a tree of topologies, in which every topology is evaluated until the simplest possible topology in each branch is found. Topologies that can be no further reduced capture key mechanisms of the underlying system, provided that these mechanisms were included in the original model. This topological filtering can substantially reduce the number of models that need to be evaluated computationally.

Fig. 1 Method flow chart and application to a small example system.

(A) Flow chart for the computational method starting from a model that incorporates all hypothetical mechanisms. The parameter space of the original model is explored by the method described in (20) to automatically generate a set of new, reduced models and to compute Bayes factors. The previous steps are then repeated for each candidate model until no new models are generated. (B) Example model based on (18) consisting of three chemical species (A, B, and C) and three reactions (solid arrows). The original model is constructed by extension with species D and with two hypothetical reactions (dashed lines) that represent complex formation of A and C to yield D. (C) Ranking of the ensemble of eight candidate models according to their model probability given the synthetic data for species A and C generated by simulation of model 4.

To describe the method’s main steps more specifically, we used a small example system that captures the enzyme-catalyzed conversion of a substrate A to a product C in a two-step mechanism [Fig. 1B, see also (18) and text S1]. Similar small systems often occur in larger biochemical networks. We generated a set of synthetic data for A and C from this model, and we extended the model by two reactions, such that the original model takes the form shown in Fig. 1B (see Materials and Methods for details). In the following example, we considered all reactions in the original model as hypothetical.

Our method considers dynamic (ODE) models, with the ith model denoted by Mi. By setting individual parameters of a model Mi to zero, a part of the model’s structure is removed, altering the model’s topology. For example, a reaction based on mass action kinetics vanishes if the corresponding rate parameter is set equal to zero; for other rate laws, more than one parameter might have to be modified. Thus, by inclusion or elimination of each of the d parameters in model Mi, we can form (at most) 2d unique candidate models.

For the original model (and subsequently for each candidate model), topological filtering identifies those regions of the parameter space for which the model’s predictions are compatible with the observational data; we termed these regions “viable regions.” To separate viable from inviable parameter points, we derived a formal viability condition (see Materials and Methods and text S2 for details). Assuming uncorrelated and normally distributed measurement noise, the viability condition depends only on the number of measurements and on the measurement noise [without these restrictive assumptions about the measurement of noise, one can use bootstrapping to compute the viability condition (19)]. With the viability condition, we used the method described in (20) to explore the parameter space, which results in a set of viable parameter points. Next, the method uses these points for model reduction as described above until no new candidate models are found (see Fig. 1A and Materials and Methods for details).

Topological filtering, when applied to the model of Fig. 1B, identified the original model M1 and seven simpler candidate models as compatible with the synthetic data (Fig. 1C). Model M4, which we used to generate the synthetic data, is one of these models. Note that parameters k1, k3, and k4 appear in all viable models, indicating that they (and the corresponding reactions) are essential to explain the data. However, parameters k2, k5, k6, and k7 are only included in half of the candidate models. That is, although k2 is actually present in the model that generated the data, there are candidate models that do not involve this kinetic parameter. However, note that the elimination of the half-maximal rate parameter k2 does not eliminate reaction r1; instead, its rate becomes constant and equal to k1 (see text S1).

To compare and rank the candidate models, we followed a Bayesian inference approach (21), where previous knowledge is encoded by prior probability distributions. By incorporating the data Y, one can use Bayes factors (22) to compute the posterior probability p(Mi|Y), which provides a measure of the “plausibility” of the model Mi, given the data Y (see Materials and Methods for details). We obtained the highest probability (around 39%) for the model that actually generated the data, M4(567), where the superscript denotes that parameters k5, k6, and k7 are deleted in this model (Fig. 1C). However, this model probability does not dominate those of all other models. Hence, to identify the correct topology even for this small system, we have to iterate between computational analyses and carefully designed experimental data generation until the support for one model is satisfactory.

Automatic topology identification for the extracellular signal–regulated kinase pathway in PC12 cells

To investigate the applicability of our method to larger models, we used it for model generation and discrimination for the extracellular signal–regulated kinase (ERK) pathway in PC12 cells. A previous, combined computational and experimental study (6) manually defined a set of four hypothetical mathematical ODE models for the pathway, M1–4. These models incorporate 27 to 30 states (that is, time-varying component concentrations), 51 to 57 parameters, and 27 to 31 reactions. All the mechanisms of M1 are included in M2–4 as well. However, the mechanisms for ERK activation through C3G (Crk SH3 domain guanine nucleotide exchanger) and Rap1 are only included in M2,4, and mechanisms for EGF receptor (EGFR) degradation are only included in M3,4. The experimental data set comprises 168 data points (including replicates) for phosphorylation-mediated ERK activation for 11 distinct conditions and perturbations (6). The predictive performance of the models was compared using Bayes factors, given the available data, and the authors concluded that M2 best represents the pathway (6). This result indicates that C3G is important in ERK activation, whereas the effect of EGFR degradation is negligible.

Because the reactions in each of the models M1–3 constitute a subset of the reactions in the original model M4, we investigated if our method, starting from M4, automatically identified viable submodels, given the experimental observations. Such a procedure does not require the manual definition of M1–3 and may reveal if any of the other model reactions are not supported by experiments. We focused on reactions that are not part of the linear phosphorylation cascade between Ras and ERK (Fig. 2A). The first five reactions represented feedback mechanisms: ERK-mediated removal of SOS (r1), activation of Rap1 by EPAC (r2), guanosine triphosphatase–activating protein (GAP)–mediated deactivation of Rap1 (r3), B-Raf–mediated activation of mitogen-activated or extracellular signal–regulated protein kinase kinase (MEK) (r4), and Ras-mediated phosphorylation of B-Raf (r5). We also included the degradation of EGFR (r6) and the mechanisms for C3G-mediated activation of Rap1 [EGFR-mediated activation of C3G (r7), deactivation of C3G (r8), and C3G-mediated activation of Rap1 (r9)]. Note that a reaction can be associated with more than one parameter, for example, when the reaction rate is modeled with a Michaelis-Menten–type rate law that requires a maximal rate constant and an affinity constant.

Fig. 2 Analysis of the ERK pathway.

(A) Schematic representation of the ERK pathway activated by EGFR, with emphasis on the reactions investigated by our method. The mitogen-activated protein kinase (MAPK) phosphorylation cascade is represented by a single module. The factors stimulating the ERK pathway (ovals) are EGF, PKA agonist (PKAA), cAMP phosphodiesterase antagonist (cilostamide), and EPAC agonist (EPACA). The proteins catalyzing reactions r2, r3, and r5 are specified within brackets. The readout of the pathway is the ratio between phosphorylated ERK (ERKPP) and the total concentration of ERK. Reactions r1 to r9 are targeted by our method (dashed lines). Reactions r6 to r9 can be removed (blue) without compromising the descriptive power of the model. (B and C) Projections of the viable parameter points identified for the original model. (B) For reaction r1 (ERK-mediated degradation of SOS), the maximal rate constant Vmax is plotted against the Michaelis-Menten constant KM. (C) The Vmax for reaction r4 (B-RafPP–mediated activation of MEK, within the MAPK module) on the x axis and the mass action constant k for reaction r6 (degradation of EGFR) on the y axis.

We first explored the viable parameter space of the original model as a basis for our automatic method. Even simple visual inspection of the viable parameter points may reveal candidate parameters for model reduction. For instance, the viable parameter points projected to the maximal rate parameter Vmax for reaction r1 (ERK-mediated degradation of SOS) on the x axis and the affinity constant KM for the same reaction on the y axis indicated that, for r1, Vmax is essential, whereas KM might be expendable if Vmax is sufficiently small (Fig. 2B). Similarly, when comparing the maximal rate constant for r4 and the degradation rate for r6 (Fig. 2C), because all viable parameter points have small values for reaction r6, but not for reaction r4, only r6 is a candidate for reduction.

From the 29 = 512 different models that can be constructed by elimination of combinations of the above nine reactions, the results of topological filtering indicated that any combination of reactions r6–9 can be removed, which coincides with all four models that were previously manually constructed by Xu et al. (6) and that corresponded to elimination of r6, r7–9, r6–9, or none of these reactions. Our systematic analysis extends these previous results in two aspects. First, although not included in the models of Xu et al., our analysis indicated that models in which only some of the reactions, r7, r8, and r9, are eliminated may be interesting candidates for experimental analysis. Although none of these reactions was apparently crucial to explain the ERK data, individual reactions with small impact on the dynamics may still be present. However, the response for some of the resulting models is indistinguishable from the model with r7–9 eliminated; for a more detailed analysis, one would have to measure additional variables or include other (chemical) perturbations of the network. Second, topological filtering showed that all the reactions, r1–5, involved in feedback regulation of ERK signaling were essential to explain the data, indicating the importance of these control loops. Overall, even for this relatively large ERK signaling network, our automatic method reduced the space of biological hypotheses to a set of 24 = 16 viable models, about 3% of the 512 possible models.

Inference of mechanisms for the dynamics of Msn2

For a novel application of the topological filtering method, we focused on the short-term dynamics of the yeast transcription factor Msn2 and its control by the PKA pathway (Fig. 3A). Glucose addition to starved cells leads to the transient accumulation of cAMP, which peaks after about 1 min with an increase by one order of magnitude from the basal cAMP concentration under nutrient stress (23) (Fig. 3B). Increased PKA activity then inhibits Msn2 activity by triggering Msn2 accumulation in the cytosol. Msn2 has at least four phosphorylation sites in the NLS region (Ser582, Ser620, Ser625, and Ser633) and one in the NES region (Ser288) that are known targets of PKA in this process (16). The stress-generated signals affect both nuclear export and nuclear import, leading to localization changes of Msn2 that are completed within minutes after the external glucose concentration has been altered (14, 24) (Fig. 3C).

Fig. 3 Msn2 model, experimental data, and reduced models that are compatible with the Msn2 translocation data.

(A) Sketch of the model for the nuclear translocation dynamics of Msn2 with reactions (arrows) and components (rectangles, proteins; ellipses, small molecules) and their localization in the cytoplasm and nucleus (inside and outside the shell, respectively). The model is used to predict and identify mechanisms for the Msn2 dynamics upon stress release after addition of glucose to cells that have been starved for 20 min. (B and C) Experimental data for cAMP (B) are the input to the model, and fractional nuclear localization of Msn2 is the readout (C). Different time courses pertain to the wild type (stars), gpa2Δ (squares), and vma5Δ (diamonds) as originally published in (26). (D) Relationship between viable models that are automatically generated from the original Msn2 model (top) based on the removed parameters (arrows with annotation of the eliminated reaction). Models are named, with the number of the model subscripted and the removed parameters superscripted. (E) Model probabilities given the Msn2 translocation data (blue), after incorporation of the phosphorylation data (yellow), and with the msn5Δ data (black). OM, original model.

Although the general topology of the signaling systems is well known, this does not apply to the detailed control mechanisms and their relative importance or to most of the quantitative and dynamic aspects of the response that are critical for modulating downstream gene expression (17). The localization of Msn2 is determined by nuclear export and import and, consequently, by the net flux (the difference of the two transport rates). Because PKA simultaneously controls the nuclear import and export of Msn2, either of the two processes could be critical for establishing changes in Msn2 localization. Moreover, it is unclear whether PKA exerts its influence primarily in the nucleus or in the cytoplasm. For instance, the hypothesis that the control of nuclear export is the most important process (24) is intuitively hard to reconcile with the observation that the cAMP-inhibited regulatory subunit of PKA (Bcy1) predominantly localizes to the nucleus in cells grown in the presence of glucose, which would predict that PKA activity in the nucleus would be low under these conditions (25). The distribution of Msn2 between compartments and, correspondingly, the control of localization are graded because, even under poor growth conditions, the majority of Msn2 remains in the cytosol [we estimated a cytosolic fraction of 57% under the starvation conditions in (26)]. Dynamic shuttling of Msn2 between the cytosol and the nucleus could enable rapid adaptation to favorable growth conditions (24).

To investigate these quantitative and dynamic aspects of Msn2 control, we focused on the short-term dynamics of the PKA-induced cytosolic and nuclear phosphorylation of Msn2. To incorporate the possible mechanisms of the PKA-mediated control of Msn2 discussed above, we developed a dynamical model with four state variables representing the fractions of nonphosphorylated Msn2 (Msn2cyt and Msn2nuc) and phosphorylated Msn2 (Msn2Pcyt and Msn2Pnuc) in the cytosol and nucleus (Fig. 3A and text S3). Because PKA targets several Msn2 phosphorylation sites, we assumed that the sites in the NLS are functionally equivalent such that one can combine them into a single state variable. The kinetics of the model are similar to (27). We used the experimentally measured cAMP dynamics (Fig. 3B) as an input signal, which allowed us to consider the system in isolation from other regulators of PKA activity: cAMP constitutes the only known direct link between the Msn2 system and the rest of the PKA pathway, such that feedback loops of the PKA pathway may affect Msn2 localization only through the concentration of cAMP. The model has nine reactions; with Michaelis-Menten kinetics for reactions r5, r6, r7, and r8 that each require two parameters and the reversible reaction r9, the model has 14 parameters. The two parameters for cAMP-induced PKA activation are necessary to trigger a response, and at least one of the two reactions for Msn2 phosphorylation is also required. Excluding the four Michaelis-Menten constants and the three necessary parameters, at most seven parameters can be eliminated simultaneously, resulting in 3 × 26 = 192 possible network topologies to be discriminated. Hence, although the model is relatively small [note that “small model” is not well defined by the community; see, for example, (28)] and the general topology is well characterized, a detailed mechanistic understanding requires the analysis of a model ensemble that is too large to handle manually.

Elimination of a majority of topologies by Msn2 localization data

We focused on the short-term dynamics of Msn2 during the first 3 min and used the corresponding quantitative data from a previously published experiment in which 100 mM glucose was added to cells that had been starved for 20 min (26). Time series data were available for the cAMP concentration (Fig. 3B) and for the localization of Msn2 in the wild type and in gpa2Δ and vma5Δ mutants (Fig. 3C). These deletions compromise separate branches of the PKA pathway (26, 29). In all strains, the addition of glucose resulted in a net translocation of Msn2 from the nucleus to the cytosol (Fig. 3C).

We first applied our method to the Msn2 model using data for wild-type cells, which generated 47 models that were consistent with this data (out of 192 possible candidate models). To further restrict the number of candidate models, we then included data for the two mutants as well. This yielded 11 viable candidate models that, together with the original model, defined an ensemble of 12 candidate models, none of which could be reduced to a model that was not already in the ensemble (Fig. 3D). To obtain independent evidence for the set of viable candidate models, we also applied a global optimization algorithm to each of the possible 192 models (see text S4 for details). The global optimization found parameter points for which a model was consistent with the Msn2 localization data in most of the candidate models obtained by topological filtering, but for none of the other models (fig. S1). This indicates that topological filtering identified the correct topologies despite otherwise difficult parameter estimation (and, consequently, model selection) problems.

Comparison of the candidate topologies (Fig. 3D) shows that there is no model that can be reduced by more than two parameters. We also note that in 9 of the 11 viable models, only parameters in a single reversible reaction were eliminated. In the remaining two models [M8(48) and M10(37), where the superscript in Mi(xy) denotes that the parameters corresponding to reactions rx and ry were eliminated and the subscript represents the model number], the eliminated parameters stem from successive reactions in the same direction in the nuclear (de)phosphorylation of Msn2 and the transport of phosphorylated Msn2. This indicates that a large and nonarbitrary part of the network topology is required to capture the Msn2 localization data.

The model probabilities for the model ensemble based on the Msn2 localization data only (blue bars in Fig. 3E and fig. S2) identify model M6(78), without reactions r7 and r8, as the most plausible topology given the experimental data. However, evidence supporting M6(78) is not much stronger than that for models Mfull(), M1(1), M5(8), and M8(48). Thus, although using the localization data greatly reduced the number of biological hypotheses, it alone was not sufficient to decide in favor of one particular model and, correspondingly, in favor of one dominant mechanism controlling Msn2 localization. Therefore, we continued with iterative model discrimination based on new model predictions and on new experiments.

Further discrimination between topologies by Msn2 phosphorylation dynamics

To design informative follow-up experiments, we used the set of 12 viable models to predict features controlling Msn2 that were not included in the original data set. All models [except for M3(2), which is not the most plausible model] predicted that Msn2 phosphorylation increases 20 to 40 s after glucose addition, and that phosphorylation reaches a stable plateau after another 15 s [see Fig. 4, A and B, for models M1(1) and M4(7)]. These predicted dynamics are unexpectedly fast because we observed much slower changes in cAMP concentration and Msn2 localization (Fig. 3, B and C).

Fig. 4 Predictions and experiments of Msn2 phosphorylation.

(A and B) Predicted dynamics of absolute amount of phosphorylated Msn2 species in model 1 (A) and model 4 (B) after addition of glucose. Purple (cytosol) and orange (nucleus) regions show the 95% confidence interval of the dynamics of phosphorylated Msn2 concentrations in these two compartments for all viable parameter points, and yellow and green lines indicate the corresponding average dynamics for the viable points. (C and D) Model predictions of total Msn2 (nuclear and cytoplasmic) phosphorylation relative to time zero for model 1 (C) and model 4 (D). Average predictions are given by red lines, blue regions show the 95% confidence interval over the predictions of the viable parameter points, and black symbols denote experimental phosphorylation data (means ± SD). The experimentally measured Msn2 phosphorylation sites are Ser633 (triangles) and Ser620 (squares) in the NLS and Ser288 in the NES (circles).

When analyzed in more detail, we found that models M1(1)and M3(2) behave similarly in that they predict low Msn2 phosphorylation under starvation and intermediate phosphorylation (25 to 75%) of total Msn2, 3 min after glucose addition (Fig. 4A and fig. S3). In contrast, models M4(7) to M11(34) predict that about half of Msn2 is phosphorylated under starvation and that the vast remainder becomes phosphorylated when glucose is added (Fig. 4B and fig. S3). Despite the few experimental data used so far for model identification, the variability in the predictions for the absolute amount of the nuclear and cytosolic phosphorylated Msn2 is rather low (colored areas in Fig. 4, A and B). The predictions of Msn2 phosphorylation dynamics by model M2(1,2) are somewhere in between those of the first group [M1(1) and M3(2)] and the second group [M4(7) to M11(34)], with the nonnegligible initial phosphorylation of Msn2 slightly increasing over time [although with large uncertainty (fig. S3)].

We reasoned that this low variability might enable accurate model discrimination by comparing the models to additional experimental data. However, to design corresponding experiments, we needed to take into account that phosphoproteomics measurements typically result in relative (with respect to a reference, in this case the phosphorylation at the start of the experiment), and not absolute, quantifications. Especially models M1(1) and M3(2) gave more variable predictions for the changes in relative Msn2 phosphorylation because variability in low initial amounts translates into larger uncertainties. For model M1(1), this can be seen from the blue regions in Fig. 4C [see fig. S4 for M3(2)]. A lower variability in models M4(7) to M11(34) with lower fold changes of Msn2 phosphorylation, however, still implied distinct predictions from the first group [see Fig. 4D for the example of M4(7), and fig. S4]. Together, these results indicated a possible discrimination between the remaining models by quantitatively measuring the fold change of Msn2 phosphorylation at high time resolution.

Confirmation of the predicted, fast Msn2 phosphorylation dynamics with targeted phosphoproteomics

We conducted a targeted phosphoproteomics experiment to determine the overall phosphorylation state of Msn2 using selected reaction monitoring (SRM) assays (30) to monitor the relative abundance of Msn2 phosphorylation at two phosphorylation sites in the NLS (Ser620 and Ser633) and on a single site in the NES (Ser288). The remaining phosphorylation sites were not suitable for SRM assays because of technical limitations (either the peptides were too long or they had other potential phosphorylation sites). We quantified dynamic changes in the phosphorylation of the three sites within the first 3 min after glucose addition to starved cells. Synthetic isotopically labeled phosphorylated peptides corresponding to each individual phosphorylation site were used as internal reference peptides for accurate relative quantification (see Materials and Methods and table S1). We also measured the phosphorylation states of a number of sites on other proteins as negative controls to validate the relevance of Msn2 phosphorylation (fig. S5). The phosphoproteomics data qualitatively confirmed the predictions of all models that the amount of phosphorylated Msn2 increases between 20 and 40 s, and then reaches a plateau (Fig. 4, C and D, fig. S4, and table S2). Moreover, the phosphorylation state of all the measured phosphorylation sites followed a similar pattern (Fig. 4, C and D).

To discriminate between the models with these new data, we computed 95% confidence intervals for the experimental phosphorylation data, which indicated that a fraction of Msn2 phosphorylation larger than 26% before glucose addition is inconsistent with the data (see Materials and Methods for details). On the basis of this information, we evaluated all 12 candidate models and excluded those parameter points for each model that predicted a higher initial phosphorylation. The updated model probabilities are shown in Fig. 3E as yellow bars (see fig. S2 for values). As a result of this analysis, M1(1) becomes the most plausible model. However, M1(1) is only slightly more than twice as probable as the original model. Hence, we needed to iterate our procedure further to fully discriminate between the candidate topologies.

A single topology after incorporating translocation data for the msn5Δ mutant

To achieve this discrimination by further model predictions and experimental data, we focused on the role of Msn5 in the control of Msn2 translocation. Briefly, PKA deactivates Msn2 by phosphorylation, and PKA-mediated Msn2 phosphorylation in the NES induces nuclear export by the exportin Msn5. Evidence for this role of Msn5 stems from the observation that Msn2 remains in the nucleus under high-glucose conditions in msn5Δ mutants (14, 31). If we assume that, without Msn5, phosphorylated Msn2 leaves the nucleus at the same rate as unphosphorylated Msn2, we can simulate the Msn2 localization dynamics in the wild type and in the Msn5 deletion mutant (see Materials and Methods for details) and compare the model predictions with experimental data on the nuclear export of Msn2.

Experimentally, we observed that nuclear export of Msn2 upon glucose addition to starved cells is largely delayed in msn5Δ cells, and that a substantial fraction of Msn2 remains in the nucleus even in the absence of stress [Fig. 5A; see also (14, 32)]. Predictions of the relative change in nuclear localization of Msn2 in the wild type and in msn5Δ suggested that only M1(1) is consistent with the data, as can be seen by comparing the overlay of model predictions (lines and shaded areas) and experimental data (symbols) in Fig. 5B for M1(1) and in Fig. 5, C and D, for M2(12) and M4(7), respectively (see fig. S6 for the other models). To quantify this consistency, we computed posterior model probabilities (black bars in Fig. 3E, and fig. S2) for all models after incorporation of the msn5Δ data. The probability that M1(1) most accurately reflects the biology underlying regulation of Msn2 localization in response to release from glucose starvation, given all the available data, is 99%. We therefore conclude that among all models that our method generated, model M1(1) best describes and predicts the available data.

Fig. 5 Localization of Msn2 in wild type and msn5Δ.

(A) Experimental data for the relative change in the localization of Msn2 for the wild type (blue) and for msn5Δ (red). The black line indicates the concentration of glucose in the system (100 mM glucose is added at time zero, denoted by the dotted line). (B to D) Predictions (red, mean; blue, 95% confidence interval) and experimental data (green circles, mean; bars, SE) for the difference between the relative (to time zero) change in localization of Msn2 in wild type and msn5Δ for model 1 (B), model 2 (C), and model 4 (D).

Dominant control of the switch in the Msn2 phosphorylation state by nuclear processes

We used the models to quantitatively understand the contributions of phosphorylation and dephosphorylation in the nucleus and cytoplasm to the changes in Msn2 localization. The most probable model, M1(1), predicted that the increased phosphorylation of Msn2 after glucose addition primarily affected nuclear Msn2 (red and gray curves) (Fig. 6A). In this model, turnover of Msn2 phosphorylation is predominantly nuclear (Fig. 6B). In most of the other models [M4(7) to M11(34)], the cytoplasmic Msn2 species respond most prominently to the cAMP peak (fig. S7). Because model 1 fits the experimental data best, we concluded that the change of Msn2 localization is primarily due to nuclear phosphorylation of Msn2 mediated by an increase in the activity of the kinases Tpk1, Tpk2, and Tpk3.

Fig. 6 Model classes predict different mechanisms controlling Msn2 translocation.

(A) Dynamic predictions of Msn2 concentrations for model 1. Unphosphorylated Msn2 in the cytosol (yellow), phosphorylated Msn2 in the cytosol (blue), unphosphorylated Msn2 in the nucleus (red), and phosphorylated Msn2 in the nucleus (gray). (B) Absolute rates of phosphorylation in the cytosol (red), dephosphorylation in the cytosol (yellow), phosphorylation in the nucleus (blue), and dephosphorylation in the nucleus (black). (C) Net phosphorylation rates (the result of the difference in the phosphorylation and dephosphorylation rates) in the cytosol (blue) and in the nucleus (gray). In all cases, shaded areas correspond to a 95% confidence interval. (D) cAMP (dashed line; WW, wet weight) induces a fast net phosphorylation in the nucleus (R8) for the first minute upon glucose addition, and subsequent nuclear export (R4) and dephosphorylation (R6) of Msn2 between 1 and 3 min.

Notably, all models predict that phosphorylation and dephosphorylation rates in a given compartment are nearly identical at all times [see Fig. 6B for M1(1) and fig. S7 for the original model and M4(7)]. Dynamic changes in the Msn2 phosphorylation states result from very small differences between these rates after glucose addition, which can be observed from the net phosphorylation rates (Fig. 6C and fig. S7). Hence, regulation of Msn2 translocation appears to be driven by subtle changes in the ratio of phosphorylation and dephosphorylation rates.

The model-based identification suggests a sequence of events controlling Msn2 localization (Fig. 6D). With highly dynamic, constitutive cycling of Msn2 between its phosphorylated and unphosphorylated forms in the cytosol and nucleus, increased cAMP abundance leads to a fast net phosphorylation of Msn2 in the nucleus. This allows for a fast initial response to the cAMP stimulus. If Msn5 is available, the transport of phosphorylated Msn2 out of the nucleus then proceeds at a much faster rate than through diffusion only. Furthermore, the importance of phosphorylation-driven nuclear export suggests that mainly the NES, and not the NLS, of Msn2 controls the short-term dynamics of the transcription factor after release from glucose deprivation stress. Note that the transport of unphosphorylated Msn2 between the nucleus and cytosol is slow or absent in model 1. These transport reactions may not be negligible for the changes in Msn2 localization that occur on longer time scales after the reintroduction of glucose. However, model M1(1) remains the best model when compared with experimental data for the nuclear abundance of Msn2 for up to 6 min after glucose addition, which is longer than the 3 min used in the training data (fig. S8).

A sensor for the rate of change of cAMP

To mechanistically and quantitatively understand how the cell interprets changes in cAMP (a proxy for glucose), we used the best model to simulate Msn2 responses to cAMP inputs that were either sharper or more extended than the measured cAMP dynamics (Fig. 7, A to C). This means that the change of cAMP concentration per time unit is simulated to be faster or slower than measured in the real system. The responses in terms of the Msn2 total nuclear phosphorylation and localization to the cytoplasm are unexpected, because it is not the integral of the cAMP signal, which increases from Fig. 7, A to C, but the rate of change of cAMP concentration, which decreases from Fig. 7, A to C, that determines the quantitative response. Specifically, the rate of change is decoded by the peak of the net nuclear phosphorylation of Msn2: The faster the rate of the change in cAMP concentration, the faster the net dephosphorylation rate (Fig. 7, D to F). This behavior is very robust to changes in model parameters as indicated by the high-confidence predictions despite values of viable parameters that spread over several orders of magnitude (fig. S9).

Fig. 7 Mechanisms controlling Msn2 localization.

(A to C) Prediction of the nuclear localization of Msn2 (yellow, mean; green, 95% confidence interval), and fraction of phosphorylated Msn2 (red, mean; blue, 95% confidence interval) for three different cAMP inputs (dotted line) using model 1. cAMP has been scaled to 1/10 of the cAMP function used in predictions for convenient visualization. (D to F) Prediction of the net rates of phosphorylation of Msn2 in the nucleus (black, mean; gray, 95% confidence interval) and in the cytosol (red, mean; blue, 95% confidence interval) for three different cAMP inputs using model 1. (G) Steady-state predictions of model 1 for the nuclear localization of Msn2 (black) and Msn2 phosphorylation (yellow) as a function of the input (cAMP). (H) Protein cycles involved in the control of Msn2 translocation that operate predominantly through export of the unphosphorylated (blue) or phosphorylated (red) forms of Msn2 and reactions with low flux in either cycle (black). (I) Steady-state predictions of the simplified model for the nuclear localization of Msn2 (black) and relative phosphorylation of total Msn2 (yellow) caused by different relations between the input and relative phosphorylation of Msn2 in the nucleus (red) and the cytoplasm (blue). Relations between input and nuclear localization of Msn2 for regimens with 0% (dashed blue) and 100% (dashed red) phosphorylation in the nucleus are shown.

Mechanistically, a classical theoretical proposal for how a rate-of-change sensor could be established in biology bore resemblance to the Msn2 system in that it postulated a diversity of thresholds for the release of “packets” of signaling molecules from different intracellular vesicles (33). Here, these signaling molecules could correspond to nuclear Msn2. To investigate the hypothesis that multiple thresholds are involved in controlling Msn2 localization, we used model M1(1) to simulate the steady-state response of Msn2 to constant amounts of cAMP, which predicted a nonmonotonic response of Msn2 localization to cAMP, with a threshold behavior at intermediary concentrations (Fig. 7G, black curve). In contrast, total Msn2 phosphorylation was predicted to increase monotonically (Fig. 7G, yellow curve). These different responses implied that Msn2 phosphorylation may either increase or decrease the amount of nuclear Msn2 and that steady-state measurements of either quantity alone do not necessarily characterize the system’s state.

To enable a more intuitive understanding of these phenomena, we developed a simplified mathematical model that captures several key features of the realistic model (text S5). The simplified model has seven parameters: Parameters k2, k3, and k4 have the same interpretation as in model M1, and the four other parameters connect cAMP to the phosphorylation of Msn2 in the nucleus and cytosol. Structurally, it allows for cycling between four protein states with constitutive nucleocytoplasmic shuttling and protein phosphorylation and dephosphorylation in each compartment (Fig. 7H). The model incorporates relations between parameters inferred from the in vivo data (fig. S9 and text S5), namely, (i) fast phosphorylation and dephosphorylation reactions compared to transport reactions; (ii) high-affinity protein modifications in the nucleus, which lead to a switch-like relation between input and relative phosphorylation state (34); (iii) a sigmoidal relation for the cytoplasmic phosphorylation state; and (iv) faster transport of phosphorylated than unphosphorylated protein species. The simplified model provides a good representation of the realistic model’s qualitative steady-state response to input (compare Fig. 7, G and I), but it cannot represent the response to very low or very high inputs or the rate-of-change sensor behavior (see text S5). Simplicity, however, allows for a mathematical analysis that reveals the following: The input-output relation is a composition of two regimens (indicated by dashed lines in Fig. 7I), each of which is sigmoidal with apparent “affinity” constants determined by transport parameters, and a switch between them. The regimens correspond to two possible cycles with flow between the protein species (Fig. 7H). One cycle predominantly exports the unphosphorylated protein form from the nucleus (blue), whereas the other cycle relies on nuclear export of the phosphorylated form (red). Their interplay generates the nonintuitive relation between input and localization of the protein.

Discussion

Here, we introduced a novel computational method, topological filtering, to discriminate between many mechanistic models of a biochemical system and to automatically identify those parts of a model that contribute to the observed dynamics. With the topological filtering approach, only models that pass the filter have to be further analyzed, thus drastically reducing the computational effort needed to infer and predict previously unknown biochemical mechanisms, and enabling systematic hypothesis testing when the sheer number of possible models is prohibitive for manual construction. Indeed, when applied to growth factor–stimulated ERK signaling and yeast glucose signaling, we found that only small subsets of hypotheses were consistent with the data, even if the mechanisms of interest were embedded in large networks as for ERK signaling.

To infer dynamic models in systems biology, typically only sparse sets of experimental data are available, which may result in large uncertainties for model parameters and for model predictions. Our method tackles this problem by systematically studying sets of parameter points that are compatible with the data. This allows us to assess the relative plausibility of a set of models by a Bayesian approach, which can be computationally efficient if the analysis of parameter spaces is based on methods like the one in (20). For general optimization, bounds on parameter spaces in “sloppy” models can cause inference methods to fail (35), but this caveat does not apply here because of our gradient-free Monte-Carlo exploration of the parameter space. Topological filtering explores the viable parameter space to generate model representations through elimination of parameters (hypotheses). Thus, our method simultaneously characterizes parameter spaces and reduces mathematical models. Local methods have been proposed to use the lack of parameter identifiability as a criterion for model reduction (36, 37). One could use topological filtering as a basis for a general model reduction method that could incorporate established reduction approaches (38), such as sensitivity analysis (39), time scale analysis (40), and lumping (41).

Note that topological filtering depends on the initial set of hypotheses incorporated by the modeler; it cannot add new hypotheses. However, established methods for the inference of coarse-grained models from data (42) could be used to establish an unbiased initial set of hypotheses. Alternatively, the identification of stochastic differential equations provides opportunities for elucidating parts of the model that need to be augmented (43). As another potential caveat, we note that our approach to model reduction is “greedy.” Although based on nested models and systematic evaluation of stepwise reductions for each candidate model, finding all viable candidates cannot be guaranteed. Also, large parameter spaces affect both parameter exploration and model reduction, such that the scalability of the method beyond the model sizes considered here needs to be further investigated. Future developments could therefore focus on computational efficiency and scalability, for instance, by exploiting advanced numerical methods (44).

Topological filtering, coupled with rationally designed experiments, enabled novel insights into the control of Msn2, including the unexpected discovery of a rapid switch in Msn2 phosphorylation after glucose addition that occurs well before the end of the cAMP peak. We identified a single model topology that was quantitatively consistent with the data. In contrast, analysis of the volume of the viable space of different models as a proxy for their robustness to parameter perturbations did not yield such a unique identification (fig. S10 and text S6). The model identified by topological filtering enabled in-depth analysis of the mechanisms controlling Msn2 localization and predicted that nuclear phosphorylation and subsequent nuclear export of Msn2 are the driving forces of the translocation, and that the control of Msn2 localization relies on subtle changes in phosphorylation rates against a background of high constitutive nucleocytoplasmic cycling of the transcription factor.

In systems terminology, the rate-of-change sensing function of Msn2 corresponds to a “differentiator,” which is a fundamental class of dynamic systems, for which no real biochemical network example is known so far (45). In the physiological context of Msn2, two properties of a differentiator are important. Signal differentiation enables fast reactions to fast changes in the cell state, which is required for stress responses to fluctuating environmental conditions. However, as a trade-off, differentiation will amplify stochastic fluctuations of the signal. Another study (46) identified genes controlled by Msn2 or its homolog Msn4 in yeast as a major “noise regulon”—a set of genes that show correlated activity fluctuations—and showed that the associated dynamic responsiveness provides survival advantages under stress conditions. In addition, on time scales of tens of minutes or longer, the same mechanisms controlling Msn2 translocation that we investigated can establish a versatile signal processing module (47). These experimental findings align well with the predicted system-level differentiator function of Msn2.

Relatively few characteristics of the mechanisms controlling localization can lead to nontrivial manifestations of dynamic phenomena at the systems level. Constitutive nucleocytoplasmic cycling, control of localization by (fast) phosphorylation, and the response to highly dynamic signals are hallmarks of other transcription factors, such as Crz1 in yeast (48) or forkhead box protein O (FoxO) (49), nuclear factor κB (NF-κB) (50, 51), and nuclear factor of activated T cells (NFAT) (52), which are transcription factors in higher organisms. In most of these cases, both the quantitative mechanisms controlling their activity and the more abstract signal processing functions are poorly characterized. Our model-based analysis suggests general design principles that could establish a temporal signal differentiator and thereby extend the repertoire of known “motifs” for decoding temporal cellular signals (53). The combination of mechanisms controlling localization is also versatile. For instance, with different settings for the reactions, it can operate as a static bandpass filter (fig. S11, text S5, and table S3).

Of course, as for Msn2, the control processes involved are more complicated than our abstracted models; they involve multiple phosphorylation events, the formation of various complexes, and regulated degradation. More detailed mechanistic insights could be obtained by including them in the model-based analysis, but direct experimental evidence for many of these details of intracellular signal processing systems is not available or technically obtainable. This is one of the reasons why our model-based method for the inference of the mechanistic underpinnings of dynamic biological systems appears promising for application to other examples where mechanisms and quantitative features of a system are substantially uncertain.

Materials and Methods

Targeted phosphoproteomics of isotopically labeled synthetic peptide standards

Crude isotopically labeled ([13C,15N]lysine-labeled or [13C,15N]arginine-labeled) synthetic phosphorylated peptides were purchased from JPT Peptide Technologies. These peptides represent unpurified products of a high-throughput Spot synthesis and lack a precise peptide concentration determination. Hence, they can be used for determination of relative phosphorylation site abundance changes but not for absolute quantification. For the target protein Msn2, the following three phosphorylated peptides were purchased: RPS[phospho]YR (NLS site, Ser620), SS[phospho]VVIESTK (NLS site, Ser633), and RFS[phospho]DVITNQFPSMTNSR (NES site, Ser288). As negative controls, synthetic, isotopically labeled phosphorylated peptides from Pfk1 (DAFLEATS[phospho]EDEIISR), Pfk2 (NAVSTKPTPPPAPEASAES[phospho]GLSSK), and Hxk1 or Hxk2 (KGS[phospho]MADVPK) were used. Optimized SRM assay parameters for each phosphorylated peptide were obtained by initial investigative SRM measurements of the synthetic peptides and include (i) selection of the most intense and selective transitions per peptide (four to nine transitions; see table S1), (ii) relative transition intensities, and (iii) retention time information. Subsequently, all synthetic peptides were spiked into the phosphorylation-enriched samples of interest in constant amounts, roughly adjusted to the endogenous peptide abundance, and served as internal references.

Cell culture and sample preparation

Cells were grown in synthetic media as described (26). Yeast cells expressing Msn2-GFP (green fluorescent protein) were grown to mid-exponential phase (optical density, 0.7 to 1), harvested, washed twice, and resuspended in synthetic complete (SC) medium without glucose. Twenty minutes after the first wash with SC medium, glucose was added to a final concentration of 2%, and samples were withdrawn at the indicated time points. Cells were harvested after quenching with trichloroacetic acid (6.25%, final) and washing with ice-cold acetone. Yeast cells were lysed in lysis buffer [8 M urea, 100 mM NH4HCO3, 5 mM EDTA, 1 mM tris(2-carboxyethyl)phosphine (TCEP), pH 8.0] with glass beads. Cell debris was removed by centrifugation, and protein content was determined with a bicinchoninic acid protein assay (Pierce). For each sample, 2.0 mg of total protein was reduced (5 mM TCEP), alkylated (70 mM iodoacetamide), digested with trypsin (Promega), and prepared for a phosphorylated peptide enrichment procedure as described previously (54). Briefly, phosphorylated peptides were enriched with titanium dioxide (GL Science), eluted with 0.3 M NH4OH (pH 10.5), and purified using C18 cartridges (C18 Micro Spin columns, The Nest Group Inc.). Finally, the phosphorylated peptide mixtures were dried, resolubilized in 0.1% formic acid, and immediately analyzed. All samples were processed in parallel.

Targeted mass spectrometry

All SRM measurements were performed on a TSQ Vantage QQQ mass spectrometer (Thermo Fischer Scientific) equipped with a nano-electrospray ion source. Typically, 1 μg of phosphorylation-enriched peptides was loaded onto a 75 μm × 10.5 cm fused silica microcapillary reversed-phase column packed with Magic C18 AQ material (200-Å pore, 5-μm diameter, Michrom Bioresources). A linear 40-min gradient from 2 to 46% solvent B (solvent A: 98% water, 2% acetonitrile, 0.1% formic acid; solvent B: 98% acetonitrile, 2% water, 0.1% formic acid) at a flow rate of 300 nl/min was applied for phospho-peptide separation. The mass spectrometer was operated in the positive ion mode using electrospray ionization with a capillary temperature of 280°C, a spray voltage of +1350 V, and a collision gas pressure of 1.5 mtorr. SRM transitions were monitored with a mass window of 0.7 half-maximum peak width (unit resolution) in Q1 and Q3. All SRM measurements were performed in scheduled mode with a retention time window of 3.5 min and a cycle time of 1.5 s. The collision energy for each transition was calculated using the formula CE = 0.034 × m/z − 0.848 for doubly charged precursor ions and CE = 0.022 × m/z + 5.953 for triply charged precursor ions (in-house optimized formula, m/z is mass-to-charge ratio of the precursor ion).

Targeted mass spectrometric data analysis

All obtained SRM traces were analyzed with the software Skyline (55). Interfered or noisy transitions were removed manually. For quantification, the ratio between a given endogenous (light) and its isotopically labeled reference peptide (heavy) was calculated from the sum of all light and heavy transition peak areas. Because the reference phosphorylated peptide amount was kept constant through all samples, endogenous abundance changes between samples could be determined. All data were normalized relative to the mean of the starved condition from three biological replicates (table S2).

The SRM data set associated with this manuscript has been deposited to the PeptideAtlas SRM Experiment Library [PASSEL; http://www.peptideatlas.org/passel/ (56)] and is accessible from ProteomeXchange with identifier PXD000236 (http://proteomecentral.proteomexchange.org/dataset/PXD000236).

Imaging

Analysis of Msn2-GFP localization was performed as described (26). In brief, cells expressing Msn2-GFP were grown to mid-log phase and loaded into commercially available microfluidic chips. Cells were subjected to glucose starvation for 20 min before readdition of glucose (2%, final), and Msn2-GFP localization was followed by time-lapse imaging. Cell segmentation and analysis of Msn2-GFP localization was performed through in-house MATLAB (MathWorks)–based software (57). Cells were segmented based on the bright-field image, and Msn2-GFP localization was quantified by calculating the coefficient of variation (CV) of pixel intensities of the 500 brightest pixels per cell in the GFP channel, which allows for a relative quantification of the changes in Msn2-GFP localization over time in different genotypes but effectively excludes signals from the vacuole. Strong nuclear localization causes high CV because of bright pixels in the nucleus and reduced intensity in the cytoplasm; uniform distribution of Msn2 between the nucleus and cytoplasm leads to a reduced CV. To calibrate the model, we estimated the relative concentrations of Msn2 in both compartments by comparing the total intensity of the nuclear pixels to the total intensity of all pixels of the cell (with the background intensity subtracted) using the CellX image analysis software (58) and a nuclear volume of 7% of the cell volume (59). For each strain, the cells from three independent experiments were averaged and plotted as means ± SEM.

Model definition

We consider ODE models for systems with n state variables (for example, component concentrations for biochemical networks), ζ reactions, and d model parameters of the formM(θ)={dxdt(t)=Nv[x(t),θ,u(t)],θ=[θ1,θ2,,θd]Tyk=h[x(tk)]+e(tk),e(tk)N(0,Sk),k=1,2,,K(1)with the state variables x(t) ∈ ℝn, the stoichiometric matrix N ∈ ℤn×ζ , the reaction rates or fluxes v(·) ∈ ℝζ, the potentially time-varying inputs u(t), and the vector of model parameters θ ∈ ℝd. We assume that the system output yk ∈ ℝl at time point tk is generated by a nonlinear function h(·) of the system state variables x(t) and an additive contribution of measurement noise e(tk) ∈ ℝl. Furthermore, we assume that this noise is normally distributed N(0,Sk) with covariance matrix Sk ∈ ℝl×l.

Exploration of parameter space

Given a set of time-discrete observational data Y=[y^1,y^2,,y^K]T, we assume that the residuals εk=(yky^k)l. That is, the deviations between model and data at any time point are Gaussian distributed. Then, the likelihood of observing these data, given the model M and the parameter point θ, is given byp(Y|θ,M)=e12εTS1ε(2π)l|S|(2)where ε=[ε1T,,εKT]T=[(y1y^1)T,,(yKy^K)T]TKl and the measurement covariance matrix S ∈ ℝ(Kl)×(Kl).

For evaluating the compatibility between model and data, we use the negative log-likelihoodfunctionE(θ,Y,Mi)=ln[p(Y|θ,M)](3)as a cost function. We derived a formal criterion to classify parameter points as viable (see text S2 for details).

Generation of candidate models

For every parameter point θ=[θ1,θ2,,θm1,θm]TVexpli, the set of all viable parameter points found in the exploration, the algorithm constructs a new set of m parameter points, which is defined by projection of θ onto each Cartesian axis. That is, the ith parameter point of this set has the ith component equal to zero and the other components are the same as in θ: θ={θ1,θ2,,θm1,θm}Vexpli{{0,θ2,,θm1,θm}{θ1,0,,θm1,θm}{θ1,θ2,,0,θm}{θ1,θ2,,θm1,0}(4)

The next step of the method is to check the viability of each projected parameter point. If the ith projected parameter point is viable, it implies that the kinetic parameter θi is not essential to explain the data. Subsequently, the algorithm checks the viability of a parameter point in which all previously identified nonessential kinetic parameters are set to zero simultaneously. If this parameter point is viable, the corresponding model is included in the set of reduced models. The complete set of reduced models generated from the original model is obtained by repeating the same procedure for every parameter point in Vexpli.

Ranking of candidate models

The method uses the prior probability distributions p(Mi) and p(θ|Mi) for the probabilities of the model Mi among the set SM and for the parameter point θ given the model Mi. Both priors are assumed to be uniform (that is, an equal probability for all possible models and parameterizations is assumed before incorporating experimental data). For model ranking, we assess the ratio of plausibility between two models, Mi and Mj, with the so-called Bayes factor, Bij, by integrating overall viable parameter points. With equal priors for the parameter points and models, Bij takes the formBijp(Y|Mi)p(Y|Mj)=Θip(Y|θ,Mi)p(θ|Mi)dθΘjp(Y|θ,Mj)p(θ|Mj)dθ=VolΘjΘip(Y|θ,Mi)dθVolΘiθjp(Y|θ,Mj)dθ(5)where Θi is the parameter space of model Mi and VolΘ is the volume of this space. If we compute the Bayes factors for all models with respect to model Mi, the posterior probability of this model is given byp(Mi|Y)=(1+j=1,ijCBji)1(6)

Small example model

The model has four states and seven parameters (see text S1). Parameters were explored in the region [10−4, 104] for each parameter over the parameter space Θ7 = k1 × k2 × … × k7. The likelihood function over the 20 time points, given by Eq. 2, is formulated asp(Y|θ,M1())=e12εTS1ε(2π)2|S|(7)where the residuals ε=[ε1T,,εKT]T40 with εk=(yky^k)2, y^k=[A^(tk),C^(tk)]T, and yk(θ,M1())=[A(tk,θ,M1()),C(tk,θ,M1())]T. The covariance matrix Sk40×40 is diagonal, which allowed us to write the viability condition in the explicit form given by expression (eq. S8).

Model of the ERK pathway in PC12 cells

Our analysis started from model 4 in (6), which incorporates all the mechanisms of models 1, 2, and 3. The model was downloaded from the supplementary materials of (6): http://stke.sciencemag.org/cgi/content/full/sigtrans;3/113/ra20/DC1, in systems biology markup language (SBML) extensible markup language (XML).

The initial parameter point used in the parameter exploration with our method was taken from the best model in (6) (model 2). Because the 17 kinetic parameters of the investigated reactions have not been experimentally measured, we decided to explore a wide region ([10−5, 105]).

The initial conditions of the state variables were set as in (6) with EGF(0) = 1000. Altogether, 11 experiments were conducted in (6). Three of these are control experiments in which EGF was not added at time point zero, which did not generate a response in the phosphorylation of ERK. In the other experiments, EGF was initially added with or without a combination of the three drugs cilostamide (a cAMP phosphodiesterase inhibitor), a PKA antagonist (PKAA), and a EPAC antagonist (EPACA), which gives a total of 23 = 8 experiments (in one experiment, only EGF was added). The initial conditions of the drugs, in experiments in which they are included, were set to cilostamide(0) = 100, PKAA(0) = 100, and EPACA(0) = 100 (0 if not included). All eight experiments were simulated for each parameter point evaluated. The viability of a parameter point was determined from the negative log-likelihood of all experiments and based on our viability criterion (see text S2).

Msn2 model

For details on the construction of the model, see text S3. The experimental readout with respect to Msn2 localization is the fraction of the total amount of Msn2 that resides in the nucleus over time: y(t)=Msn2nuc(t)+Msn2nucP(t)[Msn2cyt(t)+Msn2cytP(t)+Msn2nuc(t)+Msn2nucP(t)](8)

For phosphoproteomics analysis, we designed a scheme with the six most informative time points (40, 60, 80, 120, 140, and 180 s) based on in-house software for experimental design and two time points shortly before glucose is added (−60 and −120 s). Our approach to experimental design was to select time points for which the difference in the predictions of the candidate models is the largest, because these measurements are the most informative for model inference. We only consider the magnitude of the predicted concentrations, and not the corresponding derivatives. Because the models predicted the most change in the phosphorylation state in the beginning of the experiment, we also considered the 20-s time point.

For each evaluated parameter point, we first simulated the model until a steady state was reached. We performed extensive numerical testing of the model to exclude multistationarity (that is, to confirm that the steady-state solution is independent of the initial conditions). There are 14 parameters in the model: a,r,k1,k2,k3,k4,V5,V6,V7,V8,K5,K6,K7,K8. The reduced models were automatically identified by exploration of the parameter space of the original model. For each parameter, the region [10−4,104] was explored.

Incorporation of additional data: Phosphorylation and msn5Δ

To incorporate Msn2 phosphorylation data, we computed a 95% confidence interval for the increase in phosphorylation φ. This was computed as the inverse of the normal cumulative distribution function, so that the increase in phosphorylation (φ) is larger than ul with 95% probability:Pr(ulφ)=0.95(9)for each of the three measured phosphorylation sites (Ser633, Ser620, and Ser288) at time point 40 s when phosphorylation reached a new plateau. The values of ul at 40 s for the sites are Ser633: ul = 4.16, Ser620: ul = 8.50, and Ser288: ul = 3.81. All parameter points predicting that the initial phosphorylation of Msn2 is larger than 13.810.26 were, as supported by all measured phosphorylation sites, considered to be inviable for computing the Bayes factors and the corresponding model probabilities.

To predict the translocation dynamics of Msn2 in msn5Δ, we assumed that a phosphorylated and an unphosphorylated Msn2 molecule are equally likely to diffuse out of the nucleus during a fixed time interval. Correspondingly, we set the parameter in reaction r4 equal to the parameter in reaction r2 for each viable parameter point. The difference between the model predictions of the relative change in the msn5Δ and wild-type (WT) dynamics,yδ(t)=ymsn5Δ(t)ymsn5Δ(t0)yWT(t)yWT(t0)(10)where glucose is added at t0 and y is defined in Eq. 8, was then compared to the corresponding difference in the experimental data at time points 0, 1, 2, and 3 min.

Supplementary Materials

www.sciencesignaling.org/cgi/content/full/6/277/ra41/DC1

Text S1. Definition of the small example model.

Text S2. Derivation of the viability criterion.

Text S3. Definition of the Msn2 model.

Text S4. Global optimization of Msn2 models.

Text S5. Simplified Msn2 model.

Text S6. Viable volume of the Msn2 models.

Fig. S1. Distribution of objective function values for the global optimization of Msn2 models.

Fig. S2. Probabilities for the Msn2 models.

Fig. S3. Predictions of Msn2 phosphorylation.

Fig. S4. Model predictions of the relative phosphorylation of Msn2 for all 12 viable models.

Fig. S5. Phosphorylation dynamics of negative control phospho-sites.

Fig. S6. Model predictions of the difference in Msn2 localization between the wild-type and the msn5Δ deletion mutants.

Fig. S7. Detailed model predictions of phosphorylation dynamics.

Fig. S8. Longer-term predictions of Msn2 dynamics.

Fig. S9. Viable parameter points.

Fig. S10. Per-parameter fractional viable volume.

Fig. S11. Bandpass filter behavior of the simplified Msn2 model.

Table S1. Complete SRM data set comprising peak area values for all monitored transitions (Excel file).

Table S2. Calculated relative phospho-site abundance changes (Excel file).

Table S3. Parameter values of the simplified Msn2 model.

References and Notes

Acknowledgments: We thank U. Sauer, J. Buhmann, R. Aebersold, F. Rudolf, S. Kiyomi Aoki, and M. Peter for comments and discussions; S. Dimopoulos for help with image analysis; M. Matondo and N. Selevsek for TSQ instrument maintenance; and W. Kolch for sharing the experimental data for the ERK pathway. Funding: We acknowledge funding from the Swiss Initiative for Systems Biology SystemsX.ch (project YeastX) evaluated by the Swiss National Science Foundation. Author contributions: M.S., E.Z.-S., A.W., and J.S. conceived and designed the research. M.S., E.Z.-S., A.W., and J.S. wrote the manuscript. M.S., E.Z.-S., and A.G.B. developed the computational method and performed the computational analysis. C.L. designed and performed the targeted phosphoproteomics experiments, processed samples, and performed SRM data analysis. R.D. performed and analyzed experiments on Msn2 localization and prepared samples for SRM analysis. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The SRM data set has been deposited to the PeptideAtlas SRM Experiment Library (PASSEL, http://www.peptideatlas.org/passel/) and can be accessed from ProteomeXchange with identifier PXD000236 (http://proteomecentral.proteomexchange.org/dataset/PXD000236).
View Abstract

Navigate This Article