Modelling and mathematical analysis of the M 2 receptor-dependent joint signalling and secondary messenger network in CHO cells

The muscarinic M2 receptor is a prominent member of the GPCR family and strongly involved in heart diseases. Recently published experimental work explored the cellular response to iperoxo-induced M2 receptor stimulation in Chinese hamster ovary (CHO) cells. To better understand these responses, we modelled and analysed the muscarinic M2 receptor-dependent signalling pathway combined with relevant secondary messenger molecules using mass action. In our literature-based joint signalling and secondary messenger model, all binding and phosphorylation events are explicitly taken into account in order to enable subsequent stoichiometric matrix analysis. We propose constraint flux sampling (CFS) as a method to characterize the expected shift of the steady state reaction flux distribution due to the known amount of cAMP production and PDE4 activation. CFS correctly predicts an experimentally observable


Introduction
The muscarinic acetylcholine receptor (M 2 receptor) (encoded by the CHRM2 gene) belongs to the family of G protein-coupled receptors (GPCR) and is among other locations expressed in cardiomyocytes where it influences the heart beat rate (Brodde & Michel, 1999).It is related to negative dromotropic and negative chronotropic events.Its malfunctioning has been associated with a number of diseases, such as cardiomyopathies (Brodde & Michel, 1999).GPCRs represent one of the most important target classes of proteins for drug discovery (Zheng, 2006).Hence, development of specific agonists and antagonists for muscarinic receptors, including the M 2 receptor, is still of high interest.Iperoxo is a highly affine and efficacious muscarinic agonist (Schrage et al., 2013(Schrage et al., , 2014) ) that has recently served to elucidate the crystal structure of the active state of the M 2 receptor (Hu et al., 2010;Kruse et al., 2013).In traditional pharmacology, the ligand-binding event, second messenger concentrations, ion channel function, as well as tissue, organ or body responses are recorded.New label free, whole cell techniques nowadays are used to dissect signalling of intact cells into different components (Schröder et al., 2011).In addition, iperoxo and its derivatives turned out to be valuable tools for gaining deeper insight into structure-signal relationships (Bock et al., 2014;Antony et al., 2009).
Recent experimental work explored the cellular response to iperoxo-induced M 2 receptor stimulation in Chinese hamster ovary (CHO) cells (Kruse et al., 2013;Schrage et al., 2013Schrage et al., , 2015)).The cellular response was measured by dynamic mass redistribution (DMR), a technique to quantify the intracellular mass movement via optical density (Schröder et al., 2011).Since the DMR response can be assumed to be dependent on the M 2 receptor-dependent signalling our aim was to model and study the corresponding reaction system.The pathway consists of proteins as well as the secondary messenger cyclic adenosine monophosphate (cAMP).The respective biochemical reactions are principally well known (Pierce et al., 2002;Linderman, 2009;Sunahara & Taussig, 2002;Taylor et al., 2012), but to the best of our knowledge no effort has been taken so far to derive a mathematical model, especially for CHO cells, which are very important in pharmaceutical research and for the industrial production of recombinant protein therapeutics (De Jesus & Wurm, 2011;Walsh, 2010).
In this work, we developed a mass action based mathematical description of the M 2 receptor-dependent signalling network.Our developed model consists of 79 reactions, altogether involving 64 relevant proteins and secondary messenger molecules described in literature.In our joint signalling and secondary messenger model, all binding and (de-)phosphorylation events are explicitly taken into account in order to enable subsequent stoichiometric matrix and flux distribution analysis (Wiback et al., 2004).Although this kind of analysis is usually only employed for metabolic networks, our explicit modelling of binding and phosphorylation events enables the adaption of these techniques to a mixed signalling and secondary messenger system.The usefulness of applying stoichiometric matrix analysis techniques to signalling pathways has e.g.been demonstrated by Behre & Schuster (2009), who adapted elementary flux mode (EFM) analysis to this situation.We here show, how the known flux sampling technique (Smith, 1996) can be extended to incorporate partially available experimental information (here: cAMP production, phosphodiesterase 4 (PDE4) activation).We tested our combined modelling and data driven sampling method by predicting key signalling mechanisms known from literature, but not explicitly encoded into the model.Our proposed constraint flux sampling (CFS) technique allows for qualitative predictions of downstream stimulation effects on actin and tubulin levels, which here serve as markers for the mass redistribution effect.These qualitative predictions are in agreement with the experimental observations, which suggests CFS as a technique for model checking.This is further underlined by the possibility to combine CFS and EFM analysis yielding a statistical ranking of EFMs according to their expected biological relevance.

