Reduced modeling of signal transduction – a modular approach
 Markus Koschorreck†^{1}Email author,
 Holger Conzelmann†^{1},
 Sybille Ebert^{1},
 Michael Ederer^{1} and
 Ernst Dieter Gilles^{1}
DOI: 10.1186/147121058336
© Koschorreck et al; licensee BioMed Central Ltd. 2007
Received: 19 March 2007
Accepted: 13 September 2007
Published: 13 September 2007
Abstract
Background
Combinatorial complexity is a challenging problem in detailed and mechanistic mathematical modeling of signal transduction. This subject has been discussed intensively and a lot of progress has been made within the last few years. A software tool (BioNetGen) was developed which allows an automatic rulebased setup of mechanistic model equations. In many cases these models can be reduced by an exact domainoriented lumping technique. However, the resulting models can still consist of a very large number of differential equations.
Results
We introduce a new reduction technique, which allows building modularized and highly reduced models. Compared to existing approaches further reduction of signal transduction networks is possible. The method also provides a new modularization criterion, which allows to dissect the model into smaller modules that are called layers and can be modeled independently. Hallmarks of the approach are conservation relations within each layer and connection of layers by signal flows instead of mass flows. The reduced model can be formulated directly without previous generation of detailed model equations. It can be understood and interpreted intuitively, as model variables are macroscopic quantities that are converted by rates following simple kinetics. The proposed technique is applicable without using complex mathematical tools and even without detailed knowledge of the mathematical background. However, we provide a detailed mathematical analysis to show performance and limitations of the method. For physiologically relevant parameter domains the transient as well as the stationary errors caused by the reduction are negligible.
Conclusion
The new layer based reduced modeling method allows building modularized and strongly reduced models of signal transduction networks. Reduced model equations can be directly formulated and are intuitively interpretable. Additionally, the method provides very good approximations especially for macroscopic variables. It can be combined with existing reduction methods without any difficulties.
Background
Modeling of signaling pathways
Systems biology aims at a holistic understanding of cellular processes. Mathematical models that integrate the current state of knowledge are analyzed to understand system properties that are not apparent from the characteristics of their components.
Different approaches exist to model and analyze signal transduction systems. Qualitative modeling uses solely structural information about the network and performs no quantitative statements. Examples for qualitative modeling techniques applied to signal transduction pathways are Petri nets [1], interaction graphs and Boolean networks [2]. Qualitative models can be analyzed by computational studies and structural analysis. As an example, independent signaling paths and feedback circuits can be obtained from an interaction graph using logical elementary modes [2]. The concept of Tinvariants in Petri nets allows to identify self contained subnets that are active under a given input situation [1]. With both techniques it is possible to decompose a signaling network into smaller functional units. Qualitative analysis is especially helpful if structural information but relatively little quantitative information about the system is known. Quantitative modeling explicitly considers the quantitative nature of physical systems. Processes on the molecular level exhibit stochastic behavior. When only a low number of a specific molecule is present, stochastic modeling techniques should be used. The chemical master equations [3] describe the stochastic dynamics of chemical systems in a quantitative way. However, they are difficult to analyze and simulate. If concentrations are high enough, the system dynamics shows deterministic behavior. Thus, in many cases a simplified approach neglecting the stochastic nature of the dynamics can be used.
Inside the cell, signals propagate in time and space. Therefore, partial differential equations should be used to describe them exactly. This continuous spatial behavior can often be modeled by assuming fast distribution of all species inside a cellular compartment and transport reactions between the compartments. This leads to model equations in the form of ordinary differential equations (ODEs). If processes are modeled on a molecular level mass action kinetics are a good description of the chemical processes and are frequently used. This view is the basis for the work presented here and also used by many others [4–13].
The more information about the system is available, the more powerful quantitative modeling techniques become, as prediction accuracy of models grows.
Combinatorial variety
Modeling with ODEs is challenged by biological reality. In signal transduction, association and modification of a relatively small number of different molecules usually give rise to an enormous amount of possible protein complexes [14]. This problem of combinatorial complexity has been tackled in different ways. In 1998, a stochastic modeling tool called STOCHSIM was developed to circumvent the problem of combinatorial complexity in bacterial chemotaxis [15, 16]. This algorithm is especially suited for systems with a very low number of molecules where stochastic effects may play an important role. Up to now there is research on stochastic simulations of combinatorial systems.
Many deterministic models of large signaling networks that have been published within the last years neglect the combinatorial variety of protein complexes. Their focus is on small subsets of the occurring reactions and complexes [4–11]. It was shown that such reduced model structures may provide a good approximation to a model accounting for all feasible protein complexes and reactions [17]. However, the main difficulty is to decide which reactions and complexes can be neglected and which are essential. The suited reduced model structure highly depends on the kinetic parameters of the reaction network [17]. Even in very simple examples an apparently reasonable assumption may lead to significant approximation errors [18].
An alternative approach was suggested by Blinov et al. [19], who introduced the software tool BioNetGen. This program translates a rulebased model formulation into a complete ordinary diffierential equation (ODE) model accounting for all feasible reactions and species. Additionally, a graphical rulebased representation of signaling networks has been introduced [20, 21]. BioNetGen was used to generate complete mechanistic models of early events in Fcε RI [12] and EGF signaling [13]. In both cases only a very limited number of components and binding domains are included. Nevertheless, more than three hundred ODEs are required in the case of Fcε RI signaling and over one thousand for EGF signaling. A detailed model of a complete signaling cascade easily can grow to millions or even billions of equations as shown below for the insulin signaling pathway. Recently, a novel approach to tackle this combinatorial explosion of model equations has been presented [22], formally extended and generalized [18]. The approach bases on the view that the fundamental elements of signal transduction are domains instead of molecular species [23]. The conventional mechanistic description of all feasible multiprotein complexes is replaced by a macroscopic one, where the focus is on the state of domains like levels of occupancy or degrees of phosphorylation. For many realistic systems this method allows a very strong reduction of the complete combinatorial models. The approach proposed by Borisov et al. [22] identifies reduction opportunities based on site dependencies that are determined from rules, not from the expanded network. The main deficiency of that method is that it considers only modularity within a single multisite molecule and does not consider modularity of the reaction network as a whole.
The exact domainoriented lumping technique proposed by Conzelmann et al. [18] guarantees an accurate description of macroscopic quantities. It is also applicable to whole reaction networks. However, a prerequisite for using this method is availability of the detailed model equations. The resulting models most often still consist of a very large number of ODEs which makes them difficult to handle and to analyze. In the following the insulin signaling pathway will serve as a realistic example to exemplify combinatorial variety. Afterwards, we introduce a new approximative reduced modeling method and show its efficiency by mathematical analysis and applying it to insulin signaling models of varying levels of detail.
Insulin signaling and combinatorial variety
The insulin signaling system is of high medical interest and therefore well studied [24–27]. Insulin regulates cellular glucose uptake [24, 28] and has a major impact on metabolism [26, 29]. It promotes synthesis and storage of glycogen, proteins and lipids and negatively regulates their degradation. Additionally, it negatively regulates secretion of sugars, amino acids and fatty acids [25]. Insulin is also involved in gene expression [30], cell survival and differentiation.
Defects in the insulin signaling system give rise to insulin resistance, obesity and type II diabetes mellitus [31–33]. There are intense efforts to improve therapies of these maladies [34–37].
In a simplified form the insulin signaling system will serve as demonstration object for combinatorial complexity. The insulin receptor is a transmembrane protein that is constitutively dimerized [38].
Results and discussion
Introduction of the layer based approach
The two most important mathematical tools to tackle the enormous complexity of models describing biological reaction networks are model reduction and modularization techniques. Here, we introduce a new systematic approach which allows to create considerably reduced models of signaling networks and also provides a new modularization criterion. This new modularization criterion suggests to separate molecular processes (bindings and posttranslational modifications) depending on the types of interactions between them.
The basic idea is that three different types of interactions between two processes in signal transduction networks exist. The first interaction type is called allornone interaction. In this case, there exists a causal relationship between the two processes, which means that one process must occur before the second process. Most frequently, this kind of interaction is between binding site phosphorylation and effector binding. These two processes can be divided into four molecular events (phosphorylation, dephosphorylation, binding, dissociation). The effector can only bind if the binding site is phosphorylated and dephosphorylation is only possible in the absence of effector. This means that phosphorylation is required for binding and dissociation is required for dephosphorylation.
The second interaction type is called graded interaction. There, two processes influence each other mutually, as it is the case for ligand binding and autophosphorylation of a receptor. This means that kinetic parameters for one process are influenced by the other process. Note that this allows unidirectional interactions, where the first process influences the second, while the second does not influence the first, and bidirectional interactions, where both processes influence each other.
According to our new modularization criterion, no graded interactions are allowed between different modules. Processes of different modules, which we call layers, interact only via allornone interactions. Interestingly, the so defined layers only exchange information. The number of molecules in each layer stays constant when no synthesis or degradation is considered. This is a main difference compared to metabolic pathways where mass flows occur and represents the characteristic signal flow within a signal transduction network. Considering the fact that in signaling an extracellular signaling molecule causes transmission of a signal to the nucleus of a cell without passing the cell membrane, this difference is most obvious. The signals that are exchanged between layers correspond to a very restricted number of macroscopic variables like levels of occupancy or phosphorylation as they are used by Borisov et al. [22] and Conzelmann et al. [18].
We also contribute a new aspect to the discussion about modularization of biological networks and the optimal criterion for modularization (see also [42–45]). Modularization in our new approach is tightly linked to model reduction. In fact, modularization sets the preconditions to directly postulate the reduced model, without preceding generation of detailed model equations.
In detailed modeling binding or modification events are represented by a huge number of reactions, since the involved proteins can exist in a high number of feasible configurations. The sum of a certain subset of these reaction rates defines a gross reaction rate of binding or modification. The sum of a certain subset of species corresponds to macroscopic quantities as degrees of phosphorylation or occupancy. In the following we show that in most frequent biological scenarios gross rate kinetics can be formulated using macroscopic variables and have a quite simple structure. Mathematical analysis of detailed and reduced models shows that the dynamics of essential macroscopic quantities is highly preserved in most cases. We give qualitative and also some quantitative information about the approximation quality for varying kinetic parameters. A basic advantageous feature of models that are created with the layer based approach presented here is intuitive interpretability. Additionally, preceding generation of detailed model equations is not necessary, as reduced models can be built directly by an intuitive procedure, for which a step by step procedure is given. There exists also a mathematical formalism to derive the reduced model equations. Thus, one can access model generation intuitively or by a mathematical formalism. Since both approaches are equivalent, understanding of the mathematical part is not necessary to create reduced models.
The approach combines qualitative and quantitative system descriptions. The first, qualitative step identifies processes and their interactions to define modules. Inside these modules, processes are described highly reduced using quantitative techniques. All methods of quantitative system analysis then can be applied to the model.
A general problem in modeling of signal transduction is availability of experimental data. Especially kinetic data is difficult to measure. Obviously, the layer based approach has also to face this problem. However, the same kinetic parameters can be taken for the macroscopic quantities that are states of the model as in detailed mechanistic modeling. Therefore, the problem does not become worse when using the reduced modeling technique.
Necessary equations for modeling the insulin signaling system
Scenario  Detailed modeling  Exact lumping  Layer based reduction  Combination 

