Reduced modeling of signal transduction – a modular approach

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 rule-based set-up of mechanistic model equations. In many cases these models can be reduced by an exact domain-oriented 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.

As there is no synthesis or degradation of proteins there are three balance relations. Each of them could replace one differential equation.

Reduced model
Nomenclature of state variables is analogous to the small example system and specified as follows: The unmodified receptor is denoted as ROOO. The first O specifies the unoccupied ligand binding site, the second one the unphosphorylated regulatory phosphorylation site and the third one the unphosphorylated effector binding site. Phosphorylation of a phosphorylation site leads to P on the corresponding position, ligand binding leads to L on the corresponding position. x is the sum of all receptor species with phosphorylated binding site for E, xb is the sum of all receptor species with phosphorylated binding site and bound E, x2 is the sum of all E-species with phosphorylated binding site for F , x2b is the sum of all E-species with phosphorylated binding site and bound F (in this case only the specie XEF ). Free species are referred to as E, EP and F . Note that EP is a lumped state as it represents free and phosphorylated E with bound and unbound F . The sum of all catalytically active receptors is denoted as R activ . All receptors with regulatory phosphorylation are assumed to be catalytically active and phosphorylate receptor bound E. Only receptor bound E is phosphorylated. Rates r 1 -r 4 describe ligand binding, r 5 -r 8 regulatory phosphorylation and r 9 -r 12 binding site phosphorylation on the receptor. Binding of E, phosphorylation of receptor bound E and binding of EP is described by rates r 13 , r 14 and r 15 , respectively. Dephosphorylation of free EP and binding of F to its phosphorylated binding site on E is described by rates r 16 and r 17 . As there is no synthesis or degradation of proteins there are three balance relations. They can replace one differential equation in each layer.
Receptor layer

Relations
The basic assumption of the layer based approach is the existence of relations between the species of the detailed model. These relations arise from independency assumptions between processes. For the larger example system there exist ten independent relations that explain the difference of ten between the detailed and reduced formalisms in the number of equations.
One might wonder, why there are no equations like This equations would only hold if phosphorylation of E was independent of E-binding to the receptor or if autophosphorylation of the receptor was independent of ligand binding. Then additional reduction were possible.

Transformations
This chapter deals with conversion between states of detailed and reduced formalisms. The states E and F are equivalent in both formalisms. For clarity reasons the states Ep and EF of the detailed model are referred to as dEp and dEF . The state EP of the reduced model is referred to as rEP . With the relations above the transformation is invertible.

Transformation of selected states
There are two general possibilities to approximate the states of the detailed model from the states of the reduced model. The first is to sum up the corresponding states of the detailed model for each lumped state. Micro-states have their direct counterpart in both formalisms, so there is a one-to-one relation for each micro-state. Then all relations between the species that arise from independency assumptions, as discussed in the manuscript, are included.
DOOEp is also contained in RXEP , so we take RXEP as a starting point. The fraction of phosphorylated E that is not bound by F is x2−x2b x2 . The fraction of phosphorylated binding sites on the receptor with no regulatory phosphorylation or bound ligand is ROOP x .
The result is the same as above: · ROOP x A third possibility is shown in the transformation equations. All possibilities give identical results. This can be easily verified by inserting x2 − x2b = XEp. DXXF that is used in the optimization study represents all species with F bound to the receptor via E. DXXF is contained in XEF . The fraction of phosphorylated E that is bound to the receptor is RXEP x2 .

RXEP x2
Optimization studies Optimizations studies for this system were curried out with the lsqnonlin routine from Matlab (The Math-Works, Inc.), version 7.2.0.294. This routine minimizes the sum of squares of the passed values. Optimization parameters were optimset('MaxIter', 100000, 'MaxFu-nEvals', 100000, 'TolFun', 1e-8, 'TolX', 1e-8), which is more stringent than the standard values. The choice of this criterion is a tradeoff between accuracy and time consumption. With the choice above, optimization took about ten CPU-days. Simulation time was 300 s. Ten data points output per simulation second were processed and passed to the minimization algorithm. The values passed to the lsqnonlin routine (out) were computed as where δ is the difference between the value of the reduced model and the value computed from the states of the detailed model at the corresponding time. This form of scoring function was chosen as standard routines minimize sums of squares. However, in this special case maximization had to be performed. The positive constant θ had to be introduced to prevent division by zero if the values are identical. As a tradeoff between preventing too high values for small differences and a domination of the effect of θ in the minimization, θ was chosen to be 0.1 in the present study. Start values for parameters and initial conditions were taken from Table 1. To vary inside physiologic ranges each parameter except the regulatory factors is allowed to be between hundred times the original value and a hundredth of this value. Therefore, variation within four degrees of magnitude around literature values was allowed. The regulatory factors describe the inhibition of autophosphorylation in the absence of bound insulin (f ins ) or in the presence of regulatory phosphorylation (f reg ). They are defined to be in the interval [0,1], variation was allowed in the interval [10 −4 ,0.5]. It is assumed, that receptor dephosphorylation is not subject to any receptor modification: k −2 = k −3 = k −4 = k −5 and that effectors E and F bind with the same kinetic constants to their phosphorylated binding sites: k 6 = k 8 and k −6 = k −8 . The resulting worst case parameters for both scenarios can be found in Table 1. Figures 2 and 3 show simulation results for both worst case scenarios. Reduction quality is still very high, the error is surely within the range of measurement errors.   Table 1. The axis of abscissae is given in s, the axis of ordinates in nM .  Table 1. The axis of abscissae is given in s, the axis of ordinates in nM . State variables that are plotted in the upper rows show very fast dynamics. In grayscale printouts their curves may not be visible as they lie very close to coordinate axes. Table 1: Parameters and initial conditions Initial values were 40 nM for DOOO, 250 nM for E and 50 nM for F . All other initial values were 10 −20 nM. Initial values for the reduced model were computed to be consistent. L is set to 100 nM, which is a typical insulin concentration for stimulation. k 3 defines autophosphorylation in the presence of ligand and the absence of regulatory phosphorylation. k 2 = k 3 · f ins defines autophosphorylation in the absence of ligand and regulatory phosphorylation. k 4 = k 3 · f ins · f reg defines autophosphorylation in the absence of ligand but with regulatory phosphorylation. k 5 = k 3 · f reg defines autophosphorylation with ligand and regulatory phosphorylation. It is assumed that dephosphorylation is independent from other receptor modifications: k −2 = k −3 = k −4 = k −5 . It is further assumed that E and F bind with the same kinetic constants to their respective binding sites: k 6 = k 8 and k −6 = k −8