Robust simplifications of multiscale biochemical networks

Abstract

Background: Cellular processes such as metabolism, decision making in development and differentiation, signalling, etc., can be modeled as large networks of biochemical reactions. In order to understand the functioning of these systems, there is a strong need for general model reduction techniques allowing to simplify models without loosing their main properties. In systems biology we also need to compare models or to couple them as parts of larger models. In these situations reduction to a common level of complexity is needed.

Results: We propose a systematic treatment of model reduction of multiscale biochemical networks. First, we consider linear kinetic models, which appear as "pseudo-monomolecular" subsystems of multiscale nonlinear reaction networks. For such linear models, we propose a reduction algorithm which is based on a generalized theory of the limiting step that we have developed in 1. Second, for non-linear systems we develop an algorithm based on dominant solutions of quasi-stationarity equations. For oscillating systems, quasi-stationarity and averaging are combined to eliminate time scales much faster and much slower than the period of the oscillations. In all cases, we obtain robust simplifications and also identify the critical parameters of the model. The methods are demonstrated for simple examples and for a more complex model of NF-kappaB pathway.

Conclusion: Our approach allows critical parameter identification and produces hierarchies of models. Hierarchical modeling is important in "middle-out" approaches when there is need to zoom in and out several levels of complexity. Critical parameter identification is an important issue in systems biology with potential applications to biological control and therapeutics. Our approach also deals naturally with the presence of multiple time scales, which is a general property of systems biology models.

Figure 1

Two simple examples of exactly solvable linear kinetics. a) non-branching network without cycles. b) network with a unique sink which is a cycle. On the left, ϕ(i) map is shown for the network a). The order of kinetics parameters is shown both by integer numbers (ranks) and the thickness of arrows (faster reactions are thicker).

Figure 2

Example of calculation of the dominant approximation for a linear separated reaction network shown (1). See the text for the details. The order of kinetics parameters is shown both by integer numbers (ranks) and the thickness of arrows (faster reactions are thicker).

Figure 3

Lipniacki's model a) Testing quasistationarity: nonreduced trajectories (solid), quasi-stationarity trajectories (crosses). b) Trajectories of models in the hierarchy. c) Cytoplasmatic part of the signalling mechanism: terminal species (blue), intermediate species quasi-stationary (pink) non-oscillating (green), simple submechanisms (blue). This part of the network contains three critical parameters for the damping time. Sustained oscillations were obtained by decreasing the constant k₃ ten times with respect to the value used in [53] (equivalently, this can be obtained by decreasing k₉, or by increasing k₄).

Figure 4

Log-log sensitivity of the damping time and of the period of the oscillations with respect to variations of different parameters of the modelℳ(5, 8, 15). The parameters are multiplied by a scale s ∈ (1/50, 50). The log(timescales) are represented as functions of log(s). Period and damping time are not represented on intervals of parameter values where oscillations are over-damped (the ratio of the damping time to the period is smaller than 1.75). Damping time is infinite and not represented for intervals of parameter values where oscillations are self-sustained. The latter intervals are limited by Hopf bifurcations where the damping time diverges.

Figure 5

Correspondence between the parameters of the modelsℳ(14, 25, 28) and ℳ(5, 8, 15). Parameters of the first model are gathered into monomials that are parameters of the reduced model. The integers on the arrows connecting parameters represent the corresponding powers of the parameters in the monomial. The critical monomials are connected to the property on which they act upon (here sustained oscillations). Thus, an increase of k_21p1 = k₂k₉$k_4^{-1}$ favors significantly the oscillations.

Figure 6

Complete model ℳ(39, 65, 90) (left, top). Intermediate mechanisms for 1) Production module of p65; 2) Min-funnel production of p50:p65@csl; 3) Production module of p50.

Figure 7

Model comparison a) Trajectories of various species for the model M (39, 65, 90); quasi-stationary species have concentrations in the lower cluster. b) Production rates of mRNAIκB for two models having the same reactions and species, differing only by one kinetic law. c) Trajectories (signal applied at t = 20). Notice the different behavior of IkBa@csl in ℳ(14, 25, 28).