Transfer Functions for Protein Signal Transduction: Application to a Model of Striatal Neural Plasticity

We present a novel formulation for biochemical reaction networks in the context of protein signal transduction. The model consists of input-output transfer functions, which are derived from differential equations, using stable equilibria. We select a set of “source” species, which are interpreted as input signals. Signals are transmitted to all other species in the system (the “target” species) with a specific delay and with a specific transmission strength. The delay is computed as the maximal reaction time until a stable equilibrium for the target species is reached, in the context of all other reactions in the system. The transmission strength is the concentration change of the target species. The computed input-output transfer functions can be stored in a matrix, fitted with parameters, and even recalled to build dynamical models on the basis of state changes. By separating the temporal and the magnitudinal domain we can greatly simplify the computational model, circumventing typical problems of complex dynamical systems. The transfer function transformation of biochemical reaction systems can be applied to mass-action kinetic models of signal transduction. The paper shows that this approach yields significant novel insights while remaining a fully testable and executable dynamical model for signal transduction. In particular we can deconstruct the complex system into local transfer functions between individual species. As an example, we examine modularity and signal integration using a published model of striatal neural plasticity. The modularizations that emerge correspond to a known biological distinction between calcium-dependent and cAMP-dependent pathways. Remarkably, we found that overall interconnectedness depends on the magnitude of inputs, with higher connectivity at low input concentrations and significant modularization at moderate to high input concentrations. This general result, which directly follows from the properties of individual transfer functions, contradicts notions of ubiquitous complexity by showing input-dependent signal transmission inactivation.