1  145 156 468  145 156 468  214  214 
2  145 156 468  212  214  56 
The approach does not have to be used rigorously throughout the model. Subsystems can also be modeled in the detailed formalism, which simplifies integration of existing models.
A simple example system
D_{0} and E_{0} are the total concentrations of receptor D and effector E. We denote all receptor states of the detailed model by a leading D, while later on those of the reduced models will be denoted by a leading R. In the following we will always assume that detailed kinetic models, which build the base for our consideration of gross rates, are formulated using the law of mass action, which is a frequently used assumption.
Definitions
a) Molecules and complexes
In order to provide a complete and general description of our method, we introduce a formal nomenclature. Due to its generality this nomenclature might appear cumbersome. A simplified nomenclature can be used in most cases because many examples only comprise a small subset of all formally possible cases. In our examples we consistently use a simplified denotation.
Consider a general signaling protein R, that provides a number of distinct sites. In the general case each of these sites i can be modified from a state m_{ i, a }(e.g. not modified) to other states m_{ i, b }(e.g. phosphorylated). We denote the molecule with a certain configuration of domain modifications as D_{ R }[m_{1}...m_{ n }]. We distinguish between the molecule D_{ R }[m_{1}...m_{ n }] and its concentration D_{ R }m_{1}...m_{ n }. In a second step we also consider binding of other molecules to a certain binding site of R. Such a molecule might be D_{ E }[m_{1}...m_{ k }], which also provides a number of sites. One of them can bind to R. The result of this association is a complex consisting of R and E. This shall be denoted as {D_{ R }[m_{1}...E...m_{ n }], D_{ E }[m_{1}...R...m_{ k }]}.
The general rule for representation of complexes is as follows. All molecules within a complex are listed comma separated within curly brackets. On each occupied binding site the name of the binding partner is indicated. A complex of the molecule R binding E and F would be denoted as {D_{ R }[m_{1}...E...F...m_{ n }], D_{ E }[m_{1}...R...m_{ k }], D_{ F }[m_{1}...R...m_{ q }]}. If the same molecule occurs more often than once within the complex, indices have to be used. A graphoriented molecule representation is suited to solve the problem of correct denotation in difficult cases.
Since our method often works with lumped states comprising a number of molecular species, we introduce the symbol 'X', which is a replacement character for each possible state of a binding site. An 'X' within the site definitions of the molecules indicates all possible modifications on one site. A sequence of three dots following 'X' or before 'X' indicates that there may be additional sites and abbreviates a sequence of 'X'. A sequence of dots between site configurations other than 'X' indicates other sites with distinct configurations (not 'X'). D_{ R }[X...P] for example is a lumped state describing all molecules R with a phosphorylated site, no matter in which state all other sites are. Finally, one has to distinguish between D_{ R }[X...p]with a lowercase letter and D_{ R }[X...P] with a capital letter. While D_{ R }[X...P] comprises completely all species with a phosphorylated site, D_{ R }[X...p] only includes those phosphorylated species which do not have bound any other molecule at this site. Thus, if a phosphorylated site cannot bind any other molecule it necessarily is D_{ R }[X...p] = D_{ R }[X...P].
b) Rules and reactions
According to Blinov et al. [19], rules are a possibility to condense the description of a reaction system. They correspond to macroscopic chemical reaction equations. Each rule represents a set of chemical reaction equations that have common properties. Using the nomenclature defined above such a reaction rule describing a modification can be written as
D_{ R }[X...0....X] ⇋ D_{ R }[X....p...X], (2)
and a binding reaction as
D_{ R }[X...p....X] + D_{ E }[X...0...X] ⇋ {D_{ R }[X...E....X], D_{ E }[X...R...X]}. (3)
These rules shall be interpreted as a set of elementary reactions, which all can be modeled using the mass action law. If we consider the simple example defined above a possible reaction rule is
D[0, X] + L ⇋ D[L, X], (4)
are equivalent, since both describe the reactions d_{2} and d_{4} in Equation 1 and Figure 3. Reaction rules are a short way to formulate the complete and detailed stoichiometry of a reaction network. A short list of reaction rules may correspond to a very long list of possible complexes and elementary reaction steps and thus to large number of model equations.
Our goal is to formulate reduced order models in terms of macroscopic chemical species (denoted by a leading R) and macroscopic gross reactions (denoted by r) converting them. The phosphorylation of a binding site may be described in a reduced manner by the gross reaction
R_{ R }X..0...X ⇋ R_{ R }X...P...X, (6)
with the rate
r = r(R_{ R }X...0...X, R_{ R }X...P...X). (7)
The corresponding submodel consists of the two ODEs for the macroscopic variables R_{ R }X...0...X and R_{ R }X...P...X and thus is much smaller than the model defined by the corresponding reaction rule. The crucial point in model reduction is the derivation of suited gross reaction rates r_{ j }in dependence of the concentrations of the macroscopic species. In this work we suggest a method for deriving approximative expressions for gross reaction rates.
c) Different types of interactions
As introduced earlier, there exist three structurally different types of interactions. Now we discuss the three types by means of our example (Figure 3). In this example we consider three processes, namely ligand binding/dissociation, receptor phosphorylation/dephosphorylation and effector binding/dissociation. First, we focus on ligand binding and phosphorylation, which we describe in a separate reaction network. These processes form a reaction cycle that is shown in Figure 2A. By choosing the kinetic parameters of the occurring four reactions one determines if and how the two processes interact. Note that a process, e.g. phosphorylation, includes two molecular events, e.g. phosphorylation or dephosphorylation. The same holds for a binding process which consists of the events binding and dissociation.
We start considering the special case of noninteracting processes. If the kinetic parameters describing the reactions d_{1} and d_{3} are identical and the parameters of d_{2} and d_{4} are also identical, the two processes – in this case phosphorylation and binding of L – are completely independent. Hence, there is no interaction between the two processes. If this condition is not fulfilled, the two processes interact. Characteristics of this interaction can vary in a wide range depending on the corresponding parameters. Due to this great variability we call these interactions graded interactions, which form the first important type of interaction. A hallmark of a graded interaction between two processes is that both processes can proceed independently. However, they can influence their kinetic properties mutually.
A structurally different type of interaction is given if there is a causal relationship between the two processes, which means that one processes must occur before the second process. In many cases this probably will only be an approximation to the real interaction. However, it allows to strongly simplify the model. One can assume that receptor phosphorylation only occurs if ligand is bound to the receptor and the species DOp never occurs. The reaction network, which is shown in Figure 2B, now is reduced to a reaction chain. In this case all receptors follow the remaining reaction path and none the deleted one. Hence, we call this type of interaction allornone interaction. In most examples the assumption that phosphorylation and ligand binding interact in this way will only be a rough approximation. However, if we consider phosphorylation and effector binding the importance of allornone interactions becomes apparent. This case is depicted in Figure 2C, where DOOO is the unphosphorylated receptor without bound effector E. Phosphorylation is indicated by a p on the second position, binding of E by an E on third position. Note, that the species DOOE represents species without phosphorylation but with bound E. However, phosphorylation is an essential precondition for effector binding and effector binding prevents dephosphorylation due to steric reasons. Therefore, as indicated in Figure 2C, the species DOOE and the rates related to it do not exist in this case. Interactions, where both processes can only occur under certain preconditions are denoted as allornone interactions. A hallmark of allornone interactions is that reaction cycles degenerate to reaction chains due to nonexisting species. Most of the occurring allornone interactions are interactions between phosphorylations and associated bindings. The reaction scheme in this case can be simplified by omitting nonexisting components (Figure 2D) and the notation can be simplified as indicated in Figure 2E and introduced before. This fully simplified notation will be used from now on. Note that allornone interactions can be interpreted as a limit case of graded interactions, where certain parameters (e.g. the rate constants for d_{2} and d_{3} in Figure 2B) are set to zero and therefore reaction cycles degenerate to reaction chains. Analogously noninteracting processes realize another limit case of graded interactions, where there are equality constraints on parameters. However, we define both limit cases not to be graded interactions, since they allow model reduction.
d) Modularization: layers
We define that all binding and modification events coupled by graded interactions form a module which we call a layer. Note that the definition of layers is according to interactions of processes and not according to molecules. Roughly speaking, a layer contains processes, not molecules. However, layers are often dominated by the reactions of a single molecule.
An important feature of this modularization is that it only depends on the definition of interactions. The introduction of additional graded interactions between processes of different layers leads to disruption of the modular structure. In most cases these layers now form one larger common layer.
The terminus 'layer' results from the typical structure of the modules. As described, usually binding of a signaling molecule and its modifications (usually phosphorylation) form a layer. There is also the possibility that a layer comprises only one process, typically a binding process, as can be seen for Ebinding in Figure 4A and Fbinding in Figure 4B. One can imagine the successive layers of signaling that are traversed while the signal propagates through the pathway.
The layers are interconnected by signal flows. This means that information about macroscopic quantities of a layer is exchanged with other layers (signal flow). No mass flows cross layer boundaries. This means, that no reaction exists that transports substance from one layer to the other. Therefore, in the absence of protein synthesis or degradation, the sum of concentrations of all species for each molecule remains constant within the layers. Within layers there are mass flows defined by reaction equations and corresponding rates as in detailed modeling.
Altogether, layer based reduced modeling and modularization results in a highly structured reduced model that is characterized by mass conservation within layers and signal flow between layers.
Description of the reduction method
First, we consider the more academic case that all processes within a reaction network either do not interact with each other or provide an allornone interaction. Afterwards we will consider general networks also including graded interactions. General considerations are illustrated using the simple example defined above.
a) Networks without graded interactions
General considerations: We assume that within a reaction network all occurring processes do not interact, except domain phosphorylation and subsequent effector bindings that are characterized by allornone interactions. Phosphorylation usually can be considered as an essential precondition for binding an effector protein. The reverse modification (dephosphorylation) is prevented by binding of this molecule. Both conditions have to hold, otherwise the interaction between phosphorylation and binding will be graded.
Herein, I_{ gross }is the set of all reactions describing the binding of R and E. Observe that the lumped states of the model are capital letter species R_{ R }[X...P...X]. They occur in gross reactions of phosphorylation events. Lowercase species R_{ R }[X...p...X] occur in gross reactions of binding events. No ODE is required for lowercase species, as they are defined by conservation relations (Equation 9) between the sum of all phosphorylated binding sites R_{ R }X...P...X and the sum of all occupied binding sites {R_{ R }X...E...X, R_{ E }X...R...X}.
R_{ R }X...p...X = R_{ R }X...P...X  {R_{ R }X...E...X, R_{ E }X...R...X} (9)
Obviously, the same considerations can be made concerning modification processes. So, it is only necessary to balance capital letter species.
Gross rates of bindings as well as modifications can be formulated using the mass action formalism, as the parameters of all single elementary processes are equal. This is also discussed by Borisov et al. [22] from another point of view. The resulting model in this special case also corresponds to what the application of the domainoriented approach provides [18]. Hence, all approaches are exact and consistent in this special case.
Example: The example (Figure 3) comprises the three processes ligand binding, binding site phosphorylation and subsequent effector binding. We assume that all reactions describing ligand binding are parametrized by the same kinetic constants. The same holds true for all phosphorylation and all effector binding reactions. Structurally, the system can be dissected into three different layers.
Two binding layers describing ligand and effector binding, and one modification layer describing receptor phosphorylation. Ligand binding can be described by the reaction rule
D[0, X] + L ⇋ D[L, X], (10)
which can be expanded to the reaction rates d_{1}, d_{3} and d_{7}. Since all of these reactions are parametrized by the same kinetic parameters we combine them all to one single gross reaction
ROX + L ⇋ RLX (11)
The gross reaction
RX 0 ⇋ RXP (13)
describes binding site phosphorylation and corresponds to the reaction rule
D[X, O] ⇋ D[X, p] (14)
which can be expanded to the reaction rates d_{2} and d_{4}.
in which the lowercase letter concentration RXp occurs. However, the macroscopic variable we want to describe is RXP. Hence, one has to express RXp using macroscopic variables, which in this case is relatively simple. One can use the conservation relation
RXp = RXP  RXE. (16)
The gross reaction
RXp + E ⇋ RXE (17)
describes effector binding and corresponds to the rule
D[X, p] + E ⇋ D[X, E], (18)
Note that instead of RXP or RXE we could also choose RXp as a state variable of the reduced system. However, as degrees of phosphorylation (e.g. RXP) and levels of occupancy (e.g. RXE) are usually of interest and correspond to experimental readouts, we generally choose uppercase species as state variables. This has the additional advantage that then there is a conservation relation for each molecule that is not degraded or synthesized within each layer (see Equation 20).
b) Networks including all kind of interactions
To introduce the general reduced modeling concept, where all kinds of interactions are allowed, we start with an example to illustrate the main features.
Example: Again we consider the simplified insulin model introduced above. Now, we assume that ligand binding unidirectionally influences receptor phosphorylation which in turn is an essential precondition for effector binding. Ligand binding and effector binding do not interact directly. From this it follows that the reaction rates d_{1}, d_{3} and d_{7} are all parametrized by the same kinetic rate constants. The same holds true for the reaction rates d_{5} and d_{6}.
Observe, that these six equations are linearly dependent. RXp or RXE can be expressed by using a conservation relation
x = ROP + RLP = RXp + RXE = RXP. (23)
The connection between the two layers, i.e. the information exchange, is given by x = RXP (Equation 23) and xb = RXE. The sum of phosphorylated binding sites (x) is passed to the effector layer, the sum of occupied binding sites (xb) is passed to the receptor layer.
As already mentioned, the reaction rates that are merged together (d_{3} and d_{7} as well as d_{5} and d_{6}) have the same kinetic rate constants.
The approximation bases on the assumption that ligand and effector binding are completely independent.
Kinetic parameters and initial conditions for the small example system
Parameter  Literature value  Unit  Source 