Biological Background and Network Reconstruction
GPCR-induced signalling is well-known in common literature (Pierce et al., 2002;Linderman, 2009;Taylor et al., 2012;Sunahara & Taussig, 2002).Specifically the link to the cyclic AMP (cAMP -a secondary messenger molecule)-induced signalling is in the focus of current pharmaceutical research (Milligan & Kostenis, 2006;Hu et al., 2010).Figure 1 depicts a schematic representation of the whole set of relevant molecules and their interplay, which are considered in our model.In particular, the process of receptor-induced G protein (GP) activation is well studied, where the ligand-bound receptor changes its physical structure and the inactive associated GP interacts with the receptor and dissociates into its subunits (Pierce et al., 2002).Thereby the alpha-i/alpha-s and beta/gamma subunits are activated and are able to interact independently with other proteins like adenylyl cyclase (AC) (Sunahara & Taussig, 2002;Milligan & Kostenis, 2006).The GP subunit alpha-o has no significant influence on AC but it has an influence on the DMR (Milligan & Kostenis, 2006).AC is one of the most important proteins within the GP-mediated pathway and responsible for the secondary messenger production.The large number of AC and GP subtypes causes a highly complex sub-network with many cross-reactions (Milligan & Kostenis, 2006;Sunahara & Taussig, 2002).Also the receptor activation cycle itself is not trivial.This first step in the signalling cascade is highly interesting for pharmaceutical research and led to well-developed models for receptor activation and inhibition (Woodroffe et al., 2009;Chen, 2003;Strange, 2009;Bornheimer et al., 2004).
Besides this completely membrane bound sub-network the protein kinase A (PKA)-induced phosphorylation cascade, and the feedback loop causing cAMP degradation is well studied (Taylor et al., 2012).Cyclic AMP binds to PKA and causes its activation.But an increase of PKA activity also leads to an increase in phosphodiesterase (PDE) activity, which inactivates cAMP by degrading it to AMP (Boswell-Smith et al., 2006).Through this mechanism the cell prevents a continuous overstimulation by excessive cAMP levels.Stimulation of the receptor population via the muscarinic agonist iperoxo induced a cellular DMR response at concentrations that are far lower than the corresponding concentrationbinding relationships (Schrage et al., 2013).The same authors reported this amplification phenomenon also for other ligands, including the natural ligand ACh.The exact nature of the amplification process is not understood so far, but may at least be partially attributed to intracellular signalling events (Schrage et al., 2015).According to common literature we suppose regulators of G protein signalling (RGS) and G protein-coupled receptor kinases (GRK) to be of relevance.These proteins are closely related to the deactivation of the receptor and the GP subunits (De Vries et al., 2000;Pierce et al., 2002;Hollinger et al., 2003.In this approach, we choose RGS14, GRK6 and GRK2 as important representatives for each group. Nowadays, several readouts for the stimulation response are well established, e.g. the Ca 2+ -level or the cAMP concentration (Paredes et al., 2008;Hennen et al., 2013;Gabriel et al., 2003).In addition DMR has been introduced into the pharmaceutical field in order to specifically quantify the stimulation effect on the cytoskeleton (Fang et al., 2005;Schröder et al., 2011).The DMR is an optical biosensor based procedure and measures the shift in wavelength resulting from intracellular mass movement caused by rearrangement of cell organelles and transportation processes.The optical density is very sensitive to intracellular reorganization and morphological changes of the cell, and by comparing the optical density of unstimulated and stimulated cells one can measure the specific wavelength shift and draw conclusions about the intensity of cellular response (Schröder et al., 2011).Maximal DMR response induced by iperoxo occurs typically after approximately 10 minutes (Schrage et al., 2013).According to the timescale and common literature, we did not consider transcriptional downstream responses (Shaywitz & Greenberg, 1999;Mayr & Montminy, 2001).
In this work, we chose actin and tubulin as DMR markers.Actin and tubulin are closely related to cellular movement and we assume a strong correlation between changes in both proteins and the relative wavelength shift measured by DMR (Hammond et al., 2008;Schmidt & Hall, 1998).As described in Fang et al. (2005); Strange (2009);Schröder et al. (2011) the wavelength shift is caused by intracellular mass movement.Therefore, we took all proteins directly linked to actin and tubulin into account and assumed their activation to be correlated with the wavelength shift (Fig. 1).For further references see the supplementary material.

Mathematical modelling
All interactions shown in Fig. 1 are explicitly formalized as mass action based elementary reactions (Horn and Jackson, 1972) and all known proteins and their occurring complexes are included.Hence biological information from the available biochemical knowledge is preserved.Let x 1 , ..., x n denote concentrations of all molecules in the system.Then the concentration change of molecule i can be written as where s ij is the stoichiometric coefficient of molecule i in reaction j and v j denotes the corresponding rate of reaction j.As an example, we here show the PKA activation by cAMP (Corbin et al., 1988) Here, PKA denotes the inactive form of PKA.Let us denote the rate of both reactions by v 1 and v 2 , respectively.The stoichiometric coefficients are s 11 = s 22 = −1, s 12 = s 21 = 1 and s 32 = −2.We obtain (2.6) Altogether our modelled system contains 79 elementary reactions, which can be found in the supplementary material.The full reaction system can be represented via a stoichiometric matrix S ∈ R n×m .In this matrix, every molecule is represented by one row and every reaction is represented by one column.
There are some details that need to be mentioned: Our modelled system consists of several biochemical reaction types, namely binding, stimulation and inhibition.These biochemical events need to be represented appropriately in the reaction system and the stoichiometric matrix, respectively.This was done as follows: Protein activation via phosphorylation was modelled with help of an intermediate molecule which represents the complex of the substrate and the related kinase.The kinase binds reversibly to the substrate and forms an intermediate complex which then dissociates irreversibly into the kinase and the modified substrate.For instance, GRK2 is phosphorylated by PKA (Cong et al., 2001).For this purpose, we introduce the intermediate complex PKA : GRK2 and write the reaction system as For every phosphorylation step, we assumed a backward reaction P * → P, which dephosphorylates the phospho-protein P * with the help of an unknown phosphatase.In our example, GRK2 is dephosphorylated into GRK2 GRK2 → GRK2. (2.9) Protein inhibition by kinases is modelled in a similar manner.A kinase binds reversibly to the target protein and forms an intermediate complex which then dissociates irreversibly into the kinase and the inactive protein.The inactive protein is now able to be activated again by another kinase.We illustrate this process using the inhibition of GEF by GRK2 (Eijkelkamp et al., 2010) GRK2 + GEF GRK2 : GEF (2.10) (2.11) As shown above, these reactions can be represented in a stoichiometric matrix S. The dimension of the stoichiometric matrix can be decreased by expressing the forward and backward direction of the same reversible reaction by one row where reaction rates can be both positive and negative.This is in contrast to strictly irreversible reactions where only positive reaction rates are allowed.We also used the stoichiometric model to derive a system of ordinary differential equations based on the assumption of mass action kinetics, see supplementary material.

Conservation relationships
Since signalling events are relatively fast we can assume that for each protein the overall total amount of phosphorylated, bound and unphosphorylated proteins is approximately constant, provided that the biological system is in steady state and the model was correct.Hence, checking conservation relationships is a means to check the consistency of our model.According to Palsson (2006), conservation relationships under steady state conditions are mathematically identifiable from the null space of S T .That means conservation relationships are all those vectors g for which S T g = 0. (2.12) Each entry in g corresponds to exactly one molecule.Analysis of the entries of vectors g provides thus a means to verify whether the expected constant total concentration of each protein is fulfilled in reality.Moreover, we can also obtain insights into possibly existing constant protein concentrations within whole reaction cascades.

Stimulation of the system
We are interested in qualitative changes upon receptor stimulation.For the unstimulated system, we assume a steady state characterized by constant concentrations of all molecules.Receptor stimulation causes a perturbation of this steady state resulting in dynamic changes of molecular concentrations.However, we assume that after some relaxation time the system will attain a supposingly differentsteady state, which is characterized by the molecular concentrations in the stimulated state.In reality, the stimulated state does not need to be a dynamic equilibrium in the strict sense, but we believe it to be a useful approximation for a situation of maximum response, where all concentrations are nearly constant over time.We believe that this working hypothesis is useful to analyse qualitative changes between the unstimulated and stimulated states, which is also supported by the fast-usually milliseconds-time scale of the signalling events in comparison to the observable duration of responses to receptor stimulation.
Mathematically, all stationary reaction rates v in the steady state-so called fluxes-are given as solutions of the underdetermined system of equations Sv = 0. (2.13) The unstimulated and the stimulated state correspond to different solutions of this equation.Our strategy will be to constrain the solution space of 2.13 using experimental data.We will then use Monte Carlo Sampling (see below) to compare possible fluxes in the simulated and unstimulated state.For a qualitative comparison, we suppose the DMR response to be given as the sum of all fluxes with known influence on the wavelength shift (2.14) Here, v i denotes the activating flux related to molecule i with influence on the wavelength shift.The sum runs over all k in-fluxes into tubulin and actin, which are considered as markers of the DMR response (Hammond et al., 2008;Schmidt & Hall, 1998;Schröder et al., 2011).

Sampling the flux polytope
Since we are interested in the general behaviour of the system without incorporating additional rate parameters, steady state solutions of the system can in principle be found through Markov Chain Monte Carlo Hit-and-Run sampling (Smith, 1996;Brooks, 1998;Price et al., 2004).A single move in the Hitand-Run sampling is performed by making from a given feasible solution a uniform randomly chosen B. ENGELHARDT ET AL.
move within the unit sphere.Afterwards the step size is adjusted such that the new solution is also feasible (Smith, 1996;Megchelenbrink et al., 2014).A solution v * is called feasible, if it satisfies Sv * = 0 (2.15) with bounds α i , β i .Note, that without further constraints, fluxes could take any real value, but in reality fluxes are bounded.Hence, we set for all reversible reactions α i = −1000 and β i = 1000 as loose bounds.For irreversible reactions we set α i = 0.The flux bounds can in principle be used to incorporate experimental data.We will modify the flux bounds to qualitatively incorporate fold changes between stimulated and unstimulated cells, as explained in the following section.

Constrained flux sampling
We incorporate partially available data of experimentally measured relative (steady state) molecular concentrations into the above described flux sampling scheme in order to make qualitative predictions about flux changes upon stimulation.The approach thus does not require detailed knowledge of kinetic rate constants.
Let ṽj denote the steady state flux of the j-th reaction in the case of an unstimulated receptor.According to the law of mass action (see section 2.2) with rate parameters k j , we have where {x i } is the set of molecules taking part in the particular reaction and xi their concentrations.Note that at this point we suppose involved reversible reactions to be split into two irreversible ones.Accordingly, the flux vj for the same reaction under stimulation can be defined, now with concentrations xi .Usually, in experiments relative concentration changes (fold changes) Note that the stoichiometric coefficients s ij in most cases are 1.The equation suggests a principal two-step procedure: (1) Perform conventional flux sampling for the unstimulated situation.This yields a set ṽj .
(2) Perform flux sampling for the stimulated situation by plugging observed fold changes into Eq.(2.19) in order to constrain sampled fluxes.
In reality it may be more appropriate to consider confidence intervals [f Min i , f Max i ] for f i because fold changes are subject to uncertainty.This can be addressed straightforwardly by replacing Eq. (2.19) by an inequality (2.20) The quantity ṽj in practice needs to be estimated from the empirical flux distribution under steady state conditions.A reasonable choice is to take the mean or median of the sample distribution plus/minus the standard deviation for that purpose.

Elementary flux modes
Schuster and Hilgetag introduced EFM analysis for characterizing the geometry of the solution polytope of the equation system Sv = 0 in a biologically interpretable manner (Schuster & Hilgetag, 1994).All solution vectors occur as linear combinations of EFMs.More specifically, Schuster & Hilgetag (1994) and Llaneras & Pico (2010) consider the convex flux polyhedral cone P(S) w j e j w j ≥ 0 . (2.21) EFMs are then defined as the extreme rays or edges of the flux cone P(S).A formal assumption made in this equation is that reversible reactions are split into irreversible ones.Each EFM can be characterized as the minimal set of reactions which are required for a sub-system to exist as a functional unit (Papin et al., 2004).These sub-systems either reflect fluxes through the whole reaction system or functional cycles within the system.Thus, analysis of EFMs allows for identifying biologically functional and interpretable 'building blocks' of the biological reaction system.In case of signalling, this also implies that without stimulation there exists no EFM representing the whole network and no EFM describing the signalling flow through it.
In this article, we combine EFM analysis with CFS: after having determined the flux distributions of the overall system in stimulated and unstimulated conditions we map fluxes to each of the calculated EFMs.This is possible because each flux corresponds uniquely to one reaction.We then compute the median of all fluxes related to a specific EFM.Since with our sampling procedure we generated a large (here: 100, 000) sample of flux vectors we obtained an empirical distribution of these medians for each EFM.The significance of the difference in these distributions between stimulated and unstimulated conditions can be assessed via a Wilcoxon rank test, yielding a P-value.Because we do not only compare one but several EFMs, multiple testing correction of P-values is performed via control of the false discovery rate (FDR) (Benjamini & Hochberg, 1995).Moreover, we estimated the median fold change between stimulated and unstimulated conditions.

Data: cAMP, PDE4 and DMR
Parts of the experimental data (dose-response relationships) were taken from Schrage et al. (2013): In that article, DMR was measured at 13 different concentrations of the M 2 receptor-specific activator iperoxo, giving rise to dose-response curves (Fig. 2).In addition to these data by Schrage et al., cAMP response to iperoxo treatment was measured here (Fig. 2).This was done after 30 minutes of iperoxo incubation with a concentration of 0.1μM = 10 −7 M, which corresponds to a full DMR response (Schrage et al., 2013).
The induced cAMP fold change, calculated as the ratio between the cAMP level related to the iperoxo concentration of 0.1μM and the 95% confidence interval of the basal cAMP level (see Fig. 2) is given by [2.22, 2.71].
In addition, we measured the activation of PDE4 after 30 minutes for the iperoxo concentration of 0.1μM.The 95% confidence interval of active PDE4 level (see Fig. 2) is [3.21, 3.46].For information about the experimental details we refer the reader to the supplementary material.

Resulting conservation relationships
We uncovered 14 conservation relationships within the modelled biological system under steady state conditions.Figure 3 illustrates which sets of proteins have to maintain a constant total concentration.As illustrated in Fig. 3, we found at least one conservation relationship for each protein which reflects the   1. main endogenous signalling cycles.Thus our expectations coming from the signalling character of our modelled system are verified.

Constraint flux sampling correctly predicts DMR response under receptor activation
We applied the CFS framework described above incorporating cAMP as well as PDE4 fold changes into flux constraints.DMR response measurements were not taken into consideration at this point, but left out for independent validation.Figure 4 depicts the distributions of those selected fluxes, which according to our CFS analysis are predicted to show a statistically significant shift under stimulation (FDR <1%, Wilcoxon signed rank test with Benjamini and Hochberg's FDR control of P-values for multiple testing (Hollander, 1999;Benjamini & Hochberg, 1995)).CFS predicts a significant change of RGS 14.A slight decrease of the receptor de-activation and increase of the GP subtype alpha-i de-activation via RGS14 can be expected under stimulation according to our simulation.The receptor de-activation is compensated by an increasing receptor recycling.This phenomenon of signal regulation by RGS14 and GRK6 is well described in the literature where both proteins are known as important signal regulators (Pierce et al., 2002;Dale & Rang, 2011;Berridge, 2014).The RGS family is involved into the extinction of GP-dependent signalling (Zhang and Mende, 2011), e.g.via receptor desensitization or endocytosis (Reiter & Lefkowitz, 2006).This causes the downregulation of GP alpha-i downstream and together with cAMP forms a positive feedback loop.More specifically, the inhibition of the cAMP inhibitor GP alpha-i subunit leads to an increase of the cAMP production.In addition, significant GP alpha-i de-activation causes a significant increase of GRK6 activation.
The boxplots clearly highlight that-besides cAMP-under stimulation increasing levels of AMP production and cAMP degradation are expected, which is in agreement with current literature (Pierce et al., 2002;Sunahara & Taussig, 2002;Strange, 2009;Berridge, 2014).CFS is able to correctly predict a significant positive wavelength shift, i.e.DMR response, under receptor stimulation which is in agreement with our experimental validation data (see Fig. 2).Hence, CFS allows for a qualitative check of our pathway model.

Knock-out simulations
To further check the hypothesized relevance of RGS14 for the observed DMR response we conducted an in silico knock-out simulation.This means we restricted all fluxes going through this molecule to zero while repeating our CFS.To investigate the effect of the different molecules on the DMR response, we performed knock-out simulations for all molecules except for GP alpha-s, AMP, AC, cAMP, PKA, the receptor and the ligand.We then ranked the molecules according to their statistical significance of the influence on DMR related fluxes.GP alpha-s is not considered for the knock-out simulation because no steady state solutions are possible when constraining the fluxes through these important signalling molecules to zero.AMP, AC, cAMP, PKA, the receptor and the ligand are not considered because these are characteristic molecules for signalling and removing these molecules is unphysiological.Table 2 shows the proteins with influence on the response under stimulation.
Altogether our simulations underline our findings from section 3.3.The highest impact was found for RGS14 and GP alpha-o, followed by PDE4.Furthermore, GRK6 has a significant influence.

Combining EFMs and CFS reveals important sub-networks and regulatory mechanisms
In a last step, we applied EFM analysis to the system with ligand stimulation (see section 2.4), resulting in 63 EFMs (see supplementary material).Notably, many of these EFMs represent similar biological mechanisms.We ranked all EFMs with respect to their predicted change under stimulation by the method described in section 2.7.Table 3 shows all EFMs with an FDR lower than 0.001 and a median fold change greater than 1.Interestingly enough, four of the most significant EFMs describe the GP alpha regulation via RGS14, and also the receptor regulation via GRK6 is among the most significant EFMs (Figs 5 and 6).This is in full agreement with our previous findings and provides a possible explanation for the relevance of these molecules.Further significant EFMs are related to PDE activation and cAMP/GEF production (in agreement with our experimental data). .First the inactive G protein complex consisting of the subunits alpha-i and beta/gamma binds to the active receptor and the bound GDP (Guanosine diphosphate) is replaced by GTP (Guanosine-5'-triphosphate) while the G protein dissociates from the receptor and splits into its subunits beta/gamma and alphai.Afterwards the activated alpha subunit is deactivated by replacing GTP with GDP-mediated by RGS14.In the last step, the deactivated GDP-bound alpha-subunit again associates with the beta/gamma subunit and forms the inactive G protein.

Discussion
In this article, we presented a comprehensive mathematical model of the M 2 receptor-dependent joint signalling and secondary messenger network.The motivation for our work comes from the pharmacological relevance of the M 2 receptor and the induced cellular responses.Whereas in principle the individual parts of our studied system are well described in the biological literature, to our knowledge there have been no attempts so far to combine these information into a mathematical model.

Fig. 1 .
Fig. 1.Cartoon of the M 2 receptor-dependent signalling and secondary messenger network in CHO-hM2 cells based on the known literature.The receptor is activated by a ligand (e.g.iperoxo) and induces the membrane-bound signalling cascade including G protein (G) activation and production of cAMP by adenylate cyclases (AC).Via cAMP the signal is transferred to the PKA-induced phosporylation cascade.Detailed reactions are suppressed for simplification.The detailed reaction system can be found in the supplementary material.

Fig. 2 .
Fig. 2. (a) DMR concentration-response curve of iperoxo modified from Schrage et al. (2013).The affinity between iperoxo and the receptor (pK D ) of the iperoxo-induced ligand binding curve for intact cells (obtained from Schrage et al. (2014)) is marked by the blue line.(b) Concentration effect curve of measured iperoxo-induced G protein alpha-s mediated cAMP accumulation with standard deviations and estimated confidence interval marked by the blue line.The inactivation of G protein alpha-i was induced via a pretreatment with pertussis toxin (PTX).The inactivation of the cAMP inhibiting G protein alpha-i allows for matching the measurements with their corresponding G protein alpha-s mediated network fluxes.cAMP accumulation in the absence of test compounds was set to 0% and maximum forskolin-induced binding was set to 100%.(c) Western Blot for the total amount of PDE4 (Pan-PDE4) and active PDE4 (pUCR1) under stimulation with 0.1μM iperoxo normalized against GAPDH.(d) Fold change for normalized active PDE4 (pUCR1).

Fig. 3 .
Fig. 3. Calculated conservation relationships: each column represents one conservation relationship and each row a protein.Red cells indicate proteins involved in a concentration relationship.The sum over all marked protein concentrations per column is constant.The inactive form of each protein is indicated by the subscript 'in'.

Fig. 4 .
Fig. 4. Predicted fluxes without stimulation (yellow/left) and under stimulation (green/right).(a) Boxplots illustrating the distribution of selected fluxes under different conditions using cAMP measurements.Yellow (left) boxes indicate the unstimulated and green (right) boxes the stimulated condition.Ligand induced G protein activation is shown for the alpha-s subtype here.Boxplots for all fluxes can be found in the supplementary material.(b) Overall response given as the sum of fluxes into tubulin and actin, see 2.14.Related median fold changes are shown in Table1.

Fig. 5 .
Fig.5.Elementary flux mode for G protein (G) regulation via RGS14.First the inactive G protein complex consisting of the subunits alpha-i and beta/gamma binds to the active receptor and the bound GDP (Guanosine diphosphate) is replaced by GTP (Guanosine-5'-triphosphate) while the G protein dissociates from the receptor and splits into its subunits beta/gamma and alphai.Afterwards the activated alpha subunit is deactivated by replacing GTP with GDP-mediated by RGS14.In the last step, the deactivated GDP-bound alpha-subunit again associates with the beta/gamma subunit and forms the inactive G protein.

Fig. 6 .
Fig. 6.Elementary flux mode for the GRK-mediated receptor inactivation via phosphorylation.The ligand-bound receptor gets phosphorylated by GRK.The phosporylated and hence inactive receptor is no longer able to mediate G protein activation.After ligand-dissociation the receptor gets de-phosphorylated and is now again able to mediate G protein activation.

Table 1
Median fold changes related to Fig.4

Table 2
Proteins ranked with respect to their predicted influence on the DMR response.The influence on the response was estimated by the median fold change expected by a knock-out simulation of each protein.The statistical significance of each simulated fold change is shown in terms of FDR.A high fold change implies a strong influence of the particular protein

Table 3
Significant EFMs (FDR < 1E − 6) ranked by their median predicted fold change induced by stimulation.Note that the first three fold changes are not computable, because without stimulation there is no flux through these EFMs.For the complete table please see supplementary material