Introduction
Biochemical reaction systems are usually conceptualized as dynamical systems -systems that evolve in continuous time and may or may not receive additional input to the system. Mathematically, this can be expressed by sets of ordinary differential equations (ODE), such that rates of concentration changes correspond to mass-action kinetic parameters [1,2]. In this paper we use existing mass-action dynamical systems to propose an alternate or additional framework for modeling and interpretation of biochemical reaction systems. We provide an algebraic analysis of biochemical reaction systems as a matrix of concentrations for all species, given certain input concentrations. These concentrations correspond to steady-state amounts which are reached after a delay time, and the delay times can be measured by the system as well.
We use an arbitrary published model [3] as an example for a ODE dynamical model of biochemical reactions. The model simulates intracellular signal transduction from receptor binding to molecular targets in different cellular compartments, as an important component in the long-term regulation of protein expression implied in neural and synaptic plasticity. In striatal neurons, both a calcium-dependent pathway and a cAMP-dependent pathway are activated during the initiation of neural plasticity by NMDA/AMPA receptors and neuromodulator receptors such as dopamine D1 receptors [4][5][6]. Their effects and the integration of signaling on common targets such as kinases and phosphatases have been the subject of a number of computational models [7][8][9][10][11]. In particular, the role of the DARPP-32 protein in striatal neurons in determining the outcome of membrane signaling has been modeled by different groups, based on a common set of experimental data [3,12,13], cf. [14][15][16][17][18]. Many similar models [19] have been developed in the last 10-12 years in different areas of biology. Models with dozens or more of species have up to a 100 or more equations and are consequently complex and difficult to understand as continuous dynamical systems [20]. A transformation into a matrix-based formulation of input-output functions, even at the cost of a loss of fast dynamical modeling, promises considerable gain of insight and access to a different set of mathematical tools. Simple mass-action kinetic models may be criticized for disregarding the real complexity of spatio-temporal molecular interactions. Some alternatives use spatial grids and stochastic versions of biochemical reactions to capture this complexity [21,22]. However, certain variations, such as compartmental modeling with diffusion, altered kinetics for anchored proteins, or employing molecular kinetics as the basis for binding constants may be employed within the mass-action kinetic framework to achieve better correspondence with the biological reality. These variations can be directly transferred to the proposed model as well.
In our approach, we identify input nodes, and then precompute the outcomes for all internal species (target species) in response to biological meaningful ranges and combinations of inputs. This allows to analyze a biochemical reaction system under all possible input conditions. The analysis can be done for arbitrary ODE models [19], provided minimal requirements on conservation properties are realized (cf. section ''Methods'', [23]). The results are stored as vectors or matrices ('systemic protein signaling functions' (psfs)) and can be fitted with functional parameters. It is an important aspect of the model that computations are done systemically. In section ''Elementary Biochemical Reactions'', source-target interactions are first analyzed in isolation ('elementary psfs'). They all constitute hyperbolic saturation functions, therefore rate parameters can be uniformly translated into functional parameters for signal transmission strength. But in a systemic context, source-target interactions change because of additional influences on the species from other equations in the system (section ''Systemic PSF Analysis''). Therefore a fitted systemic psf from A to B is different from the elementary psf. The pre-computed, systemic psfs may be used to create state-change simulation models, i.e. discrete-time models, which can be compared with continuous ODE models (section ''Systemic Delay and State-change Dynamics''). What is significant and novel about our analysis is that we can extract systemic transfer functions from the complex system, and thereby dissect the system into parts. We can analyze the transmission properties of individual species, compare their minimum and maximum values, and the functional shape of their transmission strength. Specifically, we can show under which circumstances a link is functional, i.e. actually transmits information (section ''Computing Input-dependent Modularity'').
The analysis has a number of restrictions. An important restriction is that our model does not allow for analysis of fast interactions below the resolution of settling into steady-state. The requirement of conservation of mass guarantees that for each input concentrations will eventually settle to some equilibrium value, but due to the prevalence of feedback interactions, they may still produce transients or dampened oscillations. This means that fast fluctuations of input will not be adequately simulated using precomputed psf functions alone. It is then necessary to refer to the underlying ODE model. The model is most suitable for studying disease states, pharmacological interventions, genetic manipulations, miRNA interference, or any system conditions which fundamentally alter the presence or concentration of molecular species. These conditions may then be tested either in steady-state or with a sequence of sufficiently slow input constellations. A second restriction is that the model inherits parameter uncertainty from mass-action kinetic models. These parameters are derived from experimental measurements, but typically with a high degree of uncertainty [24]. Our analysis offers a clear distinction between elementary and systemic functional parameters and explains why experimental measurements are so highly dependent on the systemic context. In this paper, we have treated elementary parameters as given by the underlying model [3]. In principle, systemic functional interactions can be measured experimentally, and in the model each interaction can be adjusted separately. This may offer a novel theoretical approach towards finding adequate elementary rate parameters. A biochemical reaction system formulation for signal transduction contains two different types of reactions: 1. complex formation ½Az½B<½AB 2. enzymatic reactions ½Az½E<½AE?½AÃzE

System Definition
The system has concentrations for species A, B and E (A*, AE and AB can be calculated), a set of kinetic rate parameters k on and k off for the forward and backward binding reactions in complex formation, and k cat for the rate of enzymatic production. The equational structure, the kinetic parameters and the initial concentrations of the model [3] are reproduced in Tables 1, 2, 3, 4, with slight modifications: Equation 40 was added to 'close the loop' from AMP to ATP and thus provide for conservation of all molecules in the system. Conservation of mass is necessary in order for all species to reach equilibrium. It means that for any forward reaction there needs to be a reverse reaction, such that any species receives both input and output ('weakly reversible system', cf. [23]). This implies that pure loss reactions, like endocytosis or diffusion across the cell membrane (secretion) cannot adequately be modeled with this system, unless balancing reactions are added which make up for the loss. E.g., for endocytosis of a ligand-bound receptor, both the ligand and the receptor, possibly independently, have to be recycled, which means that bridging equations for receptorligand dissolution in the endocytosed state, and input rates for receptors and ligands have to be added. Secondly, the species PP1 and its interactions (4 equations) were left out, since they contain a complex which dissolves into its 3 components in one step: this would require a, fairly trivial, addition to the current psf system implementation. The system can be depicted as a bipartite graph with nodes for species and nodes for reactions.

