Institute of Systems Biology, Ltd, 15 Detskiy proezd, Novosibirsk 630090 Russia

Design Technological Institute of Digital Techniques, The Siberian Branch of The Russian Academy of Sciences, 6 Acad. Rzhanov Str, Novosibirsk 630090 Russia

Institut Curie, 26 rue d’Ulm, F-75248 Paris, France

INSERM U900, Paris F-75248 France

Mines ParisTech, Fontainebleau F-77300 France

Institute of Cytology and Genetics, The Siberian Branch of The Russian Academy of Sciences, 10 Acad. Lavrentyev Ave, Novosibirsk 630090 Russia

Abstract

Background

Many mathematical models characterizing mechanisms of cell fate decisions have been constructed recently. Their further study may be impossible without development of methods of model composition, which is complicated by the fact that several models describing the same processes could use different reaction chains or incomparable sets of parameters. Detailed models not supported by sufficient volume of experimental data suffer from non-unique choice of parameter values, non-reproducible results, and difficulty of analysis. Thus, it is necessary to reduce existing models to identify key elements determining their dynamics, and it is also required to design the methods allowing us to combine them.

Results

Here we propose a new approach to model composition, based on reducing several models to the same level of complexity and subsequent combining them together. Firstly, we suggest a set of model reduction tools that can be systematically applied to a given model. Secondly, we suggest a notion of a minimal complexity model. This model is the simplest one that can be obtained from the original model using these tools and still able to approximate experimental data. Thirdly, we propose a strategy for composing the reduced models together. Connection with the detailed model is preserved, which can be advantageous in some applications. A toolbox for model reduction and composition has been implemented as part of the BioUML software and tested on the example of integrating two previously published models of the CD95 (APO-1/Fas) signaling pathways. We show that the reduced models lead to the same dynamical behavior of observable species and the same predictions as in the precursor models. The composite model is able to recapitulate several experimental datasets which were used by the authors of the original models to calibrate them separately, but also has new dynamical properties.

Conclusion

Model complexity should be comparable to the complexity of the data used to train the model. Systematic application of model reduction methods allows implementing this modeling principle and finding models of minimal complexity compatible with the data. Combining such models is much easier than of precursor models and leads to new model properties and predictions.

Background

Systems biology aims to study complex interactions in living systems and focuses on analysis and modeling their properties. Mathematical modeling provides several ways to describe biological processes based on experimental information of different kind. However, creation of detailed models not supported by enough experimental data often makes their analysis and interpretation difficult

