Exact model reduction of combinatorial reaction networks
- Holger Conzelmann^{1}Email author,
- Dirk Fey^{2} and
- Ernst D Gilles^{1}
DOI: 10.1186/1752-0509-2-78
© Conzelmann et al; licensee BioMed Central Ltd. 2008
Received: 11 March 2008
Accepted: 28 August 2008
Published: 28 August 2008
Abstract
Background
Receptors and scaffold proteins usually possess a high number of distinct binding domains inducing the formation of large multiprotein signaling complexes. Due to combinatorial reasons the number of distinguishable species grows exponentially with the number of binding domains and can easily reach several millions. Even by including only a limited number of components and binding domains the resulting models are very large and hardly manageable. A novel model reduction technique allows the significant reduction and modularization of these models.
Results
We introduce methods that extend and complete the already introduced approach. For instance, we provide techniques to handle the formation of multi-scaffold complexes as well as receptor dimerization. Furthermore, we discuss a new modeling approach that allows the direct generation of exactly reduced model structures. The developed methods are used to reduce a model of EGF and insulin receptor crosstalk comprising 5,182 ordinary differential equations (ODEs) to a model with 87 ODEs.
Conclusion
The methods, presented in this contribution, significantly enhance the available methods to exactly reduce models of combinatorial reaction networks.
Background
A large problem in mathematical modeling of biochemical signal transduction networks is combinatorial complexity [1]. Due to the occurrence of receptors and scaffold proteins with high numbers of binding domains and binding partners the number of feasible molecular species is enormous [1–5]. Many existing models evade the problem of combinatorial variety by assuming reduced, heuristic network structures that focus on a restricted number of molecular species and reactions [6–14]. Faeder et al. [15] showed by simulation studies that in combinatorial reaction networks only a relatively small part of the network might be active which means that the concentration of many species is negligible low. By eliminating these species as well as the associated reactions a fairly reduced model can be obtained. Faeder et al. [15] also showed that the predictions of such a reduced model match those of the complete one for the original parameter values quite well, but that the reduced model is not predictive over a larger range of parameter values. Even very small perturbations in the parameters may result in large approximation errors. Since the real kinetic parameters are usually unknown, these findings of Faeder et al. [15] indicate that a heuristically derived model structure will mostly be insufficient to approximate the dynamics of the real signaling network. This result is confirmed by Conzelmann et al. [16] who showed, by discussing a scaffold protein with only three binding domains, that reasonable, heuristically derived, model reductions may lead to significant approximation errors even for very small reaction networks.
An alternative approach to tackle combinatorial complexity is stochastic simulation. One possibility is an agent-based approach in which each protein is considered to be an autonomous individual in the reaction network. This approach restricts the number of elements that have to be considered to the total number of protein copies, while the number of feasible multiprotein complexes, which can easily grow to billions [1, 3], often exceeds this number by far. However, the computational cost for that kind of stochastic simulations can also be extremely high. Another possibility how stochastic models can help to reduce combinatorial complexity has been presented by Lok et al. [17]. The enormous complexity of the stochastic model is reduced by a new approach which incorporates complexes and reactions only when they are needed as the simulation proceeds [17, 18]. However, it is much harder to analyze the dynamic behavior of stochastic models or to identify the model parameters from measurements [18].
Another, ODE based, approach has been introduced by Blinov et al. [19]. The modeling tool BIO NET GEN allows a rule-based model specification, which is automatically expanded to a complete mechanistic ODE model. BIO NET GEN has been used to create a number of signaling models including EGF receptor signaling and Fcϵ RI signaling [20, 21].
Borisov et al. [22, 23] proposed an alternative approach which adopts the point of view that the fundamental elements of signal transduction are domains instead of molecular species [24]. Borisov et al. suggested to substitute the common mechanistic network description that includes all individual molecular species by a macrodescription. In this context, the term micro-state is used to describe individual multiprotein complexes, whereas the term macro-state refers to large sets of micro-states sharing a certain characteristic like phosphorylation of a defined binding domain. These macro-states correspond to descriptive biological quantities like phosphorylation degrees or levels of occupancy. The goal of this approach is the generation of a dynamic ODE model describing the concentrations of these macro-states. Borisov et al. [22] show for instance that in the case of independent binding domains a reduced number of ODEs is sufficient to describe the dynamics of the macro-states accurately. This approach has been extended by Conzelmann et al. [16] who introduced a systematic procedural method. The starting point of this method is a complete mechanistic ODE model which is subsequently reduced [16]. The reduction method bases on a hierarchically structured state space transformation. However, Conzelmann et al. [16] only provide a general transformation pattern for scaffold proteins with numerous single protein ligands. If these ligands can recruit further signaling proteins or the considered scaffold forms a dimer, the transformation pattern has to be extended. In this contribution, we present these extensions that enhance the applicability of the method to all kinds of signaling networks. Furthermore, we introduce a novel approach to directly generate the reduced equations. This approach circumvents the computationally expensive set-up of a full combinatorial model and its subsequent reduction.
The Section Mathematical Prerequisites will briefly discuss the mathematical concepts that are used. The main results are discussed in three parts: Section Exact Model Reduction presents the extension of the transformation pattern introduced by Conzelmann et al. [16]. Section Reduced Order Modeling of Combinatorial Reaction Networks introduces the novel approach for directly generating reduced models. Section Example: EGF and Insulin Receptor Crosstalk discusses the methods practical applicability and its benefits by generating a fairly reduced model of EGF and insulin receptor crosstalk.
Methods
Mathematical prerequisites
where $\overrightarrow{x}$(t) denotes the n-dimensional vector of all dynamic states or variables and $\overrightarrow{u}$(t) represents the m-dimensional vector of all external input signals. The vector $\overrightarrow{y}$(t) comprises all q output variables of the system, which either correspond to measured quantities or more general to all quantities of interest. The vector field $\overrightarrow{f}$ and the vector valued function $\overrightarrow{h}$ do have appropriate dimensions.
Since we assume that the output $\overrightarrow{y}$ does include all essential quantities, we are solely interested in the systems input/output behavior. It is very essential to stress that exact reducibility is solely determined by the system's structure and the definition of input and output variables. Thus, the only way to influence exact reducibility is to modify a system or to change the input and output variables. However, in most cases neither the model structure nor the choice of input and output variables is discretionary. Exact reduction methods should not try to influence exact reducibility of a model but provide the information whether a given system is exactly reducible and how it can be reduced. The question of whether a model is exactly reducible is closely related to the control theoretical concepts of observability and controllability which will be introduced below. In control theory models that are not exactly reducible are called minimal realizations [25]. If a model is no minimal realization it comprises uncontrollable or unobservable states, which can be eliminated without changing the systems input/output behavior. The elimination of unobservable and uncontrollable states can be achieved by a formal dissection of the model's state space into observable and controllable, observable but not controllable, controllable but not observable as well as neither observable nor controllable subsystems. Such a dissection is called Kalman decomposition [26].
Control theory provides a couple of methods that facilitate the separation of controllable and uncontrollable as well as of observable and unobservable states. However, a deficiency of these methods is that they are either developed for linear systems or they are only designed for small nonlinear models [25], but they are not applicable to very large nonlinear models of biological signaling cascades.
The methods we developed and present facilitate a Kalman decomposition of combinatorial reaction networks. Interestingly, for the considered type of systems a Kalman decomposition can be achieved by a linear state space transformation that is additionally independent of the systems parametrization or the choice of input and output variables. Naturally, a change of the parameter values or the input and output variables can affect the number of observable and controllable states but the proposed transformation always provides a Kalman decomposition. This statement can be verified using the approach discussed in the Section Generality of the Method.
Mathematical characterization of processes and process interactions
According to Conzelmann et al. [28] the following process interaction types can be distinguished
• non-interacting processes
Complete independence implies that the kinetic association and dissociation constants of one domain do not change upon ligand binding on the other domain. Hence, it follows for the parameters k_{2} = k_{1}, k_{-2} = k_{-1}, k_{4} = k_{3} and k_{-4} = k_{-3}.
• unidirectionally interacting processes
The binding of one ligand, e.g. ligand L, is not influenced by binding of the other one. However, L binding does change the kinetic properties of the other domain. In this case, only the conditions k_{2} = k_{1} and k_{-2} = k_{-1} have to be fulfilled.
• mutually interacting processes
This is the most general case. Binding of a ligand has an influence on binding of the other ligand and vice versa. In this case all parameters can have different values.
In addition to these, Koschorreck et al. [3] introduce
• all-or-none interactions
All-or-none interactions are a special case of mutual interactions. An mutual interaction between two processes is called all-or-none interaction, if the reaction cycle given by the four reactions of our example degenerates to a reaction chain. In real biochemical networks such interactions usually are given by domain phosphorylation and subsequent ligand binding. Mostly phosphorylation can be considered as necessary precondition for ligand binding. On the other side ligand binding prevents dephosphorylation of the domain for steric reasons. To realize an all-or-none interaction in our example, one has to choose k_{2} = k_{-2} = k_{3} = k_{-3} = 0. In this case the species R(0, E) will not occur and the remaining reactions r_{1} and r_{4} form a reaction chain.
As can be seen from the examples provided below as well as from examples discussed in other publications [22, 23, 16, 3], the absence of interactions or the occurrence of unidirectional and all-or-none interactions facilitate model reduction and modularization.
State space transformations
In the case of nonlinear transformations, the transformed model equations can be deduced following the same procedure. However, note that the inversion of a nonlinear transformation can be extremely difficult.
Observability
in which A, B and C are constant matrices of appropriate dimensions and the initial conditions are given by ${[{\overrightarrow{x}}_{1}(0),{\overrightarrow{x}}_{2}(0)]}^{T}={[{\overrightarrow{x}}_{1,0},{\overrightarrow{x}}_{2,0}]}^{T}$. Obviously, the variables denoted as ${\overrightarrow{x}}_{2}$ do not have any influence on the output variables $\overrightarrow{y}$. Hence, any initial states whose values for ${\overrightarrow{x}}_{1,0}$ coincide result in identical outputs for arbitrary initial conditions ${\overrightarrow{x}}_{2,0}$. The differences in the states ${\overrightarrow{x}}_{2}$ cannot be observed considering these outputs.
The number of observable states is determined by the dimension of the so-called observability space. For linear systems it can be calculated as
d = rank(Q) with Q = [C, C A, . . . C A^{n-1}]^{ T }
The first d linearly independent rows of Q build a basis for the observability space $\mathcal{O}$[25]. For d = n the system is called observable.
These considerations imply that an unobservable system always can be reduced without affecting the dynamics of the output variables. In Equation 6 a reduced system would exclusively comprise the ODEs for the state variables ${\overrightarrow{x}}_{1}$. Such reductions are exact with respect to the output.
Controllability
In this case the state variables ${\overrightarrow{x}}_{2}$ cannot be influenced by the inputs. Hence, the chosen input does not allow to control the system in any desired way, why the system is said to be uncontrollable [25]. In analogy to the considerations about observability, there also exists a controllability space $\mathcal{C}$ whose dimension as well as its basis can be deduced from the matrix
P = [B A B . . . A^{n-1}B]
Uncontrollable systems also can be reduced without affecting the dynamics of the output, if an additional assumption is fulfilled. If the dynamics of the uncontrollable subsystem (${\overrightarrow{x}}_{2}$ in Example 8) already decayed, ${\overrightarrow{x}}_{2}$ can be replaced by its steady state value which in the regarded example corresponds to zero.
Results and discussion
Exact model reduction
Short review
Conzelmann et al. have introduced linear transformations that realize a Kalman decomposition for models of scaffold proteins and receptors with single protein ligands [16]. The term single protein ligand indicates that one only considers the multi domain protein and its direct binding partners but no additional binding or modification processes at these ligands.
The model reduction procedure suggested by Conzelmann et al. [16] can be divided into three essential steps.
Step 1: One generates a complete mechanistic ODE model of the considered combinatorial reaction network using e.g. BIO NET GEN or ALC [19, 29]. Furthermore, one has to define input and output variables.
In analogy to Equation 6 the states ${\overrightarrow{z}}_{2}$ are unobservable. Choosing an invertible transformation matrix assures that the system's dynamics are preserved, and the original states can be retrieved from the new ones at any time as long as none of the transformed equations are eliminated.
Scaffolds with multiprotein ligands
Many scaffold proteins or receptors often recruit other scaffolds which in turn can be phosphorylated and/or bind further ligands. A prominent example is the scaffold IRS which binds to insulin receptor and can recruit numerous other ligands like Grb2 or PI3K [30]. We will call scaffolds like IRS multiprotein ligands. Note that, in general, these multiprotein ligands can either bind single protein ligands or again other multiprotein ligands. Heterodimerization, as it occurs in the ErbB signaling network [31], also fits into this pattern. Homodimerization, on the other hand, will be excluded from the following considerations. Due to the symmetry of homodimeric complexes homodimerization has to be handled differently as will be discussed later. Scaffolds with single proteins ligands as discussed by Conzelmann et al. [16] can be considered as a special case of what we examine here. The main difference between these multiprotein and single protein ligand systems is the formation of long protein chains. For this reason, we will focus on this phenomenon and its mathematical treatment. Later we will exemplify how branched multiprotein ligand systems can be handled considering a model of EGF and insulin receptor crosstalk.
Let us focus on a receptor protein R which provides n binding domains. We take the assumption that each domain i can bind an effector protein ${E}_{1}^{i}$ which in turn can recruit another effector protein ${E}_{2}^{i}$ till finally ${E}_{m-1}^{i}$ binds ${E}_{m}^{i}$. In order to reduce the number of indices we also presume that each chain of effector proteins consists of m proteins. Finally, we only regard binding processes and neglect all domain phosphorylations. Thus, each receptor domain can be either unoccupied or occupied by a multiprotein ligand consisting of one to m proteins which results in (m + 1)^{ n }distinct receptor complexes. Furthermore, the m effectors that form the different multiprotein ligands for one single receptor domain can build $\frac{m(m+1)}{2}$ distinct complexes. According to these examinations the total number of feasible multiprotein species is ${(1+m)}^{n}+\frac{m(m+1)}{2}n$.
General Transformation Pattern
In analogy to the considerations made by Conzelmann et al. [16] we require a state space transformation which facilitates a Kalman decomposition of the reaction network. The transformation pattern introduced by Conzelmann et al. [16] can be divided into tiers that describe levels of occupency of different order. The concept of using levels of occupancy as new variables is problematic for the considered multiprotein ligand systems. The term level of occupancy implies a certain hierarchy among the signaling proteins, which is certainly given in the single protein ligand scenario where one scaffold can bind numerous other effector proteins. It is obvious that in such reaction networks the scaffold takes up a prominent position which suggests the consideration of its occupancy levels. In a system that involves numerous scaffolds a clear hierarchy is missing, and the question arises which occupancy levels should be considered. In most cases, an intuitive hierarchy will be automatically chosen. For example, in the case of insulin signaling, it is quite suggestive to choose the insulin receptor as central protein of the cascade. Due to representational reasons we also assume a hierarchy in our examples with R being the central protein. However, if one e.g. considers heterodimerization of two ErbB receptors [31] it is not apparent which receptor takes up a more prominent position than the other one. Another problem is that the definition of occupancy levels for multiprotein ligand systems is not as unique as for single protein ligand systems. The quantity [R(${E}_{1}^{1}$, *, ..., *)], which can be interpreted as an occupancy level, describes all receptor species whose first domain is occupied by the single protein ${E}_{1}^{1}$ excluding all species in which ${E}_{1}^{1}$ has bound any further effectors or scaffolds. [R(${E}_{1}^{1}$(*), *, ..., *)] on the other site represents an alternative type of occupancy level which does not exclude the previously mentioned multiprotein complexes.
A more general transformation pattern is required that avoids the implication of molecular hierarchy, and is consistent with the transformation pattern introduced by Conzelmann et al. [16]. These requirements are met by the introduction of so-called occurrence levels. Occurrence levels always refer to a certain molecular subcomplex and correspond to the sums of all multiprotein species that comprise this subcomplex. Thus, for each individual molecular species one can define a respective occurrence level. If these occurrence levels are used to replace the original model states, this defines a linear and invertible transformation. The proposed transformation pattern also preserves the hierarchical structure of the transformation matrix. The 0^{th} tier of the transformation matrix, as introduced by Conzelmann et al. [16], includes the overall concentrations of all involved signaling proteins. It corresponds to the occurrence levels of individual proteins which are a very special type of subcomplex. The 1^{st} tier includes the occurrence levels of all possible two protein subcomplexes. In the case of scaffolds with single protein ligands this directly corresponds to the levels of occupancy. According to this pattern the following tiers of the transformation matrix respectively comprise all subcomplexes consisting of three, four and more proteins. If phosphorylations occur in the considered reaction network the phosphate groups have to be treated as additional molecules. For example a scaffold with one phosphorylated domain is regarded as a two molecule complex.
For the simplified case introduced above the new transformed states can be specified as follows. The 0^{th} tier comprises the states [R(*, ..., *)] and [${E}_{j}^{i}$(*)], while the first one includes the states [R(*, ..., *, ${E}_{1}^{i}$(*), *, ..., *)] as well as $[{E}_{j}^{i}({E}_{j+1}^{i}(\ast ))]$. The occurrence levels that refer to all three molecule complexes are [R(*, ..., *, ${E}_{1}^{i}$(*), *, ..., *, ${E}_{1}^{k}$(*), *, ..., *)], [R(*, ..., *, ${E}_{2}^{i}$(*), *, ..., *)] and $[{E}_{j}^{i}({E}_{j+2}^{i}(\ast ))]$. The subsequent tiers are defined according to this pattern, and the last tier only comprises the single micro-state $[R({E}_{m}^{1},\mathrm{...},{E}_{m}^{n})]$. The fact that each individual molecular species can be uniquely linked to an associated occurrence level suggests that the transformation is invertible. This can also be proven by using the mathematical induction as described for single protein ligands in Conzelmann et al. [16] (data not shown).
Examples
Example 1: Six signaling proteins
Reaction rules describing the Example depicted in Figure 2A.
R(0, *) | + | L | ⇋ | R(L, *) | k_{1}, k_{-1} |
R(0, 0) | + | E_{1}(*) | ⇋ | R(0, E_{1}(*)) | k_{2}, k_{-2} |
R(L, 0) | + | E_{1}(*) | ⇋ | R(L, E_{1}(*)) | k_{3}, k_{-3} |
E_{1}(0, 0) | + | E_{2}(*) | ⇋ | E_{1}(E_{2}(*)) | k_{4}, k_{-4} |
R(*, E_{1}) | + | E_{2}(*) | ⇋ | R(*, E_{1}(E_{2}(*))) | k_{5}, k_{-5} |
E_{2}(0, 0) | + | E_{3}(*) | ⇋ | E_{2}(E_{3}(*)) | k_{6}, k_{-6} |
E_{1}(*, E_{2}) | + | E_{3}(*) | ⇋ | E_{1}(*, E_{3}(*)) | k_{7}, k_{-7} |
E_{3}(0, 0) | + | E _{4} | ⇋ | E_{3}(E_{4}) | k_{8}, k_{-8} |
E_{2}(*, E_{3}) | + | E _{4} | ⇋ | E_{2}(*, E_{4}) | k_{9}, k_{-9} |
Hierarchical transformation that realizes a Kalman decomposition for the example system depicted in Figure 2A.
[R(*, *)] | = | [R(0, 0)] + [R(0, E_{1})] + [R(0, E_{2})] + [R(0, E_{3})] + [R(0, E_{4})] + [R(L, 0)] + [R(L, E_{1})] + [R(L, E_{2})] + [R(L, E_{3})] + [R(L, E_{4})] |
[E_{1}(*)] | = | [E_{1}(0)] + [E_{1}(E_{2})] + [E_{1}(E_{3})] + [E_{1}(E_{4})] + [R(0, E_{1})] + [R(0, E_{2})] + [R(0, E_{3})] + [R(0, E_{4})] + [R(L, E_{1})] + [R(L, E_{2})] + [R(L, E_{3})] + [R(L, E_{4})] |
[E_{2}(*)] | = | [E_{1}(E_{2})] + [E_{1}(E_{3})] + [E_{1}(E_{4})] + [E_{2}(0)] + [E_{2}(E_{3})] + [E_{2}(E_{4})] + [R(0, E_{2})] + [R(0, E_{3})] + [R(0, E_{4})] + [R(L, E_{2})] + [R(L, E_{3})] + [R(L, E_{4})] |
[E_{3}(*)] | = | [E_{1}(E_{3})] + [E_{1}(E_{4})] + [E_{2}(E_{3})] + [E_{2}(E_{4})] + [E_{3}(0)] + [E_{3}(E_{4})] + [R(0, E_{3})] + [R(0, E_{4})] + [R(L, E_{3})] + [R(L, E_{4})] |
[E_{4}(*)] | = | [E_{1}(E_{4})] + [E_{2}(E_{4})] + [E_{3}(E_{4})] + [E_{4}(0)] + [R(0, E_{4})] + [R(L, E_{4})] |
[R(L, *)] | = | [R(L, 0)] + [R(L, E_{1})] + [R(L, E_{2})] + [R(L, E_{3})] + [R(L, E_{4})] |
[R(*, E_{1}(*))] | = | [R(0, E_{1})] + [R(0, E_{2})] + [R(0, E_{3})] + [R(0, E_{4})] + [R(L, E_{1})] + [R(L, E_{2})] + [R(L, E_{3})] + [R(L, E_{4})] |
[E_{1}(E_{2}(*)] | = | [E_{1}(E_{2})] + [E_{1}(E_{3})] + [E_{1}(E_{4})] + [R(0, E_{2})] + [R(0, E_{3})] + [R(0, E_{4})] + [R(L, E_{2})] + [R(L, E_{3})] + [R(L, E_{4})] |
[E_{2}(E_{3}(*))] | = | [E_{1}(E_{3})] + [E_{1}(E_{4})] + [E_{2}(E_{3})] + [E_{2}(E_{4})] + [R(0, E_{3})] + [R(0, E_{4})] + [R(L, E_{3})] + [R(L, E_{4})] |
[E_{3}(E_{4}(*))] | = | [E_{1}(E_{4})] + [E_{2}(E_{4})] + [E_{3}(E_{4})] + [R(0, E_{4})] + [R(L, E_{4})] |
[R(L, E_{1}(*))] | = | [R(L, E_{1})] + [R(L, E_{2})] + [R(L, E_{3})] + [R(L, E_{4})] |
[R(*, E_{2}(*))] | = | [R(0, E_{2})] + [R(0, E_{3})] + [R(0, E_{4})] + [R(L, E_{2})] + [R(L, E_{3})] + [R(L, E_{4})] |
[E_{1}(E_{3}(*)] | = | [E_{1}(E_{3})] + [E_{1}(E_{4})] + [R(0, E_{3})] + [R(0, E_{4})] + [R(L, E_{3})] + [R(L, E_{4})] |
[E_{2}(E_{4}(*))] | = | [E_{1}(E_{4})] + [E_{2}(E_{4})] + [R(0, E_{4})] + [R(L, E_{4})] |
[R(L, E_{2}(*))] | = | [R(L, E_{2})] + [R(L, E_{3})] + [R(L, E_{4})] |
[R(*, E_{3}(*))] | = | [R(0, E_{3})] + [R(0, E_{4})] + [R(L, E_{3})] + [R(L, E_{4})] |
[E_{1}(E_{4}(*)] | = | [E_{1}(E_{4})] + [R(0, E_{4})] + [R(L, E_{4})] |
[R(L, E_{3}(*))] | = | [R(L, E_{3})] + [R(L, E_{4})] |
[R(*, E_{4}(*))] | = | [R(0, E_{4})] + [R(L, E_{4})] |
[R(L, E_{4}(*))] | = | [R(L, E_{4})] |
Due to the absence of protein production and degradation the six states of the 0^{th} tier remain all constant. Thus, these six ODEs can be eliminated in the considered example. First, we will discuss the case that E_{4} binding is unidirectionally influenced by E_{3} binding. In this case, our transformation does not allow any exact reduction of the model, neither for the output variables ${\overrightarrow{y}}_{1}$ nor for ${\overrightarrow{y}}_{2}$. Interestingly, the transformed model equations can be dissected into five modules, which are all unidirectionally coupled. This model structure directly resembles the interaction pattern between the five considered binding processes. In fact each of the modules describes one of these five processes. However, the modules differ in size and structure.
The first module which describes the recruitment of L to the receptor only consists of one differential equation. The second, third, fourth and fifth module comprises two, three, four and five states, respectively. Another nice property of the transformed system is the concurrently achieved modularization of the kinetic parameters. For instance, the L binding module only contains the parameters k_{1} and k_{-1}. In addition to k_{1} and k_{-1}, the second module comprises all parameters that describe binding of E_{1} to R but no others. This special hierarchical structure is very advantageous for parameter estimation.
Measurements of the transient behavior either of the states ${\overrightarrow{y}}_{1}$ or ${\overrightarrow{y}}_{2}$ facilitate a stepwise identification of the kinetic model parameters module by module.
Taking the assumption that the association of E_{3} and E_{4} is independent of all other occurring binding processes the structure of the fifth module changes. The state E_{3}(E_{4}(*)) is not controllable any more, since the respective binding process can neither be directly nor indirectly influenced by changes in the L concentration. If ${\overrightarrow{y}}_{1}$ is the output vector of the system the output variable E_{3}(E_{4}(*)) is determined by the steady state equation of the respective ODE. The remaining four states of the fifth module are not observable and can be simply omitted. Thus, the model can be exactly reduced to ten ODEs. The situation changes if one considers the output vector ${\overrightarrow{y}}_{2}$. The choice of different output variables does not affect controllability of a system. Thus, the state E_{3}(E_{4}(*)) is still uncontrollable and the respective ODE can be replaced by its steady state equation. However, all model states are observable in this case, and thus no further equation can be eliminated. An exactly reduced model in this case would comprise fourteen ODEs.
Example 2: Domain phosphorylation
Reaction rules describing the Example depicted in Figure 2B.
R(0, *) | + | L | ⇋ | R(L, *) | k_{1}, k_{-1} |
R(0, 0) | ⇋ | R(0, P) | k_{2}, k_{-2} | ||
R(L, 0) | ⇋ | R(L, P) | k_{3}, k_{-3} | ||
R(*, P) | + | E_{1}(*) | ⇋ | R(*, E_{1}(*)) | k_{4}, k_{-4} |
E_{1}(0, 0) | ⇋ | E_{1}(0, P) | k_{5}, k_{-5} | ||
R(*, E_{1}) | ⇋ | R(*, E_{1}(P)) | k_{6}, k_{-6} | ||
E_{1}(*, P) | + | E _{2} | ⇋ | E_{1}(*, E_{2}) | k_{7}, k_{-7} |
Hierarchical transformation for the example system depicted in Figure 2B.
[R(*, *)] | = | [R(0, 0)] + [R(0, P)] + [R(0, E_{1})] + [R(0, E_{1}(P))] + [R(0, E_{2})] + [R(L, 0)] + [R(L, P)] + [R(L, E_{1})] + [R(L, E_{1}(P))] + [R(L, E_{2})] |
[E_{1}(*)] | = | [E_{1}(0)] + [E_{1}(P)] + [E_{1}(E_{2})] + [R(0, E_{1})] + [R(0, E_{1}(P))] + [R(0, E_{2})] + [R(L, E_{1})] + [R(L, E_{1}(P))] + [R(L, E_{2})] |
[E_{2}(*)] | = | [E_{1}(E_{2})] + [E_{2}(0)] + [R(0, E_{2})] + [R(L, E_{2})] |
[R(L, *)] | = | [R(L, 0)] + [R(L, P)] + [R(L, E_{1})] + [R(L, E_{1}(P))] + [R(L, E_{2})] |
[R(*, P(*))] | = | [R(0, P)] + [R(0, E_{1})] + [R(0, E_{1}(P))] + [R(0, E_{2})] + [R(L, P)] + [R(L, E_{1})] + [R(L, E_{1}(P))] + [R(L, E_{2})] |
[E_{1}(P (*)] | = | [E_{1}(P)] + [E_{1}(E_{2})] + [R(0, E_{1}(P))] + [R(0, E_{2})] + [R(L, E_{1}(P))] + [R(L, E_{2})] |
[R(L, P(*))] | = | [R(L, P)] + [R(L, E_{1})] + [R(L, E_{1}(P))] + [R(L, E_{2})] |
[R(*, E_{1}(*))] | = | [R(0, E_{1})] + [R(0, E_{1}(P))] + [R(0, E_{2})] + [R(L, E_{1})] + [R(L, E_{1}(P))] + [R(L, E_{2})] |
[E_{1}(E_{2}(*)] | = | [E_{1}(E_{2})] + [R(0, E_{2})] + [R(L, E_{2})] |
[R(L, E_{1}(*))] | = | [R(L, E_{1})] + [R(L, E_{1}(P))] + [R(L, E_{2})] |
[R(*, E_{1}(P (*)))] | = | [R(0, E_{1}(P))] + [R(0, E_{2})] + [R(L, E_{1}(P))] + [R(L, E_{2})] |
[R(L, E_{1}(P (*)))] | = | [R(L, E_{1}(P))] + [R(L, E_{2})] |
[R(*, E_{2}(*))] | = | [R(0, E_{2})] + [R(L, E_{2})] |
[R(L, E_{2}(*))] | = | [R(L, E_{2})] |
In this example, the transformed model equations can be dissected into three unidirectionally coupled modules including all five output variables, and one additional module comprising the two unobservable states [R(L, E_{1}(P))] and [R(L, E_{2})]. All model states are controllable. Thus, this example shows that the existence of all-or-none interactions facilitate significant model reductions. Although the considered system comprises the same number of molecular processes than the previously regarded one, even the complete mechanistic model already consists of a lower number of ODEs. Additionally, the system comprises two unobservable states, which allows for a further model reduction.
Homodimerization of receptors and scaffolds
Homodimerization of receptors and scaffold proteins is quite common in signal transduction networks. For instance, homodimers occur in the ErbB signaling network as described by Citri et al. [31].
Homodimerization is additionally characterized by a number of unique features having a strong impact on model reduction which justifies a separate and detailed discussion. Due to their symmetric configuration the number of distinguishable homodimers is much lower than of equally large heterodimers. If one considers a receptor monomer which forms n distinct monomeric multiprotein complexes there exist $\frac{n(n+1)}{2}$ feasible homodimers. Heterodimerization of two different receptors, which both form n monomeric species, leads to n^{2} feasible heterodimers. However, the indistinguishability of symmetric receptor dimers not only has the positive effect of reducing the number of ODEs compared to heterodimers. It also leads to non-intuitive kinetic system properties, which will be discussed below.
We will consider a receptor R with n distinct binding domains. Furthermore, we presume that R can form homodimers. These homodimers shall be decipted as R(*, ..., *).R(*, ..., *). Due to the symmetry of the dimers one cannot distinguish between R(L, 0, ..., 0).R(0, ..., 0) and R(0, ..., 0).R(L, 0, ..., 0). Hence, we reach an agreement that the receptor with more occupied domains will always be noted first.
Kinetic properties
The reason for this duplication of the k_{on} value can be elucidated considering the reaction rates. Let us take the assumption that the two concentrations [R(X_{1}, ..., X_{ n })] and [R(Y_{1}, ..., Y_{ n })] are equal. The rates for the considered two reactions comprise the terms [R(X_{1,} ..., X_{ n })]^{2} and [R(X_{1}, ..., X_{ n })]·[R(Y_{1}, ..., Y_{ n })], respectively. According to the collision theory for chemical reactions these terms are measures for the likelihood of a collision of two reactants in the system. Due to our assumption that the concentrations of both species are equal, the evaluation of both terms leads to exactly the same numerical result. However, the likelihood for the formation of a non-mirror symmetric dimer is twofold higher than that for mirror symmetric ones. This becomes apparent if one considers the collision probability for both cases. In the second case the number of molecules that may collide is twofold higher than in the first one.
The realization of process interactions either between two binding or modification processes or between dimerization and some other processes is straight forward. If dimerization has an influence on L binding, Reaction 12 will be parametrized by k_{1} and k_{-1}, while the parameters k_{2} and k_{-2} will be used for the Reactions 13 and 14. However, one still has to account for the twofold higher association constant of Reaction 13 and the twofold higher dissociation constant of Reaction 14. The neglect of these additional factors corresponds to a mutual interaction between the two ligand binding processes within a dimer.
General transformation pattern
The invertibility of the transformation matrix suggested here can again be proved using mathematical induction.
Example
Reaction rules for the considered example of EGFR dimerization.
R(0, *) | + | EGF | ⇋ | R(EGF, *) | k_{1}, k_{-1} |
R(0, *).R(0, *) | + | EGF | ⇋ | R(EGF, *).R(0, *) | 2k_{2}, k_{-2} |
R(EGF, *).R(0, *) | + | EGF | ⇋ | R(EGF, *).R(EGF, *) | k_{2}, 2k_{-2} |
R(0, X_{1}) | + | R(0, X_{1}) | ⇋ | R(0, X_{1}).R(0, X_{1}) | k_{3}, k_{-3} |
R(0, X_{1}) | + | R(0, X_{2}) | ⇋ | R(0, X_{1}).R(0, X_{2}) | 2k_{3}, k_{-3} |
R(EGF, X_{1}) | + | R(0, X_{1}) | ⇋ | R(EGF, X_{1}).R(0, X_{1}) | k_{4}, k_{-4} |
R(EGF, X_{1}) | + | R(0, X_{2}) | ⇋ | R(EGF, X_{1}).R(0, X_{2}) | 2k_{4}, k_{-4} |
R(EGF, X_{1}) | + | R(EGF, X_{1}) | ⇋ | R(EGF, X_{1}).R(EGF, X_{1}) | k_{5}, k_{-5} |
R(EGF, X_{1}) | + | R(EGF, X_{2}) | ⇋ | R(EGF, X_{1}).R(EGF, X_{2}) | 2k_{5}, k_{-5} |
R(*, 0) | ⇋ | R(*, P) | k_{6}, k_{-6} | ||
R(*, 0).R(*, 0) | ⇋ | R(*, P).R(*, 0) | 2k_{7}, k_{-7} | ||
R(*, P).R(*, 0) | ⇋ | R(*, P).R(*, P) | k_{7}, 2k_{-7} |
Hierarchical transformation for the example system.
[R(*, *).*] | = | [R(0, 0)] + [R(EGF, 0)] + [R(0, P)] + [R(EGF, P)] + 2 [R(0, 0).R(0, 0)] |
+ 2 [R(EGF, 0).R(0, 0)] + 2 [R(0, P).R(0, 0)] + 2 [R(EGF, 0).R(EGF, 0)] | ||
+ 2 [R(EGF, P).R(0, 0)] + 2 [R(EGF, 0).R(0, P)] + 2 [R(0, P).R(0, P)] | ||
+ 2 [R(EGF, P).R(EGF, 0)] + 2 [R(EGF, P).R(0, P)] + 2 [R(EGF, P).R(EGF, P)] | ||
[R(EGF, *).*] | = | [R(EGF, 0)] + [R(EGF, P)] + [R(EGF, 0).R(0, 0)] + 2 [R(EGF, 0).R(EGF, 0)] |
+ [R(EGF, P).R(0, 0)] + [R(EGF, 0).R(0, P)] + 2 [R(EGF, P).R(EGF, 0)] | ||
+ [R(EGF, P).R(0, P)] + 2 [R(EGF, P).R(EGF, P)] | ||
[R(*, *).R(*, *)] | = | [R(0, 0).R(0, 0)] + [R(EGF, 0).R(0, 0)] + [R(0, P).R(0, 0)] + [R(EGF, 0).R(EGF, 0)] |
+ [R(EGF, P).R(0, 0)] + [R(EGF, 0).R(0, P)] + [R(0, P).R(0, P)] | ||
+ [R(EGF, P).R(EGF, 0)] + [R(EGF, P).R(0, P)] + [R(EGF, P).R(EGF, P)] | ||
[R(*, P).*] | = | [R(0, P)] + [R(EGF, P)] + [R(0, P).R(0, 0)] + [R(EGF, P).R(0, 0)] + [R(EGF, 0).R(0, P)] |
+ 2 [R(0, P).R(0, P)] + [R(EGF, P).R(EGF, 0)] + 2 [R(EGF, P).R(0, P)] | ||
+ 2 [R(EGF, P).R(EGF, P)] | ||
[R(EGF, P).*] | = | [R(EGF, P)] + [R(EGF, P).R(0, 0)] + [R(EGF, P).R(EGF, 0)] + [R(EGF, P).R(0, P)] |
+ 2 [R(EGF, P).R(EGF, P)] | ||
[R(EGF, *).R(*, *)] | = | [R(EGF, 0).R(0, 0)] + 2 [R(EGF, 0).R(EGF, 0)] + [R(EGF, P).R(0, 0)] |
+ [R(EGF, 0).R(0, P)] + 2 [R(EGF, P).R(EGF, 0)] + [R(EGF, P).R(0, P)] | ||
+ 2 [R(EGF, P).R(EGF, P)] | ||
[R(*, P).R(*.*)] | = | [R(0, P).R(0, 0)] + [R(EGF, P).R(0, 0)] + [R(EGF, 0).R(0, P)] + 2 [R(0, P).R(0, P)] |
+ [R(EGF, P).R(EGF, 0)] + 2 [R(EGF, P).R(0, P)] + 2 [R(EGF, P).R(EGF, P)] | ||
[R(EGF, *).R(EGF, *)] | = | [R(EGF, 0).R(EGF, 0)] + [R(EGF, P).R(EGF, 0)] + [R(EGF, P).R(EGF, P)] |
[R(*, P).R(*, P)] | = | [R(0, P).R(0, P)] + [R(EGF, P).R(0, P)] + [R(EGF, P).R(EGF, P)] |
[R(EGF, P).R(*, *)] | = | [R(EGF, P).R(0, 0)] + [R(EGF, P).R(EGF, 0)] + [R(EGF, P).R(0, P)] |
+ 2 [R(EGF, P).R(EGF, P)] | ||
[R(EGF, *).R(*, P)] | = | [R(EGF, 0).R(0, P)] + [R(EGF, P).R(EGF, 0)] + [R(EGF, P).R(0, P)] |
+ 2 [R(EGF, P).R(EGF, P)] | ||
[R(EGF, P).R(EGF. *)] | = | [R(EGF, P).R(EGF, 0)] + 2 [R(EGF, P).R(EGF, P)] |
[R(EGF, P).R(*, P)] | = | [R(EGF, P).R(0, P)] + 2 [R(EGF, P).R(EGF, P)] |
[R(EGF, P).R(EGF, P)] | = | [R(EGF, P).R(EGF, P)] |
Due to the absence of production and degradation, the overall concentration of EGFR stays constant and the respective ODE can be eliminated. The remaining 13 transformed model equations can be dissected into three modules. The first module consists of four ODEs and describes EGF binding as well as receptor homodimerization. It comprises the model states [R(EGF, *).*], [R(*, *).R(*, *)], [R(EGF, *).R(*, *)] and [R(EGF, *).R(EGF, *)]. The second module describes receptor phosphorylation and contains six ODEs, while the remaining three ODEs for [R(*, P).R(*.P)], [R(EGF, P).R(*, P)] and [R(EGF, P).R(EGF, P)] form the third unobservable module. Since all states are controllable the model can be reduced by omitting the three unobservable states. This reduced model then comprises ten ODEs.
Generality of the method
with $A={\left(\frac{\partial \overrightarrow{f}(\overrightarrow{x},\overrightarrow{u})}{\partial \overrightarrow{x}}\right)}_{{\overrightarrow{x}}_{o},{\overrightarrow{u}}_{o}}$, $B={\left(\frac{\partial \overrightarrow{f}(\overrightarrow{x},\overrightarrow{u})}{\partial \overrightarrow{u}}\right)}_{{\overrightarrow{x}}_{o},{\overrightarrow{u}}_{o}}$ and $\overrightarrow{x}$ ∈ ℝ^{ n }, $\overrightarrow{u}$ ∈ ℝ^{ m }. If all states of the linearized model are controllable and observable this proves that at least at the considered operating point (${\overrightarrow{x}}_{o}$, ${\overrightarrow{u}}_{o}$) all states are required to accurately describe the systems behavior and that the model is not exactly reducible. The system is said to be locally controllable at the operating point if the rank of the matrix P (see Equation 9) is n [25]. Accordingly, the system is said to be locally observable if the rank of the matrix Q (see Eqnation 7) is n [25]. However, using this analysis method one has to be aware of the fact that controllable and observable nonlinear system states might loose controllability and/or observability at individual operating points [25]. A matrix rank of n proves that the considered system cannot be exactly reduced. At least at the chosen operating point all states are controllable and therefore affect the system's input/output behavior. If a model still comprises uncontrollable or unobservable states, the rank of the related matrix will be smaller than n for all considered operating points. However, if the matrix rank is lower than n this is no proof that the system is exactly reducible and still comprises uncontrollable or unobservable states. It is also possible that the system was linearized at an unpropitiously chosen operating point. Hence, it might be necessary to repeat the test at several operating points. A matrix rank smaller than n for numerous operating points highly suggests the further existence of uncontrollable or unobservable states in the nonlinear system. However, it is no proof.
After having discussed the limitations of exact reducibility in general, we also want to address the limitations of our method. A limitation of an exact reduction method is it does not facilitate the reduction of an exactly reducible model. Unfortunately, it is not possible to generally prove that our method has or has not such limitations. However, each model that has been reduced using our method can be linearized and checked for local observability and controllability. If the rank of the corresponding P and Q matrices is n all unobservable and uncontrollable system states have been eliminated.
All examples discussed in Conzelmann et al. [16] as well as the examples discussed above have been checked for further uncontrollable or unobservable states using this approach. In all cases, the reduced models proved to be minimal realizations even for varying input and output variables. However, there exists an interesting border case, in which a model can be further reduced without affecting the input/output behavior. This border case shall be discussed below.
However, this reduction is not due to the elimination of unobservable states as defined above but results from the restricted choice of initial conditions. From these considerations it can be seen that except for the case of two identical linear subsystems no example has been found for which the proposed transformation does not provide a Kalman decomposition.
Reduced order modeling of combinatorial reaction networks
In the previous section, we discussed a general and systematic method that allows for significant and exact model reductions of combinatorial reaction networks. Now, an alternative approach shall be considered that facilitates the direct generation of the exactly reduced model equations. This reduced order modeling approach is based on the close relations between controllability and observability of a model and the process interactions of the examined system.
Controllability, observability and process interactions
From the previously regarded examples it can be seen that the number of observable and controllable states highly depends on the occurring process interactions. The question is whether the qualitative information about process interactions can give clues about the observability and controllability of a reaction network, or maybe even facilitate a direct translation to reduced model equations. Controllability and observability as well as process interactions provide information about interactions within the considered system, however, at different levels of abstraction.
Controllability and observability are properties of an ODE system, and both of them characterize the ODE couplings with respect to the system inputs and outputs. All observable states exert a certain influence on at least one of the output variables. On the other hand a state is said to be controllable, if it can be influenced either directly or indirectly by one of the system's input variables.
Process interactions describe the regarded system at a different level of abstraction. However, they also provide information about which processes can influence other processes and which can be influenced by other processes. Controllability and observability are closely related to the definition of input and output variables, respectively. In accordance to this definition at the ODE level, one can also formally define input and output processes at the process level. A connection between the two abstraction levels is given by the occurrence levels we previously introduced as a state space representation for combinatorial reaction networks. These coordinates allow a direct assignment of model states to specific molecular processes.
Each occurrence level like e.g. [R(*, ..., *, E_{ i }, *..., *)] can be directly assigned to its respective process, namely E_{ i }binding to R. Analogously, occurrence levels of higher tiers like [R(E_{1}, E_{2}, E_{3}, *, ..., *)] can be linked with three different processes. All processes that are related to the chosen output variables are said to be output processes, and all processes that can be directly assigned to the input variables analogously correspond to the input processes. This direct link between model variables and processes facilitates the unique translation of all input and output variables to a set of input and output processes. Let us consider an example and presume that the concentration of E_{1} is regarded as input variable while [R(*, E_{2}, *, ..., *)] and [R(*, *, *, E_{4}, E_{5}, *, ..., *)] are output variables. In this case E_{1} binding to R is an input process, and E_{2}, E_{4} and E_{5} binding to R are output processes.
Furthermore, we can formally introduce process controllability and process observability. A process shall be called process controllable if it is either directly or indirectly influenced by one of the input processes. Analogously, a process will be called process observable if it directly or indirectly affects one of the output processes. In contrast to controllability and observability of an ODE model, the respective system properties at the process level can be analyzed in a very simple way by considering the process interaction graph. In this graph processes are regarded as nodes, while process interactions are represented as directed edges. This definition of an interaction graph is very similar to that proposed by Klamt et al. [37]. A process P is process controllable if the interaction graph comprises a directed path from one of the input processes to the process P. The same process is observable if there exists a directed path from P to one of the output processes.
A relation between the controllability and observability concepts at the different abstraction levels can also be approved. Process controllability suggests that all states that are assigned to this process are influenced and therefore controllable. Process observability on the other hand indicates that the respective occurrence level of the 1^{st} tier is observable. State variables that describe occurrence levels of higher tiers, as [R(E_{1}, E_{2}, E_{3}, *, ..., *)], are only observable if the related processes all jointly affect at least one of the output processes. Thus, we have found a way to predict whether a certain state might be observable or controllable by considering the process interactions. Note, that this technique provides a conservative appraisal which in some cases will classify states as controllable or observable although they are not.
Reduced order modeling technique
Step 1: Definition of all proteins, binding domains as well as binding and modification processes that shall be included to the model.
Example: In the considered example the model will comprise the molecules A, B, C and D with their binding domains as depicted in Figure 3A. The occurring processes are usually labelled or numbered like indicated in Figure 3A. In the example we consider eight different processes, namely binding of A to B (process 1), phosphorylation of B at different domains (processes 2, 3 and 7), binding of C to B (process 4), phosphorylation of C at two distinct domains (processes 5 and 6) as well as binding of D (process 8).
Step 2: In a second step one has to define all occurring process interactions and whether these are uni- or bidirectional. These process interactions have to be consistent with both measured kinetic data of the involved proteins and the thermodynamic constraints as discussed in Ederer et al. and Conzelmann et al. [38, 28]. These constraints for instance highly restrict the possible occurrence of unidirectional interactions. Nature can only realize unidirectional interactions at the expense of energy [28]. Since a mathematical model requires a complete definition for all interactions, mostly fragmentary knowledge has to be completed by assumptions.
Example: In Figure 3A the occurring process interactions of the example are indicated by arrows. The processes (1, 2), (1, 3), (3, 7), (4, 5) and (4, 6) are assumed to interact unidirectionally. This assumption is in accordance with the mentioned thermodynamic constraints since phosphorylation consumes energy rich ATP [28]. The processes (3, 4) and (7, 8) are regarded as all-or-none interactions which by definition are mutual ones. All other processes do not interact directly.
Step 3: The interaction pattern of the system has to be translated into an interaction graph. As mentioned before the labelled or numbered processes are nodes and the occurring interactions are represented by directed edges (arrows) pointing to the process which is influenced.
Example: The interaction graph for the considered example is depicted in Figure 3B.
Step 4: One defines input and output processes according to the considered system stimulations as well as available measurements or research interests. If for instance the external ligand concentration is considered as input the related input process is ligand binding. If one is interested in receptor phosphorylation all phosphorylation processes at the receptor have to be chosen as output processes. The goal of generating a model that accurately describes the output processes at a macroscopic level necessitates the further inclusion of all other processes that are process observable.
Example: In the example we choose process 1 an input and the processes 2, 3, 5 and 8 as output processes. The output processes are marked by grey circles in Figure 3B.
Step 5: The interaction graph can be divided into output subgraphs. An output subgraph contains all nodes from which a specific output node can be reached following the directed edges. The most simple way of finding an output subgraph is to invert the directions of all arrows in the interaction graph and to start at one output process. Following the arrows one marks all processes that are directly connected with the output process. Afterwards one marks all other processes that are directly linked with the previously found processes. This procedure is repeated until one either reaches dead end processes or the number of marked processes does not increase anymore. The so found output subgraph comprises all processes, which are process observable considering the chosen output process. If a node does not occur in any output subgraph the corresponding process cannot influence any of the output processes and can be completely omitted in the following. Finally, one has to eliminate redundant information, i.e. subgraphs which are completely comprised in other bigger subgraphs. In principle one also can analogously define input subgraphs and determine which processes are uncontrollable. In this case one does not have to invert the directions of the interaction arrows. One simply uses the original directed graph, starts at an input process and performs exactly the same procedure as described for output processes. However, uncontrollable but observable processes cannot simply be eliminated from further consideration. Uncontrollability merely allows for steady state assumptions at the ODE level.
Example: The graph shown in Figure 3B can be divided into four output subgraphs as shown in Figure 3C. In this example process six does not influence any of the considered output processes and can be omitted in the following considerations. The subgraph for output process three is completely comprised in two other subgraphs and therefore can be eliminated. In this example there only exists one input subgraph that comprises all processes.
Step 6: Each of the output subgraphs describes an autonomous signaling path, which can be modeled separately. Note, that the first five steps are independent of the used modeling approach. For instance, one can also generate stochastic models on the basis of the derived output subgraphs. However, we will focus on ODE models now. Hence, the next step is to create complete mechanistic ODE models for each subgraph. Processes not being part of a subgraph are not included in the respective model.
in which the rates r_{1} and r_{2} describe the binding of A to the scaffold protein B (process 1), and the rates r_{3} and r_{4} describe the phosphorylation of B (process 2).
Step 7: The model equations have to be transformed to new more convenient coordinates, which allow to eliminate redundant information still included in the subgraphs. This redundancy is due to the fact that some processes are comprised in several subgraphs. A suitable choice of new coordinates is given by the previously introduced occurrence levels. Note, that we highly recommend to choose always the transformation patterns discussed in the previous sections since these guarantee that the redundant information can be eliminated. However, we also want to state that there exist numerous other transformations that also would allow to eliminate the redundant information.
In this example there only occurs one further tier describing the second order occurrence levels, namely
[B (A, P)] = [B (A, P)].
Processes that are not included in the currently considered subgraph are simply omitted since they are not observable. If the sub-model still contains unobservable states these can also be eliminated at this stage of the procedure.
Step 8: The proposed transformation allows to dissect the model equations of each subgraph into modules like shown above. These modules are characterized by unidirectional communication with other modules. Processes which directly or indirectly interact mutually form one module. If some processes are included in more than one subgraph, the models of these subgraphs will contain identical modules. Multiple copies of modules can be eliminated and the remaining modules can be merged to a complete model.
Example: For instance, the transformed ODEs for the discussed smallest subgraph do have a special structure. The variables [A(*)] and [B(*, *)] are constant and equal their initial concentration. The corresponding ODEs are not required. Additionally, the ODE for [B(A, *)] does not depend on [B(*, P)] and [B(A, P)], which is due to the unidirectional process interaction between A binding to B and phosphorylation of B. Hence, the remaining three ODEs can be divided into two modules. One module only comprises the ODE for [B(A, *)], which describes the dynamics of process 1. The second module comprises the other two ODEs, which describe the dynamics of process 2. The ODEs deduced from the two remaining output subgraphs shown in Figure 3C, can be divided into six more modules as indicated in Figure 3D. Each box represents a set of ODEs. The modules are labelled with the process numbers which are described by the appropriate ODEs. Two copies of module (1) and one of module (3,4) can be eliminated here. The resulting model, which consists of only 22 ODEs, is schematically shown in Figure 3. A complete mechanistic model of the exemplified network would comprise 77 ODEs of which three can be eliminated due to mass conservation relations.
Step 9: In a last step one can take a steady state assumption for all uncontrollable states that are still part of the reduced model. Note that this step is not obligatory and in some cases can be problematic since the steady state equations have to be solved. If it is possible it is advantageous to solve these equations analytically. However, in many cases an analytic solution might be unfeasible and the steady state equations have to be solved numerically in the step of numerical integration which then necessitates a DAE (Differential Algebraic Equation system) solver.
Example: In the regarded example all states are controllable and therefore no steady state assumptions can be made.
The main advantage of the proposed method is the direct generation of a reduced, but exact, system of equations, circumventing a unsuitable large model of full combinatorial complexity. Admittedly, the number of equations that has to be set up in step six mostly include redundant information about processes which can be observed at numerous outputs. However, the absolute number of ODEs that have to be generated is usually much lower than if a complete mechanistic model is created. In the considered example one only has to set up 27 ODEs using the described procedure. A complete combinatorial model would comprise 77 states. The method has to be slightly modified if one of the output variables is a higher order occurrence level which is not contained in any of the submodels. Let us assume that one of the output variables in the example is [B(*, P(*), P(*), *)] which describes both process 2 and process 3. Since none of the three subgraphs depicted in Figure 4C comprises both processes simultaneously the quantity [B(*, P(*), P(*), *)] will not be a state of the reduced 22 ODE model. This problem can be overcome by the fusion of two subgraphs. This will necessarily increase the number of ODEs that has to be generated as well as the final size of the reduced model. However, the number of ODEs would still be smaller than 77. Furthermore, the inclusion of production and degradation into the mathematical model necessitates another extension of this method. The same holds true if the regarded system includes multifunctional protein binding domains, i.e. that certain binding domains are involved in several binding processes. Both cases shall be discussed in the following. Note, that these problems do not occur if a the method described in Section Exact Model Reduction is used.
Multifunctional protein binding domains
Multifunctional protein binding domains are domains which can recruit more than one binding partner. A typical example is the effector protein Grb2 that can either bind to several ErbB receptors as well as to the adaptor protein Shc [39]. A constellation like this can lead to problems with the reduced order modeling approach introduced above. The problem occurs if such a multifunctional binding domain is part of two or more output subgraphs as shown in Figure 3C.
The probably most simple example to illustrate the occurring problem is a scaffold protein R which provides two binding domains. Both of these domains shall recruit the effector protein E which possesses one binding domain. The binding domain of E is a multifunctional one since it can bind to both R domains. If we assume that the two binding processes of the regarded system are completely independent and that both are considered as output processes, the system can be divided into two subgraphs. These subgraphs are somehow degenerated since both only comprise a single node. According to the reduced order modeling approach both binding processes can be modeled separately. However, the problem is that the binding domain of the effector E is involved in both processes. This is a typical crosstalk situation. Since the number of effector proteins E and therefore the number of E binding domains is limited, the recruitment of E to one receptor domain reduces the concentration of unbound effector and therefore has an indirect influence on the other binding process.
One possible solution for this problem is to merge all output subgraphs that share such multifunctional binding domains. This approach has the drawback that the number of equations that have to be generated in the sixth step of the modeling procedure can be significantly increased. Alternatively, one can formulate the reaction rates for both subgraphs independently. However, all species which are simultaneously involved in both submodels have to be balanced in one joint ODE. If the first subgraph of the regarded example is translated into a reaction rate one has to consider only the rate
r_{1} = k_{1} [R(0, #)] [E] - k_{-1} [R(E, #)].
In this representation the identifier # indicates that the real scaffold protein offers further binding domains but that the resulting combinatorial complexity is neglected. The second subgraph can be described by the reaction rate
r_{2} = k_{2} [R(#, 0)] [E] - k_{-2} [R(#, E)].
This means that occurrence levels can be composed of species from both submodels like [E(*)] in the regarded example.
This simple modification or extension of the proposed modeling approach facilitates its application to a larger set of reaction systems. It will also be of great importance in modeling the crosstalk between EGF and insulin receptor discussed below.
Production and degradation
A process that has been completely neglected in the preceding considerations which, however, plays a crucial role in many real signal transduction networks is production and degradation of signaling proteins. This process increases or decreases the number of available proteins. A quite simple way of modeling production and degradation, which we will adopt here, is the assumption of a constant production rate and a degradation rate which is proportional to the concentration of the degraded species. In many cases ubiquitin is used as marker for controlled degradation as for example shown for ErbB1 receptors [31]. A still unanswered question in this context is whether the whole signaling complex is degraded or only the ErbB receptor while the associated adaptor proteins are recycled.
For the sake of simplicity we take a number of assumptions. First, the considerations shall be focused on production and degradation of the regarded receptor or scaffold protein and its complexes. The individual receptor protein R shall be produced with a constant rate and all receptor species are presumed to have a natural decay rate. All other adaptor and effector proteins are neither produced nor degraded. If a receptor complex is degraded all bound adaptor proteins shall be recycled to the cytosol. Furthermore, we take the assumption that if the receptor is marked by ubiquitination its degradation rate is modulated. This change of the degradation rate from natural decay to ordered degradation can be considered as a process interaction. Ubiquitination has a direct influence on degradation. It is quite obvious that all processes that involve one of the R binding domains are affected by the considered production and degradation.
Theoretically, degradation can be considered as a process which sets the k_{on} values of all R binding domains to zero and all k_{off} values to infinity. All other effects caused by degradation are indirect effects. Note that if one takes the assumption that a complex is degraded with all its bound adaptor proteins all processes that modify or enlarge the R complex are directly influenced. All these interactions are unidirectional ones, which can be simply introduced in the process interaction graph. Production and degradation is one additional node in this graph which is influenced by ubiquitination and affects numerous other processes.
Example: EGF and insulin receptor crosstalk
Finally, the discussed methods shall be used to generate a reduced model of EGF and insulin receptor crosstalk. We will compare a complete mechanistic description of this crosstalk and an exactly reduced version.
Model definition
In a first step it shall be defined which molecules and processes are included to the model and what assumptions are made concerning the process interactions. Since a complete mechanistic model that is still manageable shall also be generated the considerations will be limited to a small part of the real signaling network. For instance only the EGF receptor (EGFR) will be taken into account and the other three ErbB receptors shall be neglected. Similar simplifications were made by many other modelers in the past years [12, 10, 27]. In order to avoid an unmanageable combinatorial explosion of feasible EGF receptor species only two intracellular domains will be modeled. According to Schulze et al. the EGF receptor possesses among others six potential residues for Grb2 and also six residues for Shc binding [40]. Hence, we consider one binding domain for each of these two effector proteins. Concerning the insulin receptor family we will focus on the insulin receptor (IR) and exclude potential crosstalk with the insulin-like growth factor receptor (IGFR) and the insulin related receptor (IRR) [41]. Again we restrict the considerations to two intracellular IR domains, namely one for Shc and one for IRS [30].
EGFR provides an extracellular binding domain that recruits EGF [31, 42]. Furthermore, the receptor monomers can form homodimers after being activated by the ligand. This dimerization induces phosphorylation of numerous intracellular domains [43–45]. According to thermodynamic constraints [33, 28], EGF binding and receptor dimerization have to interact mutually fulfilling the Wegscheider conditions. A mutual interaction is also suggested by experimental data [46, 47].
Phosphorylation can be unidirectionally influenced like discussed for the insulin receptor by Gherzi et al. [36]. Analogously, EGFR dimerization is assumed to unidirectionally influence EGFR autophosphorylation of the regarded intracellular domains. A direct interaction between EGF binding and phosphorylation is not presumed to occur. After the two intracellular domains are phosphorylated one of them recruits Grb2 and the other Shc [40]. The interaction between receptor phosphorylation and subsequent effector binding shall be an all-or-none interaction. Furthermore, it is also known that Shc can be phosphorylated after having bound to EGFR [39]. Shc binding is thought to unidirectionally affect Shc phosphorylation. The phosphorylated Shc protein can also recruit Grb2 [39]. Grb2 possesses an additional domain which recruits the adaptor protein SOS. SOS is a guanine exchange factor (GEF) which can activate the membrane bound small G-protein Ras by effecting the exchange of GDP for GTP [48, 49].
Active RasGTP in turn initiates the MAP kinase cascade. Phosphorylated ERK which is a component of the MAP kinase cascade stimulates a serine/threonine phosphorylation of SOS resulting in dissociation of the Grb2-SOS complex [50, 49]. Thus, we take the assumption that the Grb2-SOS binding is not influenced by Grb2 association to phosphorylated EGF receptor or phosphorylated Shc. However, if SOS is phosphorylated by ERK, which is considered as additional input signal, the Grb2-SOS complex dissociates. Here we assume a mutual interaction between SOS phosphorylation and Grb2-SOS binding.
The insulin receptor consists of two constitutively dimerized monomers and is activated exclusively by ligand binding without further oligomerization [30]. Due to the dimeric structure of the insulin receptor two insulin binding domains will be included to the model. According to the thermodynamic constraints and experimental results these two domains have to interact mutually [51]. Ligand binding is assumed to unidirectionally influence the phosphorylation of the two regarded intracellular domains [36, 28]. Shc is assumed to bind with other kinetic parameters to IR than to EGFR. However, Shc phosphorylation, Grb2 binding to phosphorylated Shc etc. is parametrized like in the case of EGFR. In order to reduce the complexity of the network numerous binding domains of the scaffold IRS are neglected. The model only accounts for IRS binding to the phosphorylated insulin receptor, subsequent IRS phosphorylation and binding of the Grb2-SOS complex. In order to reduce the complexity of the model, receptor internalization and degradation is also neglected for both IR and EGFR.
Reaction rules for the considered example of EGF and insulin receptor crosstalk.
IR(0, 0, *, *) | + | Ins | ⇋ | IR(I, 0, *, *) | k_{1}, k_{-1} |
IR(0, 0, *, *) | + | Ins | ⇋ | IR(0, I, *, *) | k_{1}, k_{-1} |
IR(I, 0, *, *) | + | Ins | ⇋ | IR(I, I, *, *) | k_{2}, k_{-2} |
IR(0, I, *, *) | + | Ins | ⇋ | IR(I, I, *, *) | k_{2}, k_{-2} |
IR(0, 0, 0, *) | ⇋ | IR(0, 0, P, *) | k_{3}, k_{-3} | ||
IR(I, 0, 0, *) | ⇋ | IR(I, 0, P, *) | k_{4}, k_{-4} | ||
IR(0, I, 0, *) | ⇋ | IR(0, I, P, *) | k_{4}, k_{-4} | ||
IR(I, I, 0, *) | ⇋ | IR(I, I, P, *) | k_{5}, k_{-5} | ||
IR(*, *, P, *) | + | Shc(*) | ⇋ | IR(*, *, Shc(*), *) | k_{6}, k_{-6} |
IR(I, I, Shc(0), *) | ⇋ | IR(I, I, Shc(P), *) | k_{7}, k_{-7} | ||
IR(*, *, Shc(P), *) | + | Grb 2(*) | ⇋ | IR(*, *, Grb 2(*), *) | k_{8}, k_{-8} |
IR(*, *, Grb 2(0), *) | + | SOS(0) | ⇋ | IR(*, *, SOS(0), *) | k_{9}, k_{-9} |
IR(*, *, Grb 2(0), *) | + | SOS(P) | ⇋ | IR(*, *, SOS(P), *) | k_{10}, k_{-10} |
IR(*, *, SOS(0), *) | ⇋ | IR(*, *, SOS(P), *) | k_{11}, k_{-11} | ||
IR(0, 0, *, 0) | ⇋ | IR(0, 0, *, P) | k_{12}, k_{-12} | ||
IR(I, 0, *, 0) | ⇋ | IR(I, 0, *, P) | k_{13}, k_{-13} | ||
IR(0, I, *, 0) | ⇋ | IR(0, I, *, P) | k_{13}, k_{-13} | ||
IR(I, I, *, 0) | ⇋ | IR(I, I, *, P) | k_{14}, k_{-14} | ||
IR(*, *, *, P) | + | IRS(*) | ⇋ | IR(*, *, *, IRS(*)) | k_{15}, k_{-15} |
IR(I, I, *, IRS(0)) | ⇋ | IR(I, I, *, IRS(P)) | k_{16}, k_{-16} | ||
IR(*, *, *, IRS(P)) | + | Grb 2(*) | ⇋ | IR(*, *, *, Grb 2(*)) | k_{17}, k_{-17} |
IR(*, *, *, Grb 2(0)) | + | SOS(0) | ⇋ | IR(*, *, *, SOS(0)) | k_{9}, k_{-9} |
IR(*, *, *, Grb 2(0)) | + | SOS(P) | ⇋ | IR(*, *, *, SOS(P)) | k_{10}, k_{-10} |
IR(*, *, *, SOS(0)) | ⇋ | IR(*, *, *, SOS(P)) | k_{11}, k_{-11} | ||
Shc(0) | ⇋ | Shc(P) | k_{18}, k_{-18} | ||
Shc(P) | + | Grb 2(*) | ⇋ | Shc(Grb 2(*)) | k_{8}, k_{-8} |
Shc(Grb 2(0)) | + | SOS(0) | ⇋ | Shc(SOS(0)) | k_{9}, k_{-9} |
Shc(Grb 2(0)) | + | SOS(P) | ⇋ | Shc(SOS(P)) | k_{10}, k_{-10} |
Shc(SOS(0)) | ⇋ | Shc(SOS(P)) | k_{11}, k_{-11} | ||
Grb 2(0) | + | SOS(0) | ⇋ | Grb 2(SOS(0)) | k_{9}, k_{-9} |
Grb 2(0) | + | SOS(P) | ⇋ | Grb 2(SOS(P)) | k_{10}, k_{-10} |
Grb 2(SOS(0)) | ⇋ | Grb 2(SOS(P)) | k_{11}, k_{-11} | ||
SOS(0) | ⇋ | SOS(P) | k_{19}, k_{-19} | ||
IRS(0) | ⇋ | IRS(P) | k_{20}, k_{-20} | ||
IRS(P) | + | Grb 2(*) | ⇋ | IRS(Grb 2(*)) | k_{17}, k_{-17} |
IRS(Grb 2(0)) | + | SOS(0) | ⇋ | IRS(SOS(0)) | k_{9}, k_{-9} |
IRS(Grb 2(0)) | + | SOS(P) | ⇋ | IRS(SOS(P)) | k_{10}, k_{-10} |
IRS(SOS(0)) | ⇋ | IRS(SOS(P)) | k_{11}, k_{-11} | ||
ER(0, *, *) | + | EGF | ⇋ | ER(E, *, *) | k_{21}, k_{-21} |
ER(*, 0, *) | ⇋ | ER(*, P, *) | k_{22}, k_{-22} | ||
ER(*, P, *) | + | Shc(*) | ⇋ | ER(*, Shc(*), *) | k_{23}, k_{-23} |
ER(*, Shc(0), *) | ⇋ | ER(*, Shc(P), *) | k_{24}, k_{-24} | ||
ER(*, Shc(P), *) | + | Grb 2(*) | ⇋ | ER(*, Grb 2(*), *) | k_{8}, k_{-8} |
ER(*, Grb 2(0), *) | + | SOS(0) | ⇋ | ER(*, SOS(0), *) | k_{9}, k_{-9} |
ER(*, Grb 2(0), *) | + | SOS(P) | ⇋ | ER(*, SOS(P), *) | k_{10}, k_{-10} |
ER(*, SOS(0), *) | ⇋ | ER(*, SOS(P), *) | k_{11}, k_{-11} | ||
ER(*, *, 0) | ⇋ | ER(*, *, P) | k_{25}, k_{-25} | ||
ER(*, *, P) | + | Grb 2(*) | ⇋ | ER(*, *, Grb 2(*)) | k_{26}, k_{-26} |
ER(*, *, Grb 2(0)) | + | SOS(0) | ⇋ | ER(*, *, SOS(0)) | k_{9}, k_{-9} |
ER(*, *, Grb 2(0)) | + | SOS(P) | ⇋ | ER(*, *, SOS(P)) | k_{10}, k_{-10} |
ER(*, *, SOS(0)) | ⇋ | ER(*, *, SOS(P)) | k_{11}, k_{-11} | ||
ER(E, *, *) | + | ER(0, *, *) | ⇋ | ER_{2}(E, *, *, 0, *, *) | k_{27}, k_{-27} |
ER(0, *, *) | + | ER(0, *, *) | ⇋ | ER_{2}(0, *, *, 0, *, *) | k_{28}, k_{-28} |
ER(E, *, *) | + | ER(E, *, *) | ⇋ | ER_{2}(E, *, *, E, *, *) | k_{29}, k_{-29} |
ER_{2}(0, *, *, *, *, *) | + | EGF | ⇋ | ER_{2}(E, *, *, *, *, *) | k_{30}, k_{-30} |
ER_{2}(*, 0, *, *, *, *) | ⇋ | ER_{2}(*, P, *, *, *, *) | k_{31}, k_{-31} | ||
ER_{2}(*, Shc(0), *, *, *, *) | ⇋ | ER_{2}(*, Shc(P), *, *, *, *) | k_{32}, k_{-32} | ||
ER_{2}(*, P, *, *, *, *) | + | Shc(*) | ⇋ | ER_{2}(*, Shc(*), *, *, *, *) | k_{23}, k_{-23} |
ER_{2}(*, Shc(P), *, *, *, *) | + | Grb 2(*) | ⇋ | ER_{2}(*, Grb 2(*), *, *, *, *) | k_{8}, k_{-8} |
ER_{2}(*, Grb 2(0), *, *, *, *) | + | SOS(0) | ⇋ | ER_{2}(*, SOS(0), *, *, *, *) | k_{9}, k_{-9} |
ER_{2}(*, Grb 2(0), *, *, *, *) | + | SOS(P) | ⇋ | ER_{2}(*, SOS(P), *, *, *, *) | k_{10}, k_{-10} |
ER_{2}(*, SOS(0), *, *, *, *) | ⇋ | ER_{2}(*, SOS(P), *, *, *, *) | k_{11}, k_{-11} | ||
ER_{2}(*, *, 0, *, *, *) | ⇋ | ER_{2}(*, *, 0, *, *, *) | k_{33}, k_{-33} | ||
ER_{2}(*, *, P, *, *, *) | + | Grb 2(*) | ⇋ | ER_{2}(*, *, Grb 2(*), *, *, *) | k_{34}, k_{-34} |
ER_{2}(*, *, Grb 2(0), *, *, *) | + | SOS(0) | ⇋ | ER_{2}(*, *, SOS(0), *, *, *) | k_{9}, k_{-9} |
ER_{2}(*, *, Grb 2(0), *, *, *) | + | SOS(P) | ⇋ | ER_{2}(*, *, SOS(P), *, *, *) | k_{10}, k_{-10} |
ER_{2}(*, SOS(0), *, *, *, *) | ⇋ | ER_{2}(*, *, SOS(P), *, *, *) | k_{11}, k_{-11} |
Complete mechanistic model vs. exactly reduced model
A complete mechanistic model of the described network of EGF and insulin receptor crosstalk comprises 42,956 reactions and 5,182 ODEs. According to the assumed process interactions the complete network can be parametrized by 68 kinetic parameters which can be seen in Table 7. The exact numerical value of these parameters does not play an important role in this context. The main purpose of this model is to illustrate the possibilities offered by the new reduction methods. Hence, the model equations are normalized to relative concentrations. The overall concentration of the considered components EGFR, IR, Shc, Grb2, SOS and IRS are set to 100 percent. The kinetic parameters are chosen such that the model qualitatively shows the expected behavior. We will focus on time plots of the quantities [IR(*, *, SOS(*), *)], [IR(*, *, *, SOS(*))], [EGFR(*, SOS(*), *).*] and [EGFR(*, *, SOS(*)).*].
The complete mechanistic model can be generated by BIO NET GEN or other similar rule-based modeling tools. This example was modeled using the software tool ALC [29]. ALC allows the generation of combinatorial network models and produces output files for both MATLAB and MATHEMATICA. One simulation run with the MATLAB integrator ode15s took several hours using an Intel^{®} Xeon™ CPU with 3.06 GHz and 2 GB main memory. The simulation can be sped up by providing an analytically derived Jacobian matrix of the ODE system. All non-zero elements of the Jacobian matrix have been analytically calculated using MATHEMATICA and afterwards have been exported to MATLAB. The resulting simulation files have a size of over 13 MB, and a single simulation run still takes about half an hour. All files required to simulate the complete reaction network are provided in a ZIP-files [see Additional file 1].
Conclusion
Mathematical models of biochemical reaction networks play an increasing role in cytological research. Most of the underlying reaction networks are far too complex to facilitate an intuitive understanding. In this contribution, the focus is on ODE based dynamic modeling of receptor mediated signal transduction in mammalian cells like insulin or epidermal growth factor (EGF) signaling. These networks share some common features. Ligand binding to a receptor triggers conformational changes that facilitate receptor dimerization and/or phosphorylation of numerous residues. The subsequent formation of multiprotein signaling complexes on these receptors and their scaffolding adaptor proteins initiates a variety of signaling pathways. The main problem that occurs in modeling these networks using common modeling strategies is the enormous number of feasible multiprotein species and the high complexity of the related reaction networks. The main contribution of this work for ODE based modeling of signal transduction pathways is the extension and further development of an existing model reduction technique and the introduction of a reduced order modeling technique that allows to generate manageable reduced models accounting for the dynamic effects of combinatorial complexity.
For common in- and output signals the number of unobservable and uncontrollable model states depends on the occurring process interactions and is usually fairly high for a complete mechanistic model. The elimination of uncontrollable and unobservable state variables can be achieved by a linear and hierarchically structured state space transformation, which additionally facilitate a modularization of the model equations. Due to the enormous size of many real signaling cascades the generation of a complete mechanistic model and its subsequent reduction is not practical. An alternative approach is directly based on the process interaction pattern of the regarded system. All occurring process interactions can be integrated in an interaction graph which is subsequently dissected into independent interaction subgraphs. This exact reduced order modeling technique is used to generate a reduced model of EGF and insulin receptor crosstalk. This method allows to fairly reduce the complete mechanistic model with 5,182 ODEs to solely 87. Simulation studies show that the reduced model has exactly he same input/output behavior than the complete mechanistic model.
Thus, the results of this contribution provide new and powerful tools for dynamic modeling of combinatorial reaction networks like they occur in signal transduction. The introduced reduction techniques facilitate the generation of fairly reduced and modularized dynamic models. The modular structure of the resulting models also reduces the complexity of parameter estimation. Furthermore, the availability of an alternative reduced order modeling approach also facilitates the handling of very large and complex signaling networks. This property is of immense practical relevance since most real signaling cascades are too complex to be translated into a complete mechanistic model which is subsequently reduced.
Declarations
Acknowledgements
Special thanks to Michael Ederer and Markus Koschorreck for numerous and fruitful discussions, to Rebekka Schlatter for carefully reading the manuscript. HC, and EDG acknowledge support from the Bundesministerium für Bildung und Forschung (BMBF) within the framework of the initiatives HepatoSys, QuantPro and FORSYS. HC wants to dedicate his part of the work to his brother Jochen (born April 12, 1981 – died March 4, 2008).
Authors’ Affiliations
References
- Hlavacek WS, Faeder JR, Blinov ML, Perelson AS, Goldstein B: The complexity of complexes in signal transduction. Biotechnol Bioeng. 2003, 84 (7): 783-794.View ArticlePubMedGoogle Scholar
- Conzelmann H, Saez-Rodriguez J, Sauter T, Bullinger E, Allgöwer F, Gilles ED: Reduction of mathematical models of signal transduction networks: simulation-based approach applied to EGF receptor signalling. Syst Biol (Stevenage). 2004, 1: 159-169.View ArticleGoogle Scholar
- Koschorreck M, Conzelmann H, Ebert S, Ederer M, Gilles ED: Reduced modeling of signal transduction – A modular approach. BMC Bioinformatics. 2007, 8: 336-PubMed CentralView ArticlePubMedGoogle Scholar
- Endy D, Brent R: Modelling cellular behaviour. Nature. 2001, 409 (6818): 391-395.View ArticlePubMedGoogle Scholar
- Arkin AP: Synthetic cell biology. Curr Opin Biotechnol. 2001, 12 (6): 638-644.View ArticlePubMedGoogle Scholar
- Asthagiri AR, Lauffenburger DA: A computational study of feedback effects on signal dynamics in a mitogen-activated protein kinase (MAPK) pathway model. Biotechnol Prog. 2001, 17 (2): 227-239.View ArticlePubMedGoogle Scholar
- 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 mitogen-activated protein kinase (MAPK) and Akt pathways in heregulin-induced ErbB signalling. Biochem J. 2003, 373 (Pt 2): 451-463.PubMed CentralView ArticlePubMedGoogle Scholar
- Haugh JM, Schooler K, Wells A, Wiley HS, Lauffenburger DA: Effect of epidermal growth factor receptor internalization on regulation of the phospholipase C-gamma1 signaling pathway. J Biol Chem. 1999, 274 (13): 8958-8965.View ArticlePubMedGoogle Scholar
- Haugh JM, Wells A, Lauffenburger DA: Mathematical modeling of epidermal growth factor receptor signaling through the phospholipase C pathway: mechanistic insights and predictions for molecular interventions. Biotechnol Bioeng. 2000, 70 (2): 225-238.View ArticlePubMedGoogle Scholar
- 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): 30169-30181.View ArticlePubMedGoogle Scholar
- 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: 306-320.View ArticlePubMedGoogle Scholar
- Schoeberl B, Eichler-Jonsson 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): 370-375.View ArticlePubMedGoogle Scholar
- Sedaghat AR, Sherman A, Quon MJ: A mathematical model of metabolic insulin signaling pathways. Am J Physiol Endocrinol Metab. 2002, 283 (5): E1084-E1101.View ArticlePubMedGoogle Scholar
- Eissing T, Conzelmann H, Gilles ED, Allgöwer F, Bullinger E, Scheurich P: Bistability analyses of a caspase activation model for receptor-induced apoptosis. J Biol Chem. 2004, 279 (35): 36892-36897.View ArticlePubMedGoogle Scholar
- Faeder JR, Blinov ML, Goldstein B, Hlavacek WS: Combinatorial complexity and dynamical restriction of network flows in signal transduction. Syst Biol (Stevenage). 2005, 2: 5-15.View ArticleGoogle Scholar
- Conzelmann H, Saez-Rodriguez J, Sauter T, Kholodenko BN, Gilles ED: A domain-oriented approach to the reduction of combinatorial complexity in signal transduction networks. BMC Bioinformatics. 2006, 7: 34-PubMed CentralView ArticlePubMedGoogle Scholar
- Lok L, Brent R: Automatic generation of cellular reaction networks with Moleculizer 1.0. Nat Biotechnol. 2005, 23: 131-136.View ArticlePubMedGoogle Scholar
- Blinov ML, Faeder JR, Yang J, Goldstein B, Hlavacek WS: 'On-the-fly' or 'generate-first' modeling?. Nat Biotechnol. 2005, 23 (11): 1344-5. author reply 1345View ArticlePubMedGoogle Scholar
- Blinov ML, Faeder JR, Goldstein B, Hlavacek WS: BioNetGen: software for rule-based modeling of signal transduction based on the interactions of molecular domains. Bioinformatics. 2004, 20 (17): 3289-3291.View ArticlePubMedGoogle Scholar
- Faeder JR, Hlavacek WS, Reischl I, Blinov ML, Metzger H, Redondo A, Wofsy C, Goldstein B: Investigation of early events in Fcϵ RI-mediated signaling using a detailed mathematical model. J Immunol. 2003, 170 (7): 3769-3781.View ArticlePubMedGoogle Scholar
- 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): 136-151.View ArticlePubMedGoogle Scholar
- Borisov NM, Markevich NI, Hoek JB, Kholodenko BN: Signaling through receptors and scaffolds: independent interactions reduce combinatorial complexity. Biophys J. 2005, 89 (2): 951-966.PubMed CentralView ArticlePubMedGoogle Scholar
- Borisov NM, Markevich NI, Hoek JB, Kholodenko BN: Trading the micro-world of combinatorial complexity for the macro-world of protein interaction domains. Biosystems. 2006, 83 (2–3): 152-166.PubMed CentralView ArticlePubMedGoogle Scholar
- Pawson T, Nash P: Assembly of cell regulatory systems through protein interaction domains. Science. 2003, 300 (5618): 445-452.View ArticlePubMedGoogle Scholar
- Isidori A: Nonlinear Control Systems. Springer. 2002Google Scholar
- Weiss L, Kalman R: Contributions to linear systems theory. Int J Eng Sc. 1965, 3: 141-171.View ArticleGoogle Scholar
- Blinov ML, Yang J, Faeder JR, Hlavacek WS: Graph theory for rule-based modeling of biochemical networks. Lect Notes Comput Sci. 2006, 4230: 89-106.View ArticleGoogle Scholar
- Conzelmann H, Gilles ED: Functional Proteomics: Methods and Protocols, Humana Press. 2008, 557-576. chap. Dynamic pathway modeling of signal transduction networks – A domain-oriented approachGoogle Scholar
- Koschorreck M, Gilles ED: ALC: automated reduction of rule-based models. 2008Google Scholar
- Ottensmeyer FP, Beniac DR, Luo RZ, Yip CC: Mechanism of transmembrane signaling: insulin binding and the insulin receptor. Biochemistry. 2000, 39 (40): 12103-12112.View ArticlePubMedGoogle Scholar
- Citri A, Yarden Y: EGF-ErbB signalling: towards the systems level. Nat Rev Mol Cell Biol. 2006, 7 (7): 505-516.View ArticlePubMedGoogle Scholar
- White MF: The IRS-signalling system: a network of docking proteins that mediate insulin action. Mol Cell Biochem. 1998, 182 (1–2): 3-11.View ArticlePubMedGoogle Scholar
- Heinrich R, Schuster S: The regulation of cellular systems. 1996, Chapman & HallView ArticleGoogle Scholar
- Odaka M, Kohda D, Lax I, Schlessinger J, Inagaki F: Ligand-binding enhances the affinity of dimerization of the extracellular domain of the epidermal growth factor receptor. J Biochem. 1997, 122 (1): 116-121.View ArticlePubMedGoogle Scholar
- Lemmon MA, Bu Z, Ladbury JE, Zhou M, Pinchasi D, Lax I, Engelman DM, Schlessinger J: Two EGF molecules contribute additively to stabilization of the EGFR dimer. EMBO J. 1997, 16 (2): 281-294.PubMed CentralView ArticlePubMedGoogle Scholar
- Gherzi R, Andraghetti G, Versari G, Cordera R: Effect of insulin receptor autophosphorylation on insulin receptor binding. Mol Cell Endocrinol. 1986, 45 (2–3): 247-252.View ArticlePubMedGoogle Scholar
- Klamt S, Saez-Rodriguez J, Lindquist JA, Simeoni L, Gilles ED: A methodology for the structural and functional analysis of signaling and regulatory networks. BMC Bioinformatics. 2006, 7: 56-PubMed CentralView ArticlePubMedGoogle Scholar
- Ederer M, Gilles ED: Thermodynamically feasible kinetic models of reaction networks. Biophys J. 2007, 92 (6): 1846-1857.PubMed CentralView ArticlePubMedGoogle Scholar
- Sakaguchi K, Okabayashi Y, Kido Y, Kimura S, Matsumura Y, Inushima K, Kasuga M: Shc phosphotyrosine-binding domain dominantly interacts with epidermal growth factor receptors and mediates Ras activation in intact cells. Mol Endocrinol. 1998, 12 (4): 536-543.View ArticlePubMedGoogle Scholar
- Schulze WX, Deng L, Mann M: Phosphotyrosine interactome of the ErbB-receptor kinase family. Mol Syst Biol. 2005, 1: 2005.0008Google Scholar
- Ward CW, Lawrence MC, Streltsov VA, Adams TE, McKern NM: The insulin and EGF receptor structures: new insights into ligand-induced receptor activation. Trends Biochem Sci. 2007, 32 (3): 129-137.View ArticlePubMedGoogle Scholar
- Pinkas-Kramarski R, Alroy I, Yarden Y: ErbB receptors and EGF-like ligands: cell lineage determination and oncogenesis through combinatorial signaling. J Mammary Gland Biol Neoplasia. 1997, 2 (2): 97-107.View ArticlePubMedGoogle Scholar
- Lemmon MA, Schlessinger J: Regulation of signal transduction and signal diversity by receptor oligomerization. Trends Biochem Sci. 1994, 19 (11): 459-463.View ArticlePubMedGoogle Scholar
- Jiang G, Hunter T: Receptor signaling: when dimerization is not enough. Curr Biol. 1999, 9 (15): R568-R571.View ArticlePubMedGoogle Scholar
- Schlessinger J: Cell signaling by receptor tyrosine kinases. Cell. 2000, 103 (2): 211-225.View ArticlePubMedGoogle Scholar
- Burgess AW, Cho HS, Eigenbrot C, Ferguson KM, Garrett TPJ, Leahy DJ, Lemmon MA, Sliwkowski MX, Ward CW, Yokoyama S: An open-and-shut case? Recent insights into the activation of EGF/ErbB receptors. Mol Cell. 2003, 12 (3): 541-552.View ArticlePubMedGoogle Scholar
- Garrett TPJ, McKern NM, Lou M, Elleman TC, Adams TE, Lovrecz GO, Zhu HJ, Walker F, Frenkel MJ, Hoyne PA, Jorissen RN, Nice EC, Burgess AW, Ward CW: Crystal structure of a truncated epidermal growth factor receptor extracellular domain bound to transforming growth factor alpha. Cell. 2002, 110 (6): 763-773.View ArticlePubMedGoogle Scholar
- Li N, Batzer A, Daly R, Yajnik V, Skolnik E, Chardin P, Bar-Sagi D, Margolis B, Schlessinger J: Guanine-nucleotide-releasing factor Sos1 binds to Grb2 and links receptor tyrosine kinases to Ras signalling. Nature. 1993, 363 (6424): 85-88.View ArticlePubMedGoogle Scholar
- Chen D, Waters SB, Holt KH, Pessin JE: SOS phosphorylation and disassociation of the Grb2-SOS complex by the ERK and JNK signaling pathways. J Biol Chem. 1996, 271 (11): 6328-6332.View ArticleGoogle Scholar
- Waters SB, Chen D, Kao AW, Okada S, Holt KH, Pessin JE: Insulin and epidermal growth factor receptors regulate distinct pools of Grb2-SOS in the control of Ras activation. J Biol Chem. 1996, 271 (30): 18224-18230.View ArticlePubMedGoogle Scholar
- Wanant S, Quon M: Insulin receptor binding kinetics: modeling and simulation studies. J Theor Biol. 2000, 205 (3): 355-64.View ArticlePubMedGoogle Scholar
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.