PSF Analysis
In order to set up source-target functions, we need to select input nodes from the available species nodes. In this example, we used Da (dopamine as ligand for the D1 receptor) and Ca (extracellular calcium that diffuses through ion channels in the membrane). We use input concentrations over a specified range (e.g., between 60nM and 5mM for Da), sample over the range with e. g. n = 20 steps, and use the differential equation implementation of the system to calculate the output values for all species for each sampling step. Because of the conservation of molecules, all species reach steady-state after a sufficient period of time. We define steady-state pragmatically by relative change of less than 2% over 100s. We also use the established terminology of EC10, EC50, EC90 etc. to indicate 10, 50 or 90% of steady-state concentration value. Additionally, we calculate the delay in reaching steady-state. We store input-target concentration mappings in a vector (singleinput system), or a matrix (multiple-input system). We fit the vectors with hyperbolic or linear functions, using standard techniques in Matlab (fminsearch, [25]). In this way we derive parameters which can be analyzed and used instead of the explicit vectors. (In this paper, the fitting is only done for single-input systems, multiple-input systems require different techniques.).
All information on source-target transfer functions for the complete, complex signaling system ('systemic psf') can be stored in a static data structure. For each species, it contains its concentration range, and for each reaction, it contains the parameters of the functional fit. We gain the possibility to regard any species as source and any other species as target (they may be coupled by an arbitrary number of reactions) and obtain a systemic psf as the transfer function between them. This representation allows to analyze the complex signaling system by its parts, i.e. as a set of matrices or vectors, which is the main achievement relative to the ODE dynamical system. In addition, dynamical simulation with appropriate update times may be realized by the psf representation alone, i.e. the psf simulation is not in itself atemporal, but only discrete and fairly slow.
We visualize the data structure as a bipartite graph, and label it with the calculated numeric values. Each species node is labeled with its attainable concentration range given the input range.

Elementary Biochemical Reactions
We want to represent a biochemical reaction by a timeindependent signal transfer function, such that y~f (x) for two species x,y. We do this by designating a source species x and then calculating the steady-state value for another species, the target species y, for any value of x, given the differential equations for the biochemical reaction. For complex formation  We may now calculate the concentration values f (x)~y for a target species [AB] given a range of input values for x, e.g. the source species [A]. (Fig. 1A).
In this way we separate the calculation of the signal response magnitude, i. e. the steady-state concentration, from the calculation of the time until a steady-state value is reached, the delay. For different x, f (x) will be reached after a variable delay (Fig. 1B).
With some modification, the same transformation applies to enzymatic reactions. The kinetic rate parameters are Kd~k off kon (for <) and k cat (for ?).
Here it is required that the enzymatic reaction is reversible, i.e. a reaction is a reaction that reverses [A*]. ) The differential equations, with k cat2 for ½A Ã ?½A, are: Given concentrations for E 1 , A and kinetic rate parameters Kd, k cat and k cat2 , we may now derive a function with x as the source species ½E 1 and y as the target species ½AÃ (Fig. 1C).
In both cases the resulting curve can be fitted by a saturating hyperbolic function.
Here y min , the baseline concentration, is usually set to 0. If we choose [A] as the target of [E1], we get a negative slope psf.
We call this function the elementary protein signaling function or elementary psf. This function is somewhat related to a Hill equation [26,27]. A Hill equation is a function fitted to an experimental measurement to derive a dose-response relationship, comparable to the psf. The Hill equation allows to calculate a fractional concentration h for the target (e.g. a receptor-ligand complex) from the source concentration ½L, given Kd, and fitting a parameter n for the steepness of the curve.