Model reduction is a well-established technique in many fields of biochemical research and engineering. It has been used for many years in chemical kinetics (for reviews, see

In our investigation, we used the principles of model reduction to construct reasonably accurate minimal size approximations of two different models describing the CD95 signaling pathways

The authors of the models have evaluated concentration changes of major apoptotic molecules by Western blot analysis and represented them as relative values at given time points. Using the systematic methodology

Methods

Model reduction

Mathematical modeling of biological processes based on the classical theory of chemical kinetics assumes that a model consists of a set of species _{
1,…,}
_{
m
}) associated with a set of variables _{1}(_{
m
}(^{+} and representing their concentrations, and a set of biochemical reactions with rates _{1}(_{
n
}(

Here ^{
SS
} is a steady state of the system (1) if

Model reduction implies transformation of the ODE system to another one with smaller number of equations without affecting dynamics of variables _{1}(_{
s
}(_{
i
}
^{
exp
}(_{
ij
}) at given times _{
ij
}, _{
i
}, where _{
i
} is number of such points for the concentration _{
i
}(

where _{
min
}/ω_{
i
} is used to make all concentration trajectories to have similar importance.

Reducing kinetic model is possible when some quantities are much smaller than other quantities and can be neglected. Usually, this implies some qualitative relations (much bigger, much smaller) between model parameters. When these relations satisfy certain rules, we can approximate the detailed model by a simpler one. If the parameters are (approximately) determined, as in our case, we find an approximation specific for a region of the parameter space. However, if we want to investigate the model behavior for the entire space, we need to decompose it into regions characterized by asymptotically different behaviors of the dynamical systems. Afterwards, a specific reduction should be performed for each region (Figure

Schematic representation of the model reduction technique

**Schematic representation of the model reduction technique.** A model consisting of four species (_{A}(_{B}(_{r1}(_{r2}(_{A}(_{B}(_{r1}(_{r2}(

We define the

Reduction of mathematical model complexity is achievable by different methods

**(MR1)**
_{1} is much slower than reaction ^{-2}. When such reactions do not affect experimental dynamics of species, we neglect them. Note, that in some cases, it is sufficient to consider ^{−1} or even

**(MR2)**
^{
s
}(^{
f
}(

where ^{
f
} · ^{
f
}(

**(MR3)**

with kinetic rates _{1} = _{2} and

where _{1} = _{2}. Note, that in this example lumped variables

**(MR4)**

then we can replace one of them with another in all kinetic laws of the model.

**(MR5)**

where the enzyme concentration is a dynamic variable allowing to use the same kinetics in different regions of the phase space. If _{
S
}(_{
S
}(_{
r
}(_{
E
}(

**(MR6)**
_{1} + _{2} →

When

One way to perform the model reduction is to apply the foregoing methods in the numerical order (Figure

Flow chart of the model reduction

**Flow chart of the model reduction.** The methods of model reduction are presented in the order of their application. The last two methods are of the same type, so this is a choice of the systems biologist which one to use first.

Model analysis and comparison

The function χ^{2} is defined by the formula (3) with weights 1/
_{
ij
}
^{2} (instead of _{
min
}/_{
i
}) _{
ij
} of experimental values _{
i
}
^{
exp
}(_{
ij
}) are calculated using smoothing spline

For _{
i
} we assumed _{
i
}(_{
i,j+k
}) = _{
i
}
^{
exp
}(_{
i,j+k
}) = 0.

Models can be compared using the Akaike criterion, if they approximate the same experimental data. When the

where _{
exp
} determines the number of experimental points.

_{
i
}
^{
ss
}, _{
j
},

It is useful to compute the mean log sensitivity for all fitted species:

where _{
p
} is the number of all parameters |

Modular modeling

Considering the mathematical models of CD95 signaling pathways

Results

Preliminary analysis, problems and inconsistencies in the precursor models

The mathematical model by Bentele _{
aa
} intended for quantification of “apoptotic activity” caused by active caspases-3, -6 and -7 and calculated by the formula

where _{36act
}, _{7act
}, _{367act
} ∈ _{
aa
} specifies a process of degradation in the model introduced as an exponential decay function _{
degr
}(_{
aa
}). This function is defined by the formula _{
d
} · _{
aa
}
^{2} + _{
ds
} with _{
d
},_{
ds
} ∈ _{
d
} · _{
aa
}
^{2} for all other molecules besides cPARP (cleaved PARP), where _{
degr
}(_{
aa
}) is constant. All species in the model are degraded with the exception of cytochrome

The model by Bentele

**The model by Bentele ****and results of its reduction. A.** The original model decomposed into modules according to three steps of apoptosis: activation of caspases-8 and -9 and inactivation of PARP. The species retained after the model reduction are colored. **B.** The modular view of the model. The dashed connections were deleted during the model reduction. **C.** The minimal reduced model. **D.** The graphical notation used for representation of the models **A**-**C**.

In total, the model contains 80 reactions (Table _{
aa
} production specified above and two reactions of cytochrome

where _{
trigger
} is the time when tBid concentration reaches a value of Bcl-2/Bcl-XL, _{
release
} is the start time of release and _{
contr
} is the contraction coefficient.

**Original model**

**Reduced model**

**№**

**Reactions (Kinetics)**

**Rates (nM/min)**

**№**

**Reactions (Kinetics)**

Prototypes of parameter notations in the original Bentele’s model are listed in Additional file

**
Activation of caspase-8 induced by CD95
**

br1

CD95L + CD95R → СD95R:СD95L (MA)

10^{1}/10^{0}

br1*

CD95L → DISC (_{
LR
} · _{
CD95R
} · _{
CD95L
}, _{
CD95R
} is determined by the formula (14))

br2

FADD + СD95R:СD95L → DISC (MA)

10^{1}/10^{0}

br3

pro8 + DISC → DISC:pro8 (MA)

10^{0}/10^{-1}

br4

pro8 + DISC:pro8 → DISC:pro8_{2} (MA)

10^{0}/10^{-1}

br5

DISC:pro8_{2} → DISC:p43/p41 (MA)

10^{0}/10^{-1}

br2*

pro8 –DISC → casp8 (_{
DISC_pro8} · _{
pro8} · _{
DISC
})

br6

DISC:p43/p41 → casp8 + DISC (MA)

10^{0}/10^{-1}

**
Activation of caspase-8 by caspase-3
**

br7

pro6 -casp3 → casp6 (M-M)

10^{0}/10^{-1}

br3*

pro8 –casp3 → casp8 (_{38} · _{
casp3})

br8

pro8 -casp6 → casp8 (M-M)

10^{-1}/10^{0}

**
Inhibition of the DISC complex
**

br9

DISC + cFLIPL → DISC:cFLIPL (MA)

10^{1}/10^{0}

br4*

cFLIP + 2·DISC + pro8 → DISC:FLIP:pro8 (_{
DISC_FLIP
} · _{
FLIP
} · _{
DISC
})

br10

pro8 + DISC:cFLIPL → DISC:cFLIPL:pro8 (MA)

10^{0}

br11

DISC + cFLIPS → blocked DISC (MA)

10^{1}/10^{0}

br12

DISC:pro8 + cFLIPL → DISC:cFLIPL:pro8 (MA)

10^{-2}/10^{-3}

br5*

DISC:FLIP:pro8 → p43/p41 (_{
DFp8} · _{
DISC:FLIP:pro8})

br13

DISC:pro8 + cFLIPS → blocked DISC (MA)

10^{-2}/10^{-3}

br14

DISC:cFLIPL:pro8 → p43/p41 + blocked DISC (MA)

10^{-1}

**
Activation of caspase-9 triggered by cytochrome C
**

br15

Cyt C stored → Cyt C (_{
release
}(_{
CytCstored
})

10^{2}

br6*

pro9 → cleavage (the formula (16))

br16

Apaf-1 + Cyt C → Cyt C:Apaf-1 (MA)

10^{2}

br17

Cyt C:Apaf-1 → Apaf-1 + Cyt C (MA)

10^{1}

br18

pro9 + Cyt C:Apaf → Apop (MA)

10^{0}

br19

Apop -casp3 → casp9 + Cyt C:Apaf-1

10^{0}

(M-M)

br20

Apop → casp9 + Cyt C:Apaf-1 (MA)

10^{-2}

**
Activation of Bid
**

br21

pro2 –casp3 → casp2 (M-M)

10^{0}

br7*

Bid –casp8 → tBid (

br22

Bid –casp8 → tBid (M-M)

10^{0}

br23

Bid –casp2 → tBid (M-M)

10^{-1}/10^{0}

br8*

pro2 –casp3 → cleavage (

**
Blocking of IAP by Smac
**

br24

Smac stored → Smac (_{
release
}(_{
Smac stored
})

10^{2}

–

br25

Smac + IAP → IAP:Smac (MA)

10^{-5}

br26

IAP:Smac → Smac + IAP (MA)

10^{-5}

**
Activation of caspases-3 and −7
**

br27

pro3 –casp8 → casp3 (M-M)

10^{1}/10^{0}

br9*

pro3 –casp8 → casp3 (M-M)

br28

pro3 –casp9 → casp3 (M-M)

10^{-3}/10^{0}

br29

pro7 –casp8 → casp7 (M-M)

10^{0}/10^{-1}

br10*

pro7 –casp8 → cleavage (

br30

pro7 –casp9 → casp7 (M-M)

10^{-5}/10^{-3}

**
PARP inactivation
**

br31

PARP -casp3 → cPARP (М-М)

10^{-1}

br11*

PARP -casp3 → cPARP (

br32

PARP –casp7 → cPARP (М-М)

10^{0}/10^{-1}

**
Inhibition of caspases-3, -7 and −9
**

br33

casp3 + IAP → casp3:IAP (МA)

10^{0}

br12*

casp3 + IAP → inhibition (MA)

br34

casp7 + IAP → casp7:IAP (МA)

10^{-1}

br35

casp9 + IAP → casp9:IAP (МA)

10^{-3}

br36

casp3:IAP → casp3 + IAP (МA)

10^{-2}

br37

casp7:IAP → casp7 + IAP (МA)

10^{-4}/10^{-3}

br38

casp9:IAP → casp9 + IAP (МA)

10^{-5}/10^{-4}

**
Production of the virtual variable
**

br39

→ _{
aa
} (1 − _{
aa
}/_{367act
} + (1 − _{
aa
}) · (_{36act
} · _{
casp3} + _{36act
} · _{
casp6} + _{7act
} · _{
casp7}))

10^{0}

br13*

→ _{
aa
} (

Notations of parameters used in the paper and in the original models. **Table 2S.** Table of all reactions (excepting degradation) and rate constants in the composite model. **Table 3S.** Table of degraded species and kinetic laws of degradation in the composite model. **Table 4S.** Table of nonzero initial concentrations in the composite model. **Table 5S.** Steady state analysis of the Bentele’s model and the composite model. Species. **Table 6S.** Steady state analysis of the Neumann’s model and the composite model, **Table 7S.** Calculation of the mean sensitivity for the investigated apoptosis models. **Table 8S.** Analysis of predictions regarding apoptosis in HeLa cells as formulated by Neumann et al. **Table 9S.** Predictions of the models for SKW 6.4 cells. **Table 10S.** Parametric constraints of the models by Bentele et al. and Neumann et al.

Click here for file

The model comprises 43 species (including _{
aa
}) and 45 kinetic parameters, which the authors estimated based on the experimental data obtained by Western blot analysis for the human cell line SKW 6.4. Cells were stimulated by 5 μg/ml and 200 ng/ml of anti-CD95 (fast and reduced activation scenarios, respectively) and dynamics of several proteins (Bid, tBid, PARP, cPARP, procaspases-2, -3, -7, -8, -9, cleaved product of procaspase-8 p43/p41, and caspases-8) were measured.

We could not reproduce the dynamics of the original model using the parameter values provided by the authors. In particular, there was too rapid consumption of procaspases-2, -3, -7, -8 in the case of 5 μg/ml of anti-CD95 and procaspases-2, -9 in the case of 200 ng/ml. Degradation rates of procaspase-8 and caspases-8 was insufficient for both activation scenarios. Thus, we had to make several modifications of the original model to obtain the same dynamics as described in the original paper. Namely, we multiplied the rate constants of all the bimolecular reactions by the value 5∙ 10^{-5} and specified _{
d
} equal to 3.56 min^{-1} and 0.62 min^{-1} (instead of 0.891 min^{-1} and 0.184 min^{-1}) for fast and reduced activation scenarios, respectively.

The model by Neumann

The original model by Neumann

**The original model by Neumann ****and results of the model reduction. A.** The original model decomposed into modules according to three steps of the CD95 signaling pathways: activation of caspase-8, pro-apoptotic pathway resulting in caspase-3 activation and anti-apoptotic pathway regulated by NF-κB. The species retained during the model reduction are colored; reactions are represented by solid lines. **B.** The modular view of the model. Activation of caspase-8 and p43-FLIP (product of cFLIPL cleavage) occurs at the DISC complex and triggers simultaneous processes of cell death and survival.

Reduction of the CD95-signaling model

We started the Bentele’s model reduction by excluding the direct dependence of the virtual species _{
aa
} on caspases-6 and -7. For this purpose, we approximated the amount of caspases-7 by the linear function of caspase-3 concentration:

where

and simplified it:

After that we decomposed the model into modules (Figure

Below we provide detailed description of the reduction process of all modules. The process is based on the flow chart represented in the Figure

Firstly, we eliminated slow reactions br12 and br13 according to the method MR1. Next, we applied the quasi-steady-state analysis (MR2) to the module and removed quasi-stationary intermediates CD95R:CD95L, DISC:pro8, DISC:pro8_{2}, DISC:p43/p41 and DISC:cFLIPL. Thus, in particular, we got two main reactions instead of br1-br6: fast formation of the DISC complex (Table

Further, we noticed that the consumptions of cFLIPL and cFLIPS satisfied the same kinetic laws. Therefore, these species could be lumped (MR3) resulting in the reaction br4*:

cFLIP + 2∙ DISC + pro8 → cFLIP:DISC:pro8,

where cFLIP indicated two isoforms cFLIPL and cFLIPS so that _{
cFLIP
} = _{
cFLIPL
} = _{
cFLIPS
}.

We approximated the concentration of caspase-6 by the linear function _{
casp3},

pro8 –casp3 → casp8.

Since the reaction br8 followed the Michaelis-Menten kinetics with the constants _{68} and _{68}, then the reaction rate of br3* was provided by the kinetic law

where _{38} = 0.145 · _{68} and _{38} = _{68}. Analysis of this kinetic law for low level of CD95L ensured that _{
pro8} ≫ _{38}. In addition, in the case of the fast activation scenario the value of _{
br3*
} was much lower than the rate of caspase-8 activation mediated by CD95. Thus, we established _{
br3*
} = _{
38
}·_{
casp3
} without significant changes in the results of the model simulation (MR1, MR5).

Complementing the above reduction steps, we removed the elements that were unnecessary to fit the experimental dynamics provided by Bentele

where _{0} and _{0} denote initial concentrations of the ligand and receptor, respectively. Solving the differential equation

we obtained the analytical function for the receptor concentration:

Accordingly to Bentele

Here _{
release
}(

where

Then, we noticed that the experimental measurements of procaspase-9 concentration were presented by Bentele and his colleagues only for 200 ng/ml of anti-CD95. In this case, degradation of procaspase-9 was insignificant. Thus, neglecting it and considering the kinetic law _{
Apop
}·_{
pro9
}·_{
CytC:Apaf-1
} of the second reaction in (15), we obtained the differential equation of procaspase-9 dynamics:

Solving it, we found:

Finally, analyzing reactions of Bid activation (br21-br23), we removed the slower reaction mediated by caspase-2. The similar reaction involving caspase-8 followed the Michaelis-Menten kinetics with the constant _{28Bid
} ≫ _{
Bid
}. Therefore, we redefined the kinetics of this reaction based on the law of mass action (MR5):

_{
aa
}) and 11 reactions (excepting reactions of degradation). Six reactions (activation of caspases-3 and -7 by caspases-8 and -9, and cleavage/inactivation of PARP) were based on the Michaelis-Menten kinetics (br27-br32), four reactions reproduced the reversible inhibition of caspases-3 and -7 based on the law of mass action (br33, br34, br36, br37), and one corresponded to the production of _{
aa
} (br38), simplification of which was discussed above.

We deleted slow reactions of casp3:IAP and casp7:IAP complexes dissociation and ignored slow degradation of procaspase-7 and IAP. Then analyzing the cleavage of procaspases-3 and -7, we eliminated reactions triggered by caspase-9, which were slower in comparison to the similar reactions induced by caspase-8. For the same reason we ignored the reaction of PARP inactivation by caspase-3 (MR1).

Considering the remaining reactions with the Michaelis-Menten kinetics (br29, br32), we used the following inequalities _{78} ≫ _{
pro7} and _{376act
} ≫ _{
PARP
}, which, in combination with MR5 and (12), led to the following kinetic laws:

where _{
3act
} = 0.18 · _{
7act
}

Thus, the CD95-signaling model consisting of 43 species (including one virtual species), 80 reactions and 45 kinetic parameters was approximated by a model with 18 species, 26 reactions and 25 parameters, except the constant concentration of Cytochrome

Results of the composite model approximation to the experimental data by Bentele

**Results of the composite model approximation to the experimental data by Bentele **** et al.** The experimental data (dots) were obtained by Bentele

Species without experimental evidence (CD95L, DISC, cFLIP, DISC:pro8:cFLIP, IAP, caspase-3 and virtual species _{
aa
}) and reactions br1*-br13* directly affect the simulated dynamics of the experimentally measured concentrations (Table

**Species**

**Reactions**

**The role in the apoptosis process**

“Experimental data” in this table refers to the data obtained by Bentele

CD95L

br1*

DISC

br2*

caspase-3

br3*

br8*

br9*

br11*

cFLIP

br4*

DISC:cFLIP:pro8

br5*

br6*, br7*, br10*

IAP

br12*

_{
aa
}

br13*

The virtual variable _{
aa
} specifies the process of degradation.

Reduction of the CD95-mediated and NF-κB signaling model

As was mentioned above, Neumann

Based on the method MR6, we also took into account the inequality _{
IKK
} ≫ _{
p43−FLIP
} and simplified the kinetic law

of the reaction nr19 to the form

Therefore, we reduced the number of model species from 23 to 20, the number of reactions from 23 to 18 and the number of kinetic parameters from 17 to 15 (besides the constant concentration of IKK). These modifications did not change the fit to the experimental data provided by Neumann

Comparison of the composite model simulation results with the experimental data by Neumann

**Comparison of the composite model simulation results with the experimental data by Neumann **** et al.** The simulated concentrations of the original and reduced models (solid lines), as well as the simulated values of the composite model (dashed lines) were obtained by Neumann

Model composition

Analysis of the reduced models based on the Akaike criterion (7) confirmed that they had lower complexity than the original models (Table

**Bentele ****
et al.
**

**Neumann ****
et al.
**

**Composite model**

**Original model**

**Reduced model**

**Original model**

**Reduced model**

**
AIC
**

182.46

138.83

201.58

198.17

318.10

**Mean ****
AIC
**

1.50

1.14

0.96

0.94

0.96

We took into account that the considered experimental data were obtained with SKW 6.4

(A) initial species concentrations could vary for different cell lines;

(B) kinetic parameters of all reactions (besides the degradation rate modeled as the function _{
d
} · _{
aa
}
^{2}) have the same values for both cell lines, whereas the value of _{
d
} could be regulated by various entities and, moreover, is dependent on the initial concentration of anti-CD95

(C) type I and type II cells conform to different reaction chains of caspases-3 activation

We constructed the composite model of CD95 and NF-κB signaling in the following way. Firstly, according to the assumption (C), we replaced reaction nr15 in the Neumann’s model by a chain of three reactions:

The first reaction of this chain had the same kinetic law as br7*. The other reactions represented br6* and nr15 with kinetics modified using the law of mass action. Such modification of br6* was necessary to get the slow increase of caspase-3 concentration experimentally described by Neumann _{
casp9} for _{
casp8} .

Secondly, we supplemented the Neumann’s model with reaction of caspase-3 inhibition by IAP and evaluated parameters in this reaction together with parameters in (17) using the corresponding data by Neumann

Thirdly, we analyzed the changes of the initial concentrations and parameters of the derived model required to reproduce the Bentele’s experimental data fixing levels of procaspases-3, -8, cleaved product p43/p41, and caspase-8 (Table

**Step**

**Changed parameters and concentrations**

**Initial values, modifications**

**Reason**

**Ranges and initial values of fitting**

1

CD95R:FADD

91.266 nM, increase

[10^{2}, 10^{3}], 442.821, Bentele

2

procaspase-8

64.477 nM, increase

[10^{2}, 10^{3}], 442.821, Bentele

3

FLIPS

5.084 nM, increase

[10^{1}, 10^{2}], 65.021, Bentele

4

k10

0.121 nM^{-1}min^{-1}, decrease by an order of magnitude

0.012

5

procaspase-3

1.443 nM, refitting

[10^{0}, 10^{1}], 1.443, Neumann

6

Bid, procaspase-9

5.003 nM, 2.909 nM, replacement by Bentele’s values

231.760 (Bid), 245.101 (pro-9), Bentele

Finally, we combined the modules of caspase-8 and NF-κB activation with the modules of caspase-9 activation and PARP inactivation (Figure

Creation of the composite model

**Creation of the composite model.** The modules of caspase-8 and NF-κB activation were taken from the reduced Neumann’s model and combined with the modules of caspase-9 activation and PARP inactivation isolated from the reduced Bentele’s model. In addition, two last modules were supplemented by reactions from the modified module of caspases-3 activation

The resulting composite model (Figure

The composite model of the CD95- and NF-κB-signaling

**The composite model of the CD95- and NF-**κ**B-signaling. A.** The model integrating pro- and anti-apoptotic machinery described by the models by Neumann **B.** The modular view of the model. Activation of caspase-8 triggered by CD95 leads to NF-κB activation and cell survival, on the one hand, and PARP cleavage resulted in the cell death, on the other hand.

Calculation of the mean

The reduced models are equivalent to the precursor models with respect to their biological predictions

Analyzing the models constructed by Neumann

The second group of predictions characterized behavior of the models after changing concentrations of some species, such as CD95L, CD95R, procaspase-8, and inhibitors cFLIPL, cFLIPS, and IAP. Analyzing the predictions of the models by Neumann

**№**

**The composite model behavior**

**Predictions by Neumann ****
et al.
**

All the predictions marked with an asterisk were experimentally tested by Neumann

1*

The concentration of anti-CD95 required for the apoptosis induction (the apoptotic threshold), is within the range of 30–100 ng/ml. This range remains the same for CD95 decreased by about 12-fold. The simulation time, which we used to reproduce this prediction, was 60 hours.

2*

The decreased receptor number results in impairment of both CD95- and NF-κB-signaling pathways. To test this prediction, Neumann

3*

Along with increasing the concentration of anti-CD95 from 500 ng/ml to 1500 ng/ml, p43/p41 peaks earlier, while there is almost no difference for p43-FLIP.

4

Increased concentrations of cFLIPS inhibit both apoptotic and NF-κB pathways, although p43-FLIP generation is inhibited at a lower threshold than p43/p41 generation.

5*

Increasing the concentration of cFLIPL leads to a steep increase in p43-FLIP generation until it reaches a maximum, after which the curve drops. Lowered levels of cFLIPL result in very little p43-FLIP but almost unchanged levels of p43/p41.

At very high concentrations of cFLIPL no p43-FLIP is generated. This drop-off was not observed experimentally by the authors.

6*

Only an intermediate level of cFLIPL promotes NF-kB activation. Decreased levels of procaspase-8 lead to a significantly lower amount of p43-FLIP and, subsequently, NF-κB. The figures show of logarithmic dependence of the maximal NF-κB concentration on the initial values of procaspase-8 and cFLIPL

7*

High cFLIPL or low procaspase-8 concentrations cause suppression of apoptosis. The figures show the same dependence as considered in the previous prediction, but with caspases-3 instead of NF-κB.

**№**

**The models behavior**

**Experimental observations by Bentele ****
et al.
**

The simulation time was 2880 min (2 days) in all predictions. The apoptotic threshold in the prediction 1 is the concentration of anti-CD95 after which cPARP amount exceeds 10% of the initial PARP level.

1

**Experiments by Bentele ****
et al.
**

− for 1 ng/ml of anti-CD95, PARP cleavage was not observed;

− the measured death rate for 10 ng/ml of anti-CD95 was 20-30%.

**Original model (red):**

− the apoptotic threshold is 1.9 ng/ml;

− cPARP concentration rises dramatically within an extremely narrow interval of anti-CD95 levels overcoming the apoptotic threshold.

**Composite model (blue):**

− the apoptotic threshold is 3.5 ng/ml;

− cPARP concentration rises in a smooth manner along with the increase of anti-CD95 level.

2

**Experiments by Bentele ****
et al.
**

**Original (red) and composite (blue) models:**

− the apoptotic threshold is highly sensitive to the concentration of cFLIP;

− decreasing the initial concentration of cFLIP by more than 51% and 49% for the original and composite models, respectively, leads to cell death upon stimulation by 1 ng/ml of anti-CD95.

3

**Experiments by Bentele ****
et al.
**

**Original (red) and composite (blue) models:** anti-CD95 concentrations which are slightly above the apoptotic threshold result in caspase-8 activation after a delay of many hours.

The figure shows peak times of caspase-8 concentration exceeding 0.1% of the initial procaspase-8 level.

4

**Original model:**

− low concentrations of IAP (less than 1 nM) result in complete cell death;

− high concentrations of IAP prevent a significant increase of caspase-3 even for high concentrations of the ligand.

**Composite model:**

− low concentrations of IAP (less than 1 nM) block apoptosis for CD95L less than 0.3 nM;

− high concentrations of CD95L lead to cell death.

The figures show logarithmic dependence of the maximal caspases-3 concentration on initial values of IAP and CD95L.

Model composition modifies some properties of the precursor models and leads to new predictions

Considering the predictive ability of the composite model, we found that in the case of HeLa cells it was completely preserved (Table

In the case of SKW 6.4 cell line, we estimated the value of degradation rate parameter _{
d
} based on the experimental observations by Bentele _{
d
} depending on whether the initial concentration of CD95L was within the range 1–100 ng/ml, within the range 100–1000 ng/ml or greater than 1000 ng/ml (Additional file

Analyzing the model predictions, we detected differences in the behavior of the Bentele’s and the composite models (Table

The next prediction, which we observed for SKW 6.4 cells, considers a delay of caspase-8 activation caused by low concentrations of anti-CD95 (Table

Limitations of the composite model

We analyzed a set of parametric constraints which allowed us to reduce the precursor models (Additional file

Discussion

In this paper, we considered two approaches to development of mathematical models of cell fate decisions. The first concerns the methodology of model reduction and involves approximation of one model by another one of lower dimension without affecting dynamics of experimentally measured species. The second implies composition of the models and aims at reproducing experimental dynamics of all precursor models.

There are many advantages of working with low-dimensional models

Model composition aiming at getting a single model from several ones is useful because in such a case a computational biologist is able to investigate the composite model behavior under different conditions that cannot be performed in the precursor models separately. For example, our model allows studying the role of p43-FLIP or IAP in the type I SKW 6.4 or type II HeLa cells, respectively, that might become a task for the future work. In other words, we constructed the model that describes pro- and anti-apoptotic signal transduction in different cell types with reasonable accuracy instead of a couple of different models. Whenever necessary, some reaction chains and parameters can be switched off giving opportunity to simulate a certain type of cells. In addition, the composite model covers experimental data obtained from all precursor models, each of which separately satisfies its own data only.

Model composition sometimes causes modification of some properties of the initial models, resulting in new testable predictions. In our case, such predictions were related to SKW 6.4 cells, and some simulation results were different from the corresponding results of the original Bentele’s model. Nonetheless, the composite model behavior remained in good agreement with experimental data used in modeling by Bentele

One of the questions that can be asked concerns the need of reducing models in order to combine them. An attempt to construct the composite modular model of apoptosis was made in

choice of elements in the overlapping parts of different models;

incomparable sets of parameters;

the lack of experimental data, as well as inability to use the data obtained for different cell lines or in different ways (e.g. single-cell or cell culture measurements);

inability to make accurate predictions.

Model reduction allows solving some of these problems. In particular, reducing the number of model elements, we reduce the overlapping parts as well. This may be essential for direct combining of models. In addition, the complexity of the reduced model become comparable with the complexity of available experimental data. Therefore, the risk of model overfitting is decreased.

In the case when two models are directly related (for example, one model was emerged from the other), their composition may be significantly easier in comparison to composition of quite different models. In our work, the models by Bentele

Another quite useful principle that we used in modeling was modular structure of the developed model. This principle provides flexibility for future extensions. Thus, we are planning to extend the composite model and supply it with modules and data from studies of TRAIL

Characterising the models used in our work, we can say that we saved our time working with the model by Neumann

Finally, we believe that the composite model itself may be useful for further investigation of apoptosis.

Conclusions

Mathematical modeling provides a powerful tool for studying the properties of biological processes. Methods of model reduction allowed us to take a first step towards validation of the modular model of apoptosis

The models by Bentele

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

EK performed reduction, integration and analysis of the apoptosis models. AZ coordinated the model reduction part of the study. RS and FA coordinated the study of the apoptosis machinery. FA also led implementation of the analysis toolbox in the BioUML software. All authors have participated in writing and editing the manuscript, as well as read and approved the final manuscript.

Acknowledgements

This work was supported by FP6 grant 037590 "Net2Drug", FP7 grant 090107 "LipidomicNet", EU FP7 (project APO-SYS) and Agence Nationale de la Recherche (project ANR-08-SYSC-003 CALAMAR). AZ is a member of the team “Systems Biology of Cancer,” Equipe labellisée par la Ligue Nationale Contre le Cancer.