k _{1}  0.001  nM ^{1} s ^{1}  [52] 
k _{1}  4·10^{4}  s ^{1}  [52] 
k _{2}  0  s ^{1}  ass. 
k _{2}  0.00385  s ^{1}  [51] 
k _{4}  0.0231  s ^{1}  [50] 
k _{4}  0.00385  s ^{1}  [51] 
k _{5}  0.033  nM ^{1} s ^{1}  [53] 
k _{5}  0.113  s ^{1}  [53] 
Summary: Only one variable has actually been eliminated by model reduction, but the example serves to illustrate the main elements of the method. These are:

Two kinds of real interactions between two processes exist: allornone interactions and graded interactions. The third possibility is that these two processes do not interact.

Modularization is achieved by analyzing interactions between processes. No graded interactions are allowed between modules which are called layers. Layers are only connected by allornone interactions.

Each layer is modeled independently of the others. Gross reactions are formulated that correspond to reactions of macrostates.

Dephosphorylation reactions of binding sites need special attention, as approximation of lowercase species is necessary. All other reactions can be formulated as in detailed mechanistic modeling, however, using also macroscopic species.

Concentrations of macroscopic species (or sums of them) are transfered between the layers. Between layers there is only signal flow, but no mass flow.

Combinatorial complexity is decreased, as allornone interactions between layers reduce the number of binding events and bound species by introducing a macroscopic description of the processes.