h~½
L n Kdz½L n The concentration of the other compound of the complex is not used (assumed large), and the absolute magnitude of the target is not calculated. An equivalent for enzymatic reactions is not defined. The parameter n allows to measure the effect of competing binding reactions (n = 1 if none are present), which in our terminology translates into a systemic psf with multiple binding partners for a single target compound. Systemic psfs are a more general concept than Hill equations, but they relate to the same type of data, namely dose-response functions in steady-state.
We have seen that signal transmission strength is uniformly characterized by saturating hyperbolic functions. This means that it is highest for low x and diminishes as x increases (Fig. 1D). For instance, in Fig. 1D, a 100% signal increase leads to 100%, 18% or only 6% increase in the target depending on the source concentration. For enzymatic reactions, absolute concentration changes have different effects for sources and targets of a signaling interaction. Signal transmission strength depends on the absolute concentration of the source, the target concentration is irrelevant. This is an important observation, since protein signaling systems are subject to long-term regulation of concentrations. In the context of disease states or other sources of protein expression up2/downregulation, independence of transmission strength from target concentration may be an important conservative property.
Signal transmission is strongest if a source species is expressed at a low concentration. We need to bear in mind however, that reaction velocity operates inversely to signal transmission strength: a low source concentration means a slow reaction (Fig. 1B). A functioning signaling system would therefore have to use an intermediate range to maximize signal transmission within time constraints. Our analysis opens a new way of analysis for a signaling system: Optimization techniques could find a best source range for both time and signal transmission constraints.

Systemic PSF Analysis
A source-target psf can be derived for any pair of species in a complex biochemical reaction system. For a complex system, or set of equations, we define a set of input nodes, and compute the output values for each possible input configuration. The analysis gives us the output concentration range (notwithstanding transients in a dynamic context, s. below) for each species, as well as a (fitted) function, or matrix of input-output correspondences. A biochemical reaction will produce a different psf, when it is elementary or when it is embedded in a context, where the participants of the elementary reaction also participate in other reactions. This is true for both protein complex formation and enzymatic reactions. We therefore call source-target functions 'systemic psfs', when they are derived from the context of a specific signaling system.
We have provided this analysis for the example system. We show the concentration ranges and the signal transmission functions for the whole system [3] as a weighted dynamic graph for Da as a single input (Fig. 2). We label each species node with its concentration range, determine source and target species nodes for each reaction node, and provide fits for the systemic psf, the transfer function that characterizes each reaction. We see from (Fig. 2) that a number of systemic psfs can be fitted well with a linear function (y~mxzb), showing that systemic psfs sometimes consist only of a short section from a full mapping of concentration values. Also, many species have only small concentration ranges, which means they don't have much response to Da input.
It is an obvious advantage of the psf analysis that we are able to dissect the complex system and extract local properties, such as concentration ranges of individual species, and transfer functions for individual reactions under input stimulation. This allows to critically analyze a model, compare these properties with biological data, and adjust or improve the model in a detailed manner.
In Fig. 3, the concentration ranges for some target species are given. We see, for instance, that among DARPP-32 phosphorylation variants, pThr75 is always more abundant that pThr34, by an order of magnitude. This is an example of a high-level property, which could be related to biological data. As another example, we notice that the active receptor conformation (Da-D1R) remains below 160 nM even under stimulation with 1mM Da and more. With a D1R total concentration of 500 nM, we could adjust the ligand binding coefficient to produce more or less active receptors. Finally, the analysis shows a very low maximal PKAc level (12nM) in spite of a total PKA concentration of 1:2mM. In the original model [3], blind parameter adjustment has probably generated a very low level of PKAc in order to achieve high signal transmission for phosphorylation of the target species pThr34, which is experimentally required, but which could be achieved in other ways (e.g. PP2B) as well.
With our analysis, properties of individual species become apparent, and they can be compared to biological data, tested and adjusted on a localized basis. Even more interestingly, we could look for principles of 'rational system design', for instance question the transmission of a seven-fold increase of cAMP in the mM range to a maximal three-fold increase in the 109s of nM for PKAc, and analyze given biological systems from this perspective.
In addition to the concentration ranges, we also have access to the functional mapping between species in the model. The systemic psfs, like the elementary psfs, are stored as vectors, which are matched by functional parameters. The advantage of the psf analysis is that we can probe a complex system on a single reaction level because the influence of the cellular context is encoded in the systemic psf. Thus we can compare the elementary psf with its transformation as a systemic psf for individual reactions. Fig. 4A shows elementary and systemic psfs for G-protein activation of AC5 and the calcium-activated complex AC5Ca. We see that the systemic psfs are somewhat deflected, compared to the elementary psf, which is what we expect from the parallel activation of AC5Ca and AC5 by the same species. We may specify a desired psf using only functional parameters, and adjust elementary parameters to match the psf (Fig. 4B). A local change to the Kd binding coefficient between AC5Ca and GoaGTP allows a change in the systemic function. Since other systemic psfs may be affected by such a change -this can be detected by re-computing the weighted dynamic graph -more adjustments of elementary rate parameters may be indicated, possibly by an iterative process (cf. [28]).
Sampling for multiple inputs yields a transfer function matrix, which can be analyzed for dependence of the target concentration on each input separately. This can be done by standard matrix analysis such as principal component analysis (PCA). For our example, we show how species which are poised to integrate signals from two different sources do this under the numeric conditions (Fig. 5). For cAMP production (AC5), we find that AC5GoaGTP is dependent only on Da (Fig. 5A), AC5Ca is only dependent on Ca (Fig. 5B), while AC5CaGoaGTP is almost not activated at all (Fig. 5C). Even though a link of reactions (Ca-AC5Ca-AC5CaGoaGTP-cAMP) exists, signal integration of Ca and Da on AC5 fails because of the weak transmission from GoaGTP to AC5Ca. Signal integration between Ca and Da occurs for cAMP degradation by calcium-dependent calmodulin regulation of PDE1. PP2A with the two variants (calciumactivated) PP2Ac and (PKAc-activated) PP2Ap is another potential source of signal integration (Fig. 5D,E,F). The psf analysis shows when signal integration occurs (here: Da having influence on PP2Ac), and when this effect is negligible (here: Ca not having influence on PP2Ap). This may now be studied for correspondence with the biological situation. These results emphasize the necessity for numeric analysis of input-dependence, beyond the mere existence of links.

Systemic Delay and State-change Dynamics
We would like to be able to use systemic psfs with their simple and transparent mathematical structure for dynamical simulations. This allows direct experimental testing and fitting by time series measurements beyond dose-response relationships. In order to do this, we need to compute the systemic delays, i.e. the reaction time until a steady state is formed. Then we can build a state-change dynamical model from systemic psfs alone, using the appropriate delays for the input and the update of a system state.
Systemic delays depend on the absolute size of the signal and also the direction (increase/decrease) of signaling. Delays for species in the example system in response to input are shown in Table 5. For the computation of target concentrations, we only need a ratio such as kon koff~K d (binding coefficient) or kcat kcat2 for forward and backward enzymatic reactions. For the delays, the difference between k on and k off or k cat and k cat2 defines reaction times for synthesis and degradation. Therefore, delay computations are fairly complex, but the results are often within a fairly narrow range for each reaction (Table 5). For discrete state-change simulations we may use maximal delays for each species. From a biological perspective, this table provides an important test on the validity of the model. In many cases, systemic delays can be measured. For instance, the delay for PKAc at 150-250 s rather than 30-60 s, as measured in [29] (cf. [9]), seems large and may be an indication for a revision of the underlying parameters. From the theoretical perspective, this system seems to operate on separate time scales: 1-10 s, 150-300 s and 450-600 s. Such a separation of reactions by their characteristic delay times is interesting, since it could lead to simulation models with different discrete time scales. Here we may calculate psf values for fast species with 10 s time resolution, for intermediate species with 300 s time resolution, and for slow species with 600 s time resolution, i.e. system update time for state changes (Fig. 6A,B). It is an empirical question, whether separate time scales rather than a continuum of delay values will prove to be an organizing principle in protein signaling systems [30]. A general study, for instance, using models from the BioModels Database [19], might give answers to this question. Time scale separation may provide a conservative property of a signaling system against fluctuations of concentrations. If total concentrations in the system change, e.g. by protein expression up-or down regulation, miRNA interaction, or diffusional processes across compartments, the relevant interactions will continue within each time slice. Concerted regulation of protein expression levels may set a clock for the rapidity of signal transduction.
Systemic dynamics, in contrast to elementary reaction dynamics, need not follow a hyperbolic curve. If there are feedbacks in the system, the dynamics may contain transients, i.e. the concentration may be higher or lower before it settles into its steady-state value [31]. The dynamic response of target species to input are shown in Fig. 6A,B. For a species without a transient response, the actual value of a species at a shorter delay is always bounded by the steady-state value, and all possible concentrations in a continuous-time dynamical system are bounded by the psf concentration range [30]. However, if there are transients, a psfbased dynamical simulation will miss these transients and plot a simplified trajectory. This means that results from a psf analysis with slow inputs cannot be extrapolated to much faster input dynamics -in contrast to continuous-timed dynamical systems where arbitrary time units can be chosen. This restriction may capture a biological reality: steady-state behavior provides the framework and may operate according to rules and principles which are separate from the effects of short term fluctuations.
The psf system allows to generate a dynamical system as a sequence of states defined by fluctuations of input. Fig. 6C,D shows an overlay of a differential equation simulation and psf state change simulation for a sequence of inputs with 10 s duration. Accordingly, the psf approximation is excellent for all species with a delay time of v10s. If we plot psf values for slow species at intervals corresponding to their maximal delays, we may linearly interpolate between points, and in this case achieve a good approximation of the continuous-time model.
A psf-based dynamical system is an important tool in order to generate a time-series simulation from a calculated model system for comparison with experimental data. The psf system utilizes parameters which are uniform and have linear error ranges (cf. Fig. 1A,C), and therefore should improve interaction of the model with the experimental reality. The psf model will also allow to predict the optimal stimulation times for different inputs such that responses can be measured in steady-state.