Reduced model equations can be obtained directly without previous generation of detailed mechanistic model equations.

Model equations can be understood and interpreted intuitively. Model variables are macroscopic quantities that are converted by rates following simple kinetics.
General considerations: In reaction networks including graded interactions the formulation of gross reaction rates as defined above becomes more difficult than only with allornone interactions. We again start by dissecting the whole network into layers. Processes coupled by graded interactions are merged into one layer. All binding and modification processes within a layer must be directly or indirectly linked by graded interactions, while the different layers are only coupled by allornone interactions. The coupled layers only exchange information about macroscopic variables, like phosphorylation degrees and levels of occupancy.
Now we assume that the processes of each layer form an isolated network and formulate a detailed kinetic model of solely these processes. Since processes from other layers are neglected in these considerations combinatorial variety is highly reduced. Each state of the submodel can be interpreted as a sum of states of the complete model, and each single reaction in the submodel represents a number of reactions in the complete network. Interestingly, all reactions of the complete network corresponding to a certain reaction in the subnetwork are parameterized by the same kinetic parameters. The reason for this can be found in the definition of layers. Since there are no graded interactions between layers, alterations in other layers do not change the kinetic properties inside the considered layer.
Observe, that we define the reactions of the isolated partial network as gross reactions. We stated earlier that if all reactions forming a gross reaction have the same kinetic parameters it can be formulated using the law of mass action. However, one has to accommodate these rates by including mass conservation relations to eliminate lowercase species from the description (as in Equation 9). As processes within a layer can be modeled separately each layer has to fulfill certain conservation relations. The receptor layer e.g. is characterized by a conservation relation for R. Due to these conservation relations within each layer the gross reactions for phosphorylations have to be formulated with capital letter concentrations like RX...P...X, which also comprise all phosphorylated species with additionally bound effectors. Gross rate kinetics are nevertheless characterized by lowercase concentrations. However, a mass conservation relation as Equation 9 does usually not allow to replace the lowercase concentrations by macroscopic variables. This results from the presence of graded interactions within the layers. If binding site phosphorylation is influenced by other processes of the same layer, e.g. ligand binding, more than one gross rate is needed to describe binding site phosphorylation. Therefore, the mass conservation relation as Equation 9 does only provide information about the sum of concentrations of species that may be dephosphorylated. However, we can formulate approximative equations
R_{ R }M_{1}...p...M_{ n }≈ c_{ I }·R_{ R }M_{1}....P....M_{ n }, 0 ≤ c_{ I }≤ 1 (29)
It is assumed that this fraction is identical for all species R_{ R }[M_{1}....P....M_{ n }]. Note that the factor c_{ I }is time dependent as the fraction of phosphorylated binding sites that is unoccupied may change over time. A detailed discussion of approximation quality will be given in the mathematical background section. The considered gross rates for the phosphorylation of R can be formulated using the law of mass action, if the resulting rate equations are accommodated using correction terms c_{ I }(Equations 29 and 30).
Consistent initial values
When comparing the reduced and detailed models the initial conditions have to be chosen that transformation equations hold for the starting point. This is guaranteed for the example system (Figure 5) if Equation 27 holds. If there is no detailed model, choice of consistent initial values also is of great importance. Without caring about this problem simulation can end up with negative concentrations, e.g. if initial concentrations satisfy RXE > ROP + RLP. Effector binding requires previous receptor phosphorylation and the receptor cannot be dephosphorylated while an effector is bound. Therefore, the number of phosphorylated receptors always has to be greater than or equal to the number of receptoreffector complexes. In general, for each connection between layers there is one inequality constraint, e.g. RXE ≤ ROP + RLP. Ignoring them results in physically infeasible systems. So there are two general possibilities to avoid infeasible systems. The first is to start with initial conditions for the detailed model and transform them into initial conditions for the reduced model (Equation 22). The other possibility is to ensure that the inequality constraints hold for each connection between layers.
Another aspect of choosing proper initial conditions arises from approximation of lowercase species (Equation 29). As there is division by R_{ R }X...P...X to calculate c_{ I }(Equation 30), all R_{ R }X...P...X, where P represents phosphorylation of a binding site, have to be larger than zero. However, if R_{ R }X...P...X equals zero, R_{ R }X...p...X also equals zero. This results in an 0/0 case, for which the limit is one. As this may be problematic during model simulation we recommend to choose initial conditions for all R_{ R }X...P...X larger than zero. Values very close to zero can be taken, if zero is the desired initial condition.
Mathematical background
In this section we analyze the presented reduced modeling method from a mathematical point of view. First, we introduce some general mathematical considerations about model reduction. Afterwards we will show that the layer based reduction method also fits into this general procedure. This will help us to evaluate the method and make statements about approximation errors.
General considerations
The relevant states
Now, we show that the layer based reduction method fits into the previously introduced general pattern of model reduction. Hence, we first have to define the two set of states ${\overrightarrow{z}}_{k}$ and ${\overrightarrow{z}}_{n}$ implying the transformation $\overrightarrow{\phi}(\overrightarrow{x})$. Second, we have to define algebraic equations that can be used to approximate the states ${\overrightarrow{z}}_{n}$. The definition of the states ${\overrightarrow{z}}_{k}$ is quite simple, since the layer based reduction method makes clear statements about the states of the reduced model (see step by step procedure). The state vector of the reduced model ${\overrightarrow{z}}_{k}$ comprises all states of the reduced model that are defined by an ODE.
which transforms the states ${\overrightarrow{z}}_{n}$ of Equation 32 to new states ${\overrightarrow{z}}_{n}$. Since I is the identity matrix, ${\overrightarrow{\tilde{z}}}_{k}={\overrightarrow{z}}_{k}$.
If one now replaces ${\overrightarrow{\tilde{z}}}_{n}$ in Equation 35 by what is given in Equation 36 and observes ${\overrightarrow{\tilde{z}}}_{k}={\overrightarrow{z}}_{k}$ the resulting reduced model again is ${\dot{\overrightarrow{z}}}_{k}={\overrightarrow{g}}_{1}({\overrightarrow{z}}_{k},\overrightarrow{\psi}({\overrightarrow{z}}_{k}))$. This means that the reduced model structure is completely determined by defining the states ${\overrightarrow{z}}_{k}$ plus the definition of dim(${\overrightarrow{z}}_{n}$) algebraic equations in order to approximate the neglected states. These equations are the ratio equations like Equation 27, transformed to $\overrightarrow{z}$coordinates. They are necessary to reconstruct the states of the detailed model. Alternatively, for building the reduced model directly it is not necessary to postulate these equations explicitly. They are already contained implicitly in the reduced model equations (confer Equation 25).
Approximation of neglected states
Borisov et al. showed that if this equation is fulfilled at a point of time t_{0} it will be fulfilled for all times t > t_{0}. This equation can be simplified by elementary transformations to
DOO·D 11  D 1O·DO 1 = 0, (38)
which is equivalent to the ratio Equation 27.
In both the examples here and in Figure 3 the four species occurring in the ratio equations form a reaction cycle. These cycles consist of two processes in different layers, which do not influence each other directly, e.g. ligand and effector binding. It is quite obvious that each reaction network that is decomposed into layers includes such independent cycles. All pairs of processes being located in distinct layers do not interact directly. Indirect interactions with at least one allornone interaction in between are possible. For each of these cycles one can formulate a ratio equation like that described above. Note that each ratio equation lowers the number of states of the reduced model by one.
The assumption of rapid equilibrium [47] for each pair of reactions describing the same process also leads to the same ratio equations if the equilibrium constant is eliminated from the equations. From this it follows that these ratio equations are at least approximately fulfilled if the rapid equilibrium assumption is a good approximation. Therefore, the approximation error vanishes completely if the system reaches thermodynamic equilibrium [47], as illustrated below. The ratio equations are used to approximate the states ${\overrightarrow{z}}_{n}$ as discussed above.
Approximation quality
g = DOO·DLp  DOp·DLO. (39)
with
a = k_{1} + k_{1} + k_{2} + k_{2}
and
u(t) = J_{1}·DLp  J_{2}·DOp  J_{3}·DLO + J_{4}·DOO
where the rates d_{1} and d_{3} are parametrized by k_{1} and k_{1} and the rates d_{2} and d_{4} are parametrized by k_{2} and k_{2}.
These equations show that both the steady state error as well as the maximal dynamic error decrease for increasing values of a, which corresponds to increasing values of the kinetic parameters k_{1}, k_{1}, k_{2} and k_{2}, and is zero for one of these values going to infinity. Note, that even for a large error g the reduced dynamic model equations may provide very good approximations for the states ${\overrightarrow{z}}_{k}$. This is the case if the erroneous states ${\overrightarrow{z}}_{n}$ that are approximated only have very little influence on ${\overrightarrow{z}}_{k}$, which is closely related to the system theoretical concept of observability [48]. However, since in most relevant cases even the states ${\overrightarrow{z}}_{n}$ are approximated quite well, this shall not be discussed in detail. Additional explanations and detailed calculations can be found in Additional files 1 and 2.
The modeling procedure – step by step
The following procedure provides a guide to build reduced models step by step. This procedure is used to generate a large model of insulin signaling that covers all processes discussed in the introduction (Additional file 3).
 1.
Identify all processes (Figure 4, white boxes) and their interactions (Figure 4, green and red lines).
 2.
Define layers: all processes that are coupled by graded interactions are within the same layer. Layers are coupled by allornone interactions or do not interact (Figure 4, blue boxes).
 3.
Model each layer individually.
 (a)
Define all sums of phosphorylated binding sites x_{ i }(e.g. x = ROP + RLP, see Equation 26b) and all sums of occupied binding sites x_{ i }b (e.g. xb = RXE). The terminal 'b' indicates that there is something bound to the binding site.
 (b)
Define the concentrations of all unoccupied phosphorylated binding sites that are needed as binding partners within the considered layer (e.g. RXp = x  xb). These species act as binding partners for binding reactions (see Figure 5).
 (c)
Define rules and reactions (including dephosphorylation of binding sites) as if there were no other layers and in particular as if there was no binding of effectors (see Figure 5). Use (uppercase) 'P' to indicate phosphorylation of binding sites and use another denotation (e.g.lowercase 'p') to indicate phosphorylation of regulatory sites.
 (d)
Translate each reaction into the corresponding rate by using the desired kinetic law. This step is analogous to detailed mechanistic modeling. For dephosphorylation reactions of binding sites multiply the expression describing dephosphorylation with (x_{ i } x_{ i }b)/x_{ i }using the appropriate x_{ i }. This ensures that only unoccupied binding sites are dephosphorylated (see Equation 26a).
 (e)
Optional: for each molecule that is not degraded or synthesized a conservation relation can be formulated (e.g. ROO = R_{0}  RLO  ROP  RLP, see Equation 26b).
 (f)
Construct ODEs as a sum of rates for each species that is used in this layer and not defined by an algebraic equation (see Equation 26c).
 4.
Additional information transfer between layers is allowed, as long as no additional graded interactions are introduced (e.g. R_{ activ }in Additional file 4).
This procedure also outlines how automation of the modeling procedure can be achieved. Steps 1) to 3c) most probably will be performed by the user, whereas expansion of rules and generation of rates and ODEs could be automated. The modeling procedure then remains in close similarity to automated rule based building of detailed mechanistic models by BioNetGen [19].
A larger example system
Optimization study of the larger example system
An optimization study was performed to analyze the worst case scenario within physiologic parameter ranges. Measures for reduction quality are the errors in XEF and DXXF. XEF corresponds to all species of effector F, that are bound to phosphorylated effector E, be E bound to the receptor or not. DXXF corresponds to all species of F that are bound to the receptor via E. In contrast to XEF, DXXF is not a state of the reduced model. It is the sum of all approximated states of the detailed model, where F is bound to the receptor via E. Depending on the scenario, both XEF and DXXF can correspond to activated effector F and therefore have physiologic relevance. Reconstruction of states (or sums of states as it is DXXF) of the detailed model is extensively discussed in Additional file 4. Starting from literature values, parameters were varied over four orders of magnitude to maximize the error in XEF and DXXF. This parameter interval should cover physiologic parameter ranges. Reduction quality is still very high and even the error of DXXF is within the range of measurement errors in typical experiments. Worst case parameters and simulation results are shown in Additional file 4.
Synthesis, degradation and transport of proteins
Up to now it was assumed that there is no protein synthesis or degradation. However, synthesis and degradation of free unmodified proteins is easy to handle. The rate of synthesis or degradation just has to be included in the differential equation for the free species. Degradation or synthesis of complexes is also possible though often not being as easy to realize. If a scaffold protein or even the receptor is to be degraded one has to observe that there exist lumped states in several layers that correspond to complexes with this effector or the receptor. In this case, the rate of degradation or synthesis has to be considered in different layers to guarantee consistency. Now more communication between the layers is necessary. The same rates that are modified by different correction terms to reflect complex composition occur in distinct layers. Transport between different compartments can be handled as degradation in one compartment and synthesis of the same complex in the other. Therefore, synthesis, degradation and transport of complexes is possible with the layer based formalism. Concise modular structure of the model with a minimum of information transfer between the layers is yielded when these processes are limited to free protein species.
Outline: application on insulin signaling
A model for the insulin signaling system which includes all events that were mentioned in the introduction can be built with 64 + 128 + 4 + 11 + 5 + 2 = 214 differential equations instead of 1.5·10^{8} in the detailed case. The 214 equations for the reduced model are composed as follows. 2^{6} = 64 equations in the receptor layer describe binding of two insulin molecules and phosphorylation of four binding sites (two for Shc and two for IRS). The 2^{7} = 128 equations of the IRS layer derive from six binding sites, each of them can be phosphorylated and unphosphorylated. IRS can be bound to the receptor and unbound. 4 equations are needed for the Shc layer (Shc binding to the receptor and becoming phosphorylated). SOS and Grb2 are merged into one layer. This allows SOS binding to Grb2 influencing Grb2 binding to IRS and Shc. The corresponding layer is described by 11 equations (six equation describing binding of complexes of Grb2 and SOS to IRS and Shc, five equations for free species, remember that SOS can be phosphorylated). The PI3K layer contains 5 differential equations (binding to four binding sites) and the SHP2 layer 2 (binding to IRS). Free species are considered. The two binding sites on the receptor for Shc and IRS in each case are assumed to be equivalent. Other underlying assumptions are specified in Table 1. Layer based reduced modeling of this system is demonstrated in Additional files 3 and 6.
Reduction ratios and combination with domainoriented lumping
As shown above for insulin signaling combinatorial complexity increases as the effector is subject of additional modification and binding events. The reduction potential of the layer based reduced modeling method strongly increases with increasing combinatorial complexity. For the small example system the number of equations is reduced by 20%, for the larger example system by 48%. The fraction of necessary equations for the signaling system as described in the introduction is reduced by 99.9999 %. This illustrates that even large systems can be described with a number of equations that can be handled by manual modeling. However, molecules with many binding sites or regulatory phosphorylation sites are still difficult to handle. Here combination with the domainoriented lumping technique [18] is possible. An example scenario is analyzed, underlying assumptions are specified in Table 1. There, case 1) is a very general setting, while case 2) is one of the most simple realistic special cases. In this case it is assumed that the phosphorylation state of binding sites does not influence phosphorylation of other binding sites on the same molecule. The layer based formalism requires still 214 equation, domainoriented lumping requires 212 equations, the combination of both methods yields a model with only 56 equations. To combine both approaches, first the layer based model is generated. Then the domainoriented approach is applied to each layer separately (see Additional files 7 and 8).
So, combination of layer based reduced modeling with the domainoriented approach [18] is possible. Under specific conditions, this results in an additional decrease of the fraction of necessary equations. The ideal scenario for combining both methods is the occurrence of molecules with many sites, where these sites do not or sparely interact. Then domainoriented lumping allows further strong reduction of the layer based reduced model. Reduction quality stays the same as the domainoriented approach provides exact lumping.
Conclusion
We present a reduced modeling approach which allows to tackle the problem of combinatorial complexity in signal transduction and regulation networks. For physiologic systems combinatorial complexity is dramatically decreased, as demonstrated on insulin signaling. Similarly to Pawson and Nash [23], Borisov et al. [22] and Conzelmann et al. [18] our method aims at a more macroscopic description of a system using levels of occupancy and phosphorylation.
A modularization principle is introduced. There, one has to distinguish between graded interactions and allornone interactions. Modules, which correspond to layers of signal transduction, are chosen such that they only interact via allornone interactions. All processes inside layers are connected directly or indirectly via graded interactions. For each layer we formulate a detailed reaction scheme comprising all processes of the current layer but neglecting all other processes. In the subsequent modeling step one has to formulate gross reaction rates for all of these reactions. A step by step procedure for building reduced modular models is given.
A potential drawback of the method is that even small changes to the assumptions of the model may lead to merging of layers which results in lower decrease of combinatorial complexity. This is the case if additional graded interactions between processes of different layers are introduced.
Mathematical analysis as well as simulation and optimization studies using insulin signaling models show that the approximation quality is excellent for relevant parameter settings. In the considered examples it is even possible to reconstruct the eliminated microstates with high accuracy. As we also showed the method allows an enormous reduction of the number of necessary ODEs compared with detailed combinatorial models.
The method can be combined with the domainoriented lumping technique [18], which can result in further reduction of the model.
Methods
Initial conditions
The total concentration of insulin receptor in hepatocytes was reported to be 10^{5} receptors per cell [49]. A hepatocyte is assumed to be a ball with diameter 20 μ m. This ball has a volume of 4.2·10^{12} l. With one mol standing for 6.0221415·10^{23} molecules there are 1.66·10^{19} mol receptor per cell. This is a concentration of 40 nM. The total concentration for IRS (E) was assumed to be the same as for Shc, which is 250 nM [9]. The total concentration for PI3K (F) was taken from literature [9]. It is assumed that the signalling system is totally inactive at the beginning. Initial conditions other than zero but close to zero (10^{20} nM) for the states which should start with an initial condition of zero were chosen. This prevents division by zero in correction terms of the rates of the reduced models and in the computations for consistent initial values. The states of the reduced model were chosen consistently such that transformation equations hold.
Parameter values
Insulin receptor autophosphorylation in vitro has a halflife of about 0.5 min (Figure 1 in [50]). Assuming first order kinetics, this corresponds to a rate constant of 0.0231 s^{1}. Insulin receptor dephosphorylation on the plasma membrane in vitro has a halflife of about 3 min (Figure 2 in [51]). Assuming first order kinetics, this corresponds to a rate constant of 0.00385 s^{1}. Parameters describing binding of IRS to the insulin receptor were originally reported to describe the binding of the p85 subunit to IRS [50].
Abbreviations
 nM:

nano molar (10^{9} mol·l^{1})
 ass.:

assumption
 ODE:

ordinary differential equation
Declarations
Acknowledgements
The authors acknowledge support from the german federal ministry of education and research (Bundesministerium für Bildung und Forschung, BMBF). This work was funded by the Network Systems Biology HepatoSys.
Authors’ Affiliations
References
 Sackmann A, Heiner M, Koch I: Application of Petri net based analysis techniques to signal transduction pathways. BMC Bioinformatics. 2006, 7: 482PubMed CentralView ArticlePubMed
 Klamt S, SaezRodriguez J, Lindquist JA, Simeoni L, Gilles ED: A methodology for the structural and functional analysis of signaling and regulatory networks. BMC Bioinformatics. 2006, 7: 56PubMed CentralView ArticlePubMed
 Gillespie DT: A rigorous derivation of the chemical master equation. Physica A. 1992, 188: 404425.View Article
 Asthagiri A, Lauffenburger D: A computational study of feedback effects on signal dynamics in a mitogenactivated protein kinase (MAPK) pathway model. Biotechnol Prog. 2001, 17: 227239.View ArticlePubMed
 Hatakeyama M, Kimura S, Naka T, Kawasaki T, Yumoto N, Ichikawa M, Kim JH, Saito K, Saeki M, Shirouzu M, Yokoyama S, Konagaya A: A computational model on the modulation of mitogenactivated protein kinase (MAPK) and Akt pathways in heregulininduced ErbB signalling. Biochem J. 2003, 373 (Pt 2): 451463.PubMed CentralView ArticlePubMed
 Haugh J, Schooler K, Well A, Wiley H, Lauffenburger D: Effect of epidermal growth factor receptor internalization on regulation of the phosphoslipase Cgamma1 signaling pathway. J Biol Chem. 1999, 274: 89588965.View ArticlePubMed
 Haugh J, Well A, Lauffenburger D: Mathematical modeling of epidermal growth factor receptor signaling through the phospholipase C pathway: mechanistic insights and predictions for molecular interventions. Biotechnol Bioeng. 2000, 70: 225238.View ArticlePubMed
 Kholodenko BN, Demin OV, Moehren G, Hoek JB: Quantification of Short Term Signaling by the Epidermal Growth Factor Receptor. J Biol Chem. 1999, 274 (42): 3016930181.View ArticlePubMed
 Moehren G, Markevich N, Demin O, Kiyatkin A, Goryanin I, Hoek JB, Kholodenko BN: Temperature dependence of the epidermal growth factor receptor signaling network can be accounted for by a kinetic model. Biochemistry. 2002, 41: 30620.View ArticlePubMed
 Schoeberl B, EichlerJonsson C, Gilles ED, Müller G: Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors. Nat Biotechnol. 2002, 20 (4): 370375.View ArticlePubMed
 Sedaghat AR, Sherman A, Quon MJ: A mathematical model of metabolic insulin signaling pathways. Am J Physiol Endocrinol Metab. 2002, 283 (5): E1084101.View ArticlePubMed
 Faeder JR, Hlavacek WS, Reischl I, Blinov ML, Metzger H, Redondo A, Wofsy C, Goldstein B: Investigation of early events in Fc(epsilon) RImediated signaling using a detailed mathematical model. J Immunol. 2003, 170 (7): 376981.View ArticlePubMed
 Blinov ML, Faeder JR, Goldstein B, Hlavacek WS: A network model of early events in epidermal growth factor receptor signaling that accounts for combinatorial complexity. Biosystems. 2006, 83 (2–3): 136151.View ArticlePubMed
 Hlavacek WS, Faeder JR, Blinov ML, Perelson AS, Goldstein B: The complexity of complexes in signal transduction. Biotechnol Bioeng. 2004, 84 (7): 783794.View Article
 MortonFirth CJ, Bray D: Predicting temporal fluctuations in an intracellular signalling pathway. J Theor Biol. 1998, 192: 117128.View ArticlePubMed
 Shimizu TS, Novere NL, Levin MD, Beavil AJ, Sutton BJ, Bray D: Molecular model of a lattice of signalling proteins involved in bacterial chemotaxis. Nat Cell Biol. 2000, 2 (11): 792796.View ArticlePubMed
 Faeder JR, Blinov ML, Goldstein B, Hlavacek WS: Combinatorial complexity and dynamical restriction of network flows in signal transduction. IEE Systems Biology. 2005, 2: 515.View ArticlePubMed
 Conzelmann H, SaezRodriguez J, Sauter T, Kholodenko B, Gilles E: A domainoriented approach to the reduction of combinatorial complexity in signal transduction networks. BMC Bioinformatics. 2006, 7: 34PubMed CentralView ArticlePubMed
 Blinov ML, Faeder JR, Goldstein B, Hlavacek WS: BioNetGen: software for rulebased modeling of signal transduction based on the interactions of molecular domains. Bioinformatics. 2004, 20 (17): 328991.View ArticlePubMed
 Faeder JR, Blinov ML, Hlavacek WS: Graphical rulebased representation of signal transduction networks. Proc ACM Symp Appl Computing. 2005, 133140.
 Blinov ML, Yang J, Faeder JR, Hlavacek WS: Graph theory for rulebased modeling of biochemical networks. Proceedings of BioCONCUR 2005, A workshop on concurrent models in molecular biology. 2005
 Borisov NM, Markevich NI, Hoek JB, Kholodenko BN: Signaling through receptors and scaffolds: independent interactions reduce combinatorial complexity. Biophys J. 2005, 89 (2): 951966.PubMed CentralView ArticlePubMed
 Pawson T, Nash P: Assembly of cell regulatory systems through protein interaction domains. Science. 2003, 300 (5618): 445452.View ArticlePubMed
 Khan A, Pessin J: Insulin regulation of glucose uptake: a complex interplay of intracellular signalling pathways. Diabetologia. 2002, 45 (11): 147583.View ArticlePubMed
 Saltiel AR, Pessin JE: Insulin signaling pathways in time and space. Trends Cell Biol. 2002, 12 (2): 6571.View ArticlePubMed
 Saltiel A, Kahn C: Insulin signalling and the regulation of glucose and lipid metabolism. Nature. 2001, 414 (6865): 799806.View ArticlePubMed
 Taniguchi CM, Emanuelli B, Kahn CR: Critical nodes in signalling pathways: insights into insulin action. Nat Rev Mol Cell Biol. 2006, 7 (2): 8596.View ArticlePubMed
 Chang L, Chiang SH, Saltiel AR: Insulin signaling and the regulation of glucose transport. Mol Med. 2004, 10 (7–12): 6571.PubMed CentralPubMed
 Plum L, Belgardt BF, Brüning JC: Central insulin action in energy and glucose homeostasis. J Clin Invest. 2006, 116 (7): 17616.PubMed CentralView ArticlePubMed
 Mounier C, Posner BI: Transcriptional regulation by insulin: from the receptor to the gene. Can J Physiol Pharmacol. 2006, 84 (7): 71324.View ArticlePubMed
 Leng Y, Karlsson HKR, Zierath JR: Insulin signaling defects in type 2 diabetes. Rev Endocr Metab Disord. 2004, 5 (2): 1117.View ArticlePubMed
 Chakraborty C: Biochemical and molecular basis of insulin resistance. Curr Protein Pept Sci. 2006, 7 (2): 11321.View ArticlePubMed
 Ludvigsson J: Why diabetes incidence increases – a unifying theory. Ann N Y Acad Sci. 2006, 1079: 37482.View ArticlePubMed
 Musi N, Goodyear LJ: Insulin resistance and improvements in signal transduction. Endocrine. 2006, 29: 7380.View ArticlePubMed
 Salsali A, Nathan M: A review of types 1 and 2 diabetes mellitus and their treatment with insulin. Am J Ther. 2006, 13 (4): 34961.View ArticlePubMed
 Stumvoll M, Goldstein BJ, van Haeften TW: Type 2 diabetes: principles of pathogenesis and therapy. Lancet. 2005, 365 (9467): 133346.View ArticlePubMed
 Hirsch IB: Insulin analogues. N Engl J Med. 2005, 352 (2): 17483.View ArticlePubMed
 Luo R, Beniac D, Fernandes A, Yip C, Ottensmeyer F: Quaternary structure of the insulininsulin receptor complex. Science. 1999, 285 (5430): 107780.View ArticlePubMed
 White M: The IRSsignalling system: a network of docking proteins that mediate insulin action. Mol Cell Biochem. 1998, 182 (1–2): 311.View ArticlePubMed
 Gual P, MarchandBrustel YL, Tanti JF: Positive and negative regulation of insulin signaling through IRS1 phosphorylation. Biochimie. 2005, 87: 99109.View ArticlePubMed
 Pirola L, Johnston A, Obberghen EV: Modulation of insulin action. Diabetologia. 2004, 47 (2): 17084.View ArticlePubMed
 Hartwell L, Hopfield J, Leibler S, Murray A: From molecular to modular cell biology. Nature. 1999, 402 (6761 Suppl): C4752.View ArticlePubMed
 SaezRodriguez J, Kremling A, Conzelmann H, Bettenbrock K, Gilles ED: Modular Analysis of Signal Transduction Networks. IEEE Contr Syst Mag. 2004, 24 (4): 3552.View Article
 SaezRodriguez J, Kremling A, Gilles ED: Dissecting the puzzle of life: Modularization of signal transduction networks. Comput Chem Eng. 2005, 29 (3): 619629.View Article
 Ederer M, Sauter T, Bullinger E, Gilles ED, Allgöwer F: An Approach for Dividing Models of Biological Reaction Networks into Functional Units. Simulation. 2003, 79 (12): 703716.View Article
 Borisov NM, Markevich NI, Hoek JB, Kholodenko BN: Trading the microworld of combinatorial complexity for the macroworld of protein interaction domains. Biosystems. 2006, 83 (2–3): 15266.PubMed CentralView ArticlePubMed
 Heinrich R, Schuster S: The Regulation of Cellular Systems. 1996, Chapman & HallView Article
 Isidori A: Nonlinear control systems. 1995, Springer, 77105. 3, chap. Global decomposition of control systems
 Rother K, Imai Y, Caruso M, Beguinot F, Formisano P, Accili D: Evidence that IRS2 phosphorylation is required for insulin action in hepatocytes. J Biol Chem. 1998, 273 (28): 174917.View ArticlePubMed
 Faure R, Baquiran G, Bergeron J, Posner B: The dephosphorylation of insulin and epidermal growth factor receptors. Role of endosomeassociated phosphotyrosine phosphatase(s). J Biol Chem. 1992, 267 (16): 1121521.PubMed
 Drake P, Bevan A, Burgess J, Bergeron J, Posner B: A role for tyrosine phosphorylation in both activation and inhibition of the insulin receptor tyrosine kinase in vivo. Endocrinology. 1996, 137 (11): 49608.PubMed
 Wanant S, Quon M: Insulin receptor binding kinetics: modeling and simulation studies. J Theor Biol. 2000, 205 (3): 35564.View ArticlePubMed
 Felder S, Zhou M, Hu P, Ureña J, Ullrich A, Chaudhuri M, White M, Shoelson S, Schlessinger J: SH2 domains exhibit highaffinity binding to tyrosinephosphorylated peptides yet also exhibit rapid dissociation and exchange. Mol Cell Biol. 1993, 13 (3): 144955.PubMed CentralView ArticlePubMed
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.