Computing Input-dependent Modularity
Since we are able to define signal transmission capacity, we have a tool to investigate modularity of a signaling system. As species saturate or return to basal levels, they act as inactive links, i.e. they are stuck at the same concentration value, and cannot transmit further increases or decreases of inputs. We hypothesize that this effect is actually important in many protein signaling systems. We may define an inactive connection as a species node which has only limited (e.g. v10%) signal transmission capacity. The interconnectedness of a system is then proportional to the number of inactive connections, and a module is a part of the system with few or no active connections to the rest of the system.
In the following we discuss the activation/inactivation of links with respect to input increases. I.e. given a certain level of extracellular signaling, what happens if this level is raised and then kept at the higher level for some time, sufficient for the system to settle into a new steady state? Which species nodes will respond to the increase and transmit it to downstream targets, and which species will become saturated and only respond with their saturated value? It is also clear that species which have become saturated (inactive) will not respond to fast extracellular signals anymore. This analysis is therefore useful both for the steady-state context and for understanding fast input fluctuations. There is an input level for each target species, where the species ceases to be responsive to further input increases. If that input level has been reached, the species can be considered to have become an inactive connection, i.e. a node which does not transmit signals. We may define systemic psfs for 'input-target psfs' from the input to any target species (Fig. 7). We notice how number of steps in the computation of a species concentration relates to a lower cut-off  value for signal transmission (e.g. DaD1R, AC5GoaGTP, cAMP vs. PKAc, pThr34, pThr75). In other words, earlier steps in the sequence saturate at a higher input level than later steps. Two parallel targets (pThr34, PP2Ap) of an intermediate step (PKAc) may saturate at very different input levels (3 3mM for pThr34 vs. 1:5mM for PP2Ap). This mechanism demonstrates the effect of a sequence of saturating functions, and constitutes a general principle in the construction of a signal transduction system. The model allows to study whether a node responds to a specific input with any change of steady-state. In Fig. 8 and figures S1, S2, we show modularization of the system under various Da and Ca input conditions. With Da input and high Ca ( Figure 8B), species which are proximal to Da input, the receptor-ligand complex DaD1R, the signaling complex through G proteins and adenylyl cyclase AC5, as well as cAMP, are most responsive to Da over a large range of input. Species in the 'integration zone' between Da and Ca, such as DARPP-32, PP2A cease responding to Da increases at lower levels and become inactive links at higher levels. Species in the system with no significant change in concentration at any input level are for the most part proximal to Ca instead of Da, and thus highlight modularity among pathways. In this case, we see that Da inputs are transmitted to distant targets only up to an intermediate range and that there are a significant number of species which do not react to Da at all. Above that range, even though closely coupled targets still respond to the input, the signal increase is not registered beyond cAMP production and synthesis. With Ca low ( Figure S1) there is widespread responsivity to Da up to 1-2mM and only a few of the calcium-related species do not respond at all. The Da signal is therefore able to influence the calcium-related pathway, provided calcium is low.
When Ca is the input (Fig. 8A), the calcium-responsive proteins like calmodulin, CaMKII, calcium-activated PP2B (calcineurin) transmit signals, while the GPCR pathway remains almost completely unresponsive. There is some signal integration with calmodulin-activated PDE1 for cAMP. Other than that, we see that PP2Ac and the pThr34 variant of DARPP-32 respond strongest to calcium while PP2Ap and the pThr75 variant of DARPP-32 have less or no responsiveness to calcium. With Da low, as in Figure S2, again, most of the GPCR-related species do not respond to Ca at all, or cease responding at 10-25% of maximal Ca (1-4mM). However Ca-related species like calmodulin, CamkII, PP2B remain responsive. The difference between the high and low Da condition is small. Here the few existing links from Ca to GPCR (such as Ca regulation of AC5) are not strong enough to influence the GPCR pathway significantly in any condition. In this case, it is not just the saturation that matters, the input signal has limited reach in influencing distant species in general.
We are analyzing a biological system with two 'pathways', the cAMP and the calcium pathway, which are cross-linked in a convergence zone of species which are influenced by both pathways. It is remarkable how clearly three different modules appear: the GPCR/cAMP pathway, the Ca pathway, and the signal integration zone.
There is also a general observation to be made about signal transmission in a protein signaling system: Signal integration is strongest when inputs are low. This is a direct consequence of the effect of coupling saturating nodes. It means that there are few saturating species in the system which impose modularity, and signals spread further. Widespread interconnectedness is only possible at low input levels. Because activation curves for biochemical reactions are mostly uniform hyperbolic (saturating), a stronger modularization with many inactive links results from higher input levels. This may also have implications on transient responses. For instance, phosphatase and kinase response often differs with kinases being saturated only at higher input levels. If and where that is the case, we may observe transient responses that only reflect kinase activity, since phosphatases are saturated and do not participate in signaling, they do not add or subtract and in this way obscure the kinase signal. This shows that the steady-state input-response system may work well as the framework wherein fluctuating signaling operates.

Discussion
Continuous dynamical systems -systems that use change over time as a system primitive -are notoriously difficult to analyze and may not be the best choice of a tool for signal transduction systems of moderate or large size. Steady-state matrix computations are simple and fast, scale well to very large sizes, and offer multiple opportunities for analysis. By calculating transfer functions from a systemic dynamical model, we also gain the opportunity to extract and analyze parts of the system. This may help in creating re-usable system parts. This paper demonstrated a transformation of a mass-action kinetic biochemical reaction model implemented by a set of differential equations into an input-response transfer function model. The transformation is done by calculating steadystate concentrations for each species in response to a range of input values, and then analyzing the resulting vectors (matrices) as the basis of a transfer function ('psf'). For small, toy-like networks of few components, such dose-response relationships have been investigated by [32]. They analyse kinetic models in the same way, however, they don't make a distinction between parameters in elementary interactions, and the actual parameters in a systemic context. They also do not address the question of temporal embedding of the dose-response relations. In our approach, we use dose-response functions, similar to Hill functions, but we are extending the concept. By using rectangular signals, i.e. constant signaling levels, we calculate the response as the steady-state value. In addition, we calculate the time to steady-state from the underlying dynamical model. Only because of this additional computation can one attempt to create a discrete dynamical model -something that is not within the purview of a Hill equation model. Hill equations attempt to fit or create systemic parameters, which are different from elementary parameters, i.e. they do recognize the dependence of the transfer function on the systemic context. The psf model is an approach of making Hill equations (dose-response relations) work in large-scale modeling.
There are numerous attempts to simplify dynamical models in the temporal domain by creating hybrid models (e.g. [33]). Sometimes slow interactions are regarded as constant, and only fast interactions are dynamically modeled. Sometimes, fast interactions are replaced by time-independent values and only slow interactions are dynamically modeled (many examples, for instance [34]). There are also attempts to replace an ODE model by a delay-differential equation (DDE) model, i.e. to compute and use explicit time delays and eliminate many intermediate species, simplifying the size of the model [35].
In our approach, time and concentration are regarded as separate, which makes it different from any hybrid approach. The main restriction is the assumption of a constant signaling level over periods of time sufficient to induce steady-state. The approach is therefore best suited to check the limiting conditions of a dynamical model, e.g. in drug development applications, where multiple dose-response relations derived from the model can be crosschecked and used to locally improve the model. However, the simple, atemporal transfer functions can also be applied directly: (a) in experimental settings where signal duration can be controlled, and (b) in physiological settings, when fluctuations occur around different mean values, and we model the step changes for the mean extracellular signaling level, but not the short-term fluctuations. It could be shown that psfs are sufficient to create discrete time dynamical models, with certain restrictions on fast dynamical inputs below the time resolution of the system. The primary focus of the analysis was on single-input systems, where signaling functions can be matched by parameterized hyperbolic functions. We showed how the analysis can be extended to multiple-input systems by computing and analyzing psf transfer matrices. The psf functions allow to cut through the complexity of the model and examine interactions in a localized way. If continuous dynamical models are being used with the goal of adequately simulating cellular processes, this kind of analysis is an indispensable tool to check for the consequences of the modeling assumptions in terms of dose-response relationships. The analytical tools of the psf approach may also be used to systematically investigate the effect of localized changes to the system [28,36], and to offer local corrections to the complex system in a transparent way.
The extraction of input-target psfs with characteristic hyperbolic saturating properties allows to determine input level dependent inactivation of a species node. Accordingly, we can define the limits of signal transmission by the distribution of inactive nodes. In the example system, a strong distinction of calcium-and cAMPdependent pathways and a signal integration zone were revealed. We could see that at low levels of input, widespread interactions are possible, while at higher levels of input, many species enter into a state of a constant function value and become inactive links. This corresponds to biological results and expectations, and provides a foundation for the concept of pathways in signal transduction. For transient responses, the steady-state input-response system provides the boundary values to which the system eventually settles. This is especially interesting if we have saturating response levels, which are unresponsive to further signaling (e.g. phosphatases), and allow transient signals (e.g. kinases) to emerge.
Signal transduction may also be analyzed from the perspective of rational system design. Such work is still in its infancy. We may for instance investigate the effect of negative feedback links on concentration ranges and times to steady state. Another question would be the optimization of a signal transduction system for the trade-off between speed and efficacy of signal transmission. With this choice of model, many new questions can be raised, and old problems like parameter dependency, modularity or signal integration can be addressed in a novel way. Figure S1 Graph representation of signal transmission. Input concentration for Da is 20nM … 5mM, for Ca is 100nM. There is widespread interconnectedness for low inputs. (EPS) Figure S2 Graph representation of signal transmission. Input concentration for Da is 100nM, for Ca is 100nM … 10mM. The Da/GPCR pathway is clearly separated from Ca inputs. (EPS)