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  <cit> .
several aspects of the same process can be modeled using different levels of abstraction involving different reaction chains, chemical kinetics, and incomparable sets of parameters.
such models are difficult to merge.
meanwhile, merging is an important approach for creation of complex models.
thus, development of efficient methods and software, allowing us to combine models, is the object of intense study in systems biology.
in our work, we focused on the fact that, generally, complexity of models is not comparable to the volume of experimental data used to adjust their parameters.
due to this fact, we turn to the methods of model reduction allowing us to minimize model’s complexity without affecting the model simulation dynamics.
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  and has already found multiple applications in systems biology, including discrete modeling  <cit>  and modeling of metabolic pathways  <cit> .
the principles of this technique have been applied in computational biology  <cit>  and implemented as a part of widely used pathway simulators such as biouml  <cit> , biocham  <cit> , copasi  <cit> , and ginsim  <cit> .
model reduction led to new insights in mechanisms of translation regulation by micrornas  <cit>  and was applied for analysis of such signaling pathways as jak-stat  <cit> , nf-κb  <cit> , and egfr  <cit> .
in our investigation, we used the principles of model reduction to construct reasonably accurate minimal size approximations of two different models describing the cd <dig> signaling pathways  <cit> .
the first model explores pro-apoptotic properties of cd <dig> after stimulation by its natural ligand cd95l or by agonistic antibodies anti-cd <dig> implying formation of the death-inducing signaling complex   <cit> .
disc consists of oligomerized cd <dig>  death domain-containing adaptor fas-associated molecule , procaspases- <dig> and − <dig>  and two isoforms of cellular flip  protein .
caspase- <dig> leads to activation of effector caspase- <dig> directly  or via stimulation of cytochrome c release from the mitochondria   <cit> .
the latter step requires formation of the apoptosome complex and activation of caspase- <dig>  once activated, caspase- <dig> cleaves poly polymerase , thereby making the apoptotic process irreversible.
the second model describes the state when cd <dig> not only activates the pro-apoptotic pathway, but also induces transcription factor nf-κb that is an important regulator of cell survival functioning  <cit> .
this is possible due to cflipl cleavage in the disc complex.
the cleaved p43-flip directly interacts with the ikk complex and activates it.
the activated ikk performs phosphorylation of the iκb protein and thereby frees nf-κb.
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  <cit>  implemented in the biouml software, we reduced the models so that they still satisfy these data.
this allowed us to simplify the overlapping components of the models and find a way for their composition.
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 s =  associated with a set of variables c = , …, cm) depending on time t ∈  , t ∈  r+ and representing their concentrations, and a set of biochemical reactions with rates v = , …, vn) depending on a set of kinetic constants k. reaction rates are modeled by mass-action or michaelis-menten kinetics.
a system of ordinary differential equations  representing a linear combination of reaction rates is used to describe the model behavior over time:
dctdt=n·vc,k,t,c0=c <dig>
here n is a stoichiometric matrix of n by m. we say that css is a steady state of the system  if
n·vcss,k,t= <dig> limt→∞cit=ciss.
model reduction implies transformation of the ode system to another one with smaller number of equations without affecting dynamics of variables c <dig>  …, cs, which is fixed by a set of experimental points
ciexp at given times tij, j =  <dig>  …, ri, where ri is number of such points for the concentration ci, i =  <dig>  …,s.
to check dynamics preservation in the course of model reduction, we consider the function of deviations defined as a normalized sum of squared differences  <cit> :
fdistc <dig> k=∑i=1s∑j=1riωminωi·citij−ciexptij <dig>
where ωmin=miniωi, weights ωi=ri−1·∑jciexptij <dig> are calculated as mean squared values of experimental concentrations for each variable and the normalizing factor ω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  between model parameters.
when these relations satisfy certain rules, we can approximate the detailed model by a simpler one.
if the parameters are  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 .
we define the minimal complexity model as a reduced model with the minimal number of elements  so that the deviation function value  calculated for the reduced model is different from the original model no more than 20%.
such threshold is explained by the fact that in this work we considered experimental data obtained by bentele et al.
<cit>  and neumann et al.
<cit>  using western blot technology with the standard deviation 15-20%.
reduction of mathematical model complexity is achievable by different methods  <cit> .
description of the methods directly used in our work is provided below.
removal of slow reactions.
we say that reaction r <dig> is much slower than reaction r <dig> vr1t≪vr2t if maxtvr1t<k·maxtvr2t, where k=10- <dig>  when such reactions do not affect experimental dynamics of species, we neglect them.
note, that in some cases, it is sufficient to consider k = 10− <dig> or even k =  <dig>
quasi-steady-state approximation <cit>  assumes that all variables c of the model are split in two groups: basic variables cs and quasi-stationary variables cf for which we could write:
dcstdt=ns·vsc,k,t,dcftdt=1ε·nf·vfc,k,t,
where ɛ is small parameter.
we consider the approximation ɛ →  <dig> resulting in the equation nf · vf =  <dig>  thus, the chain of reactions s → i → p with a quasi-stationary species i can be reduced to the form s → p.
lumping analysis refers to reducing the number of model variables by grouping some species.
in particular, this can be illustrated by reactions
r1:a+s1→p <dig> r2:a+s2→p <dig>
with kinetic rates vr1t=k1·cat·cs1t and vr2t=k2·cat·cs2t.
if k1 = k <dig> and cs1t=cs2t, then the system  can be replaced by a single reaction
2a+s→p,
where cst=cs1t=cs2t, cpt=cp1t=cp2t and k = k1 = k <dig>  note, that in this example lumped variables cs1t, cs2t and cp1t, cp2t are linearly dependent: however, this is not required in general case.
for conditions on lumpability in monomolecular reaction networks, see  <cit> .
removal of approximately linearly dependent variables.
if variables a and b  are approximately linearly dependent:
at≈k·bt,t∈ <dig> t,k=const,
then we can replace one of them with another in all kinetic laws of the model.
simplification of the michaelis-menten kinetics.
consider a reaction r of the form s –e → p, where an enzyme e converts a substrate s into product p. reaction rate of such reaction is frequently defined by the michaelis-menten formula
vrt=k·cst·cetkm+cst
where the enzyme concentration is a dynamic variable allowing to use the same kinetics in different regions of the phase space.
if km ≫ cs, then we can reduce this formula to the form: vrt=k·cst·cetkm.
on the other hand, if cs ≫ km, then  can be replaced by the equation vr = k · ce.
simplification of equations based on the law of mass action when one reactant dominates others.
consider a reaction r of the form s1 + s2 → p with the kinetic law
vrt=k·cs1t·cs2t.
when cs1t≫cs2t for t ∈  , the formula  can be replaced by the linear equation vrt=k·cs10·cs2t.
the time t of validity of such pseudo-lineary approximation is defined as a period during which the relative change of cs1t does not exceed some ε, i.e.
k∫0tcs2tdt<ε.
one way to perform the model reduction is to apply the foregoing methods in the numerical order .
note that the methods mr <dig> and mr <dig> are of the same type, so it does not matter which one to use first.
this way assumes that we calculate the value of the distance function  at each reduction step to test whether approximation of the experimental data is still within allowable limit or not.
model analysis and comparison
model comparison using akaike information criterion .
this criterion  <cit>  defines the relative complexity of a model based on the goodness of experimental data fit and the number of model parameters |k| :
aic=χ2+2·k.
the function χ <dig> is defined by the formula  with weights 1/
σij <dig>   <cit> , where mean deviations σij of experimental values ciexp are calculated using smoothing spline  <cit> :
σij2=14∑k=−22citi,j+k−ciexpti,j+k <dig>
for j + k <  <dig> and j + k > ri we assumed ci = ciexp =  <dig>
models can be compared using the akaike criterion, if they approximate the same experimental data.
when the aic value of one model is less than of others, we say that this model is simpler in terms of this criterion.
model comparison with the mean aic.
when we want to compare models approximating different sets of experimental data, we could calculate the mean aic coefficients by the formula
aicmean=aicnexp
where nexp determines the number of experimental points.
steady-state analysis finds values of species concentrations according to the rules .
it is desirable to reduce a model so that steady-state concentrations of all experimentally measured chemical species have not changed significantly.
sensitivity analysis reveals steady-state concentrations response to parameter perturbations.
the local sensitivities ciss, i =  <dig>  …, m, for all parameters pj, j =  <dig>  …, |k|, are calculated via finite difference approximations:
sij=∂ciss∂pj=cisspj+Δpj−cisspjΔpj
it is useful to compute the mean log sensitivity for all fitted species:
smean=∑i,jlnsijnp·s,sij≠ <dig> i= <dig> …,s,j= <dig> …,np,
where np is the number of all parameters |k| or the number of parameters retained after the model reduction.
modular modeling
considering the mathematical models of cd <dig> signaling pathways  <cit> , we decomposed them into functional modules.
this step allowed us to identify overlapping components of the models and simplify their analysis.
we defined a module as a submodel  with input, output and contact ports.
the first two types of ports characterized variables calculated in one module and passed to another through a directed connection.
the contact ports declared common variables of modules via undirected connections .
RESULTS
preliminary analysis, problems and inconsistencies in the precursor models
the mathematical model by bentele et al.
consists of  <dig> species including a virtual variable xaa intended for quantification of “apoptotic activity” caused by active caspases- <dig>  - <dig> and - <dig> and calculated by the formula
dxaadt=1−xaa·k36act·ccasp3+k36act·ccasp6+k7act·ccasp7km367act+1−xaa−kdd·xaa,
where k36act, k7act, km367act ∈ k. the variable xaa specifies a process of degradation in the model introduced as an exponential decay function fdegr.
this function is defined by the formula kd · xaa <dig> + kds with kd,kds ∈ k for all active caspases and complexes containing their cleaved products, and by the formula kd · xaa <dig> for all other molecules besides cparp , where fdegr is constant.
all species in the model are degraded with the exception of cytochrome c and smac stored in the mitochondria, whose concentrations are constant.
in total, the model contains  <dig> reactions  including  <dig> reactions based on mass action kinetics,  <dig> reactions taken with kinetics of michaelis-menten,  <dig> reactions of degradation, the reaction of xaa production specified above and two reactions of cytochrome c and smac release from the mitochondria modeled using a discrete event.
the latter implies complete cytochrome c and smac release within  <dig> minutes as soon as tbid reaches a certain level in comparison to bcl-2/bcl-xl.
since the authors did not provide exact form of the release function, we proposed the sigmoid function which has the expected behavior:
freleaset=1−11+expkcontr·−t+ttrigger+ <dig> ·trelease
where ttrigger is the time when tbid concentration reaches a value of bcl-2/bcl-xl, trelease is the start time of release and kcontr is the contraction coefficient.
prototypes of parameter notations in the original bentele’s model are listed in additional file 1: table 1s.
abbreviations used in the table represent Мa – mass action kinetics and m-m – michaelis-menten kinetics.
the column with rates includes rate orders for the fast/reduced activation scenarios.
the model comprises  <dig> species  and  <dig> kinetic parameters, which the authors estimated based on the experimental data obtained by western blot analysis for the human cell line skw  <dig> .
cells were stimulated by  <dig> μg/ml and  <dig> ng/ml of anti-cd <dig>  and dynamics of several proteins  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- <dig>  - <dig>  - <dig>  - <dig> in the case of  <dig> μg/ml of anti-cd <dig> and procaspases- <dig>  - <dig> in the case of  <dig> ng/ml.
degradation rates of procaspase- <dig> and caspases- <dig> 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- <dig> and specified kd equal to  <dig>  min- <dig> and  <dig>  min- <dig>  for fast and reduced activation scenarios, respectively.
the model by neumann et al.
includes  <dig> species,  <dig> reactions constructed according to the law of mass action, and  <dig> kinetic parameters .
it reproduces experimental measurements of  <dig> entities  obtained using western blot technology for hela cells stably overexpressing cd95–gfp and treated with three different concentrations of agonistic anti-cd <dig> antibodies .
this is the reduced model suggested by the authors for description of the cd95-mediated pathways of cell death and nf-κb activation.
model reduction was based on intuition about what is important in the biochemical mechanisms.
it was performed as long as the parameter estimation procedure  was able to provide a good fit of the experimental data.
particularly, the details of some complex formations were omitted as well as some non-measured species were removed.
the authors submitted the model to the biomodels database  <cit> .
thus, we had no problems reproducing it.
reduction of the cd95-signaling model
we started the bentele’s model reduction by excluding the direct dependence of the virtual species xaa on caspases- <dig> and - <dig>  for this purpose, we approximated the amount of caspases- <dig> by the linear function of caspase- <dig> concentration:
ccasp7≈a·ccasp <dig>
where a =  <dig>  for both activation scenarios.
in addition, we considered inequalities holding for the parameters in the formula :
k7act·ccasp7≫k36act·ccasp3>k36act·ccasp <dig>
km367act≫1−xaa,
and simplified it:
dxaadt=k36act+ <dig> ·k7actkm367act·1−xaa·ccasp3−kdd·xaa.
after that we decomposed the model into modules  corresponding to three biological steps: activation of caspase- <dig>  cytochrome c-induced activation of caspase- <dig> and parp cleavage down-regulated by these caspases.
below we provide detailed description of the reduction process of all modules.
the process is based on the flow chart represented in the figure  <dig>  if it is enough for a reader to have a general idea of model reduction without going into full details, we suggest to go directly to the last two paragraphs of this section, where you find the summary of the model reduction procedure.
caspase- <dig> activation module includes  <dig> species and  <dig> reactions , which could be divided into three groups: cleavage of procaspase- <dig> at the disc complex, its activation triggered by caspase- <dig>  and inhibition of complexes containing disc by cflipl and cflips .
firstly, we eliminated slow reactions br <dig> and br <dig> according to the method mr <dig>  next, we applied the quasi-steady-state analysis  to the module and removed quasi-stationary intermediates cd95r:cd95l, disc:pro <dig>  disc:pro <dig>  disc:p43/p <dig> and disc:cflipl.
thus, in particular, we got two main reactions instead of br1-br6: fast formation of the disc complex  and slower activation of caspase- <dig> in this complex .
further, we noticed that the consumptions of cflipl and cflips satisfied the same kinetic laws.
therefore, these species could be lumped  resulting in the reaction br4*:
cflip + 2∙ disc + pro <dig> → cflip:disc:pro <dig>
where cflip indicated two isoforms cflipl and cflips so that ccflip = ccflipl = ccflips.
we approximated the concentration of caspase- <dig> by the linear function b · ccasp <dig>  b =  <dig>  , and merged reactions br <dig> and br <dig> of caspases- <dig> activation induced by caspases- <dig> and mediated by caspases- <dig> into one reaction br3*:
pro <dig> –casp <dig> → casp <dig>
since the reaction br <dig> followed the michaelis-menten kinetics with the constants k <dig> and km <dig>  then the reaction rate of br3* was provided by the kinetic law
vbr3*=k38·cpro8·ccasp3/km38+cpro <dig>
where k <dig> =  <dig>  · k <dig> and km <dig> = km <dig>  analysis of this kinetic law for low level of cd95l ensured that cpro <dig> ≫ km <dig>  in addition, in the case of the fast activation scenario the value of vbr3* was much lower than the rate of caspase- <dig> activation mediated by cd <dig>  thus, we established vbr3* = k38·ccasp <dig> without significant changes in the results of the model simulation .
complementing the above reduction steps, we removed the elements that were unnecessary to fit the experimental dynamics provided by bentele et al.
these elements were degradation reactions of cd95l, cd95r and flip, and the blocked disc complex .
as a result, we introduced the conservation law:
ccd95l–ccd95r=l0–r0=const,
where l <dig> and r <dig> denote initial concentrations of the ligand and receptor, respectively.
solving the differential equation
dccd95rdt=−klr·ccd95r·ccd95l,
we obtained the analytical function for the receptor concentration:
ccd95r=r0−l01−l0/r0·exp−klr·r0−l0·t
the module of caspase- <dig> activation consisted of  <dig> species and  <dig> reactions  including activation of caspase- <dig> triggered by cytochrome c, cleavage/activation of bid by caspases- <dig> and - <dig>  iap inhibition of caspase- <dig> activity and binding of inhibitors by smac .
based on the method mr <dig>  we eliminated the last two reversible reactions, as well as the slow reaction br <dig> of the apoptosome complex dissociation.
further, since this complex was a quasi-stationary intermediate and the rate of br <dig> was linearly dependent on a difference of br <dig> and br <dig> rates, we reduced the chain of procaspase- <dig> activation , using the methods mr <dig> and mr <dig>  to the form:
cytcstored→cytc:apaf− <dig> pro9–cytc:apaf−1→casp <dig>
accordingly to bentele et al., we have the following rule of cytochrome c release from mitochondria:
ccytcstored=ccytcstored0·freleaset,
dccytcreleaseddt=−dccytcstoreddt
here frelease satisfies the formula .
taking into account  and ignoring degradation of the complex cyt c:apaf- <dig>  we got
dccytc:apaf−1dt=−c·dccytcstoreddt
where c =  <dig>  was the coefficient of linearity.
then, we noticed that the experimental measurements of procaspase- <dig> concentration were presented by bentele and his colleagues only for  <dig> ng/ml of anti-cd <dig>  in this case, degradation of procaspase- <dig> was insignificant.
thus, neglecting it and considering the kinetic law kapop·cpro9·ccytc:apaf- <dig> of the second reaction in , we obtained the differential equation of procaspase- <dig> dynamics:
dcpro9dt=− <dig> ·1−freleaset·ccytcstored0·kapop·cpro <dig>
solving it, we found:
cpro9=cpro90·freleaset·expkcontr·t−1+1expkcontr·t <dig> ·ccytcstored0·kapopkcontr.
finally, analyzing reactions of bid activation , we removed the slower reaction mediated by caspase- <dig>  the similar reaction involving caspase- <dig> followed the michaelis-menten kinetics with the constant km28bid ≫ cbid.
therefore, we redefined the kinetics of this reaction based on the law of mass action :
vbr7*=k8bid/km8bid·ccasp8·cbid.
the module of parp inactivation contained  <dig> species  and  <dig> reactions .
six reactions  were based on the michaelis-menten kinetics , four reactions reproduced the reversible inhibition of caspases- <dig> and - <dig> based on the law of mass action , and one corresponded to the production of xaa , simplification of which was discussed above.
we deleted slow reactions of casp3:iap and casp7:iap complexes dissociation and ignored slow degradation of procaspase- <dig> and iap.
then analyzing the cleavage of procaspases- <dig> and - <dig>  we eliminated reactions triggered by caspase- <dig>  which were slower in comparison to the similar reactions induced by caspase- <dig>  for the same reason we ignored the reaction of parp inactivation by caspase- <dig> .
considering the remaining reactions with the michaelis-menten kinetics , we used the following inequalities km78 ≫ cpro <dig> and km376act ≫ cparp, which, in combination with mr <dig> and , led to the following kinetic laws:
vbr29=k78km78·ccasp8·cpro <dig> vbr32=k3actkm367act·ccasp3·cparp,
where k3act =  <dig>  · k7act
thus, the cd95-signaling model consisting of  <dig> species ,  <dig> reactions and  <dig> kinetic parameters was approximated by a model with  <dig> species,  <dig> reactions and  <dig> parameters, except the constant concentration of cytochrome c .
figure  <dig> shows the comparison of simulated concentrations of these models with experimental dynamics obtained by bentele et al.
for the  <dig> species mentioned above.
since this dynamics was expressed in arbitrary units, we uniquely translated it into precise values for all proteins with non-zero initial values.
however, for such molecules as caspase- <dig>  p43/p <dig>  tbid and cparp, precise values directly depended on the concentration levels at the end or in the middle of time series.
accordingly, we computed these levels during the model reduction and recalculated experimental dots if necessary.
species without experimental evidence  and reactions br1*-br13* directly affect the simulated dynamics of the experimentally measured concentrations .
thus, further reduction of the model by removing these elements is impossible.
regarding the degradation process, we agree with bentele et al.
that experimental data cannot be matched if this process is fully ignored.
actually, if we remove any of the retained reactions of degradation, then excepted concentration dynamics will be impaired.
thus, we found the minimal complexity approximation of the original model.
“experimental data” in this table refers to the data obtained by bentele et al.
<cit> .
reduction of the cd95-mediated and nf-κb signaling model
as was mentioned above, neumann et al.
simplified their model in order to reduce the large number of free parameters.
the authors noted that further simplification of the model was not possible because the model fit significantly decreased.
however, we removed slow reactions  using the method mr <dig>  and removed the reaction of nf-κb degradation , which was not significant for reproducing the experimental data of the original model.
note that we retained the slow reactions nr <dig> and nr <dig>  the former of which regulates apoptosis inhibition in the case of high levels of cflipl and procaspase- <dig> in accordance with the precursor model prediction , and the latter of which will be required to combine this model with the bentele’s model hereafter.
based on the method mr <dig>  we also took into account the inequality cikk ≫ cp43−flip and simplified the kinetic law
vnr19=kp43−flip_ikk·cikk·cp43−flip
of the reaction nr <dig> to the form
vnr19*=kp43−flip_ikk·cikk0·cp43−flip.
therefore, we reduced the number of model species from  <dig> to  <dig>  the number of reactions from  <dig> to  <dig> and the number of kinetic parameters from  <dig> to  <dig> .
these modifications did not change the fit to the experimental data provided by neumann et al.
.
model composition
analysis of the reduced models based on the akaike criterion  confirmed that they had lower complexity than the original models .
the difference between the mean aic coefficients  of the models decreased by 60%, relative to the initial value.
in addition, the reduced neumann’s model better approximated experimental data.
therefore, we used it as the basis for the model composition.
we took into account that the considered experimental data were obtained with skw  <dig>   <cit>  and hela  <cit>  cells, which were shown to behave as type i  <cit>  and type ii  <cit>  cells, respectively.
in this regard, we assumed that:
initial species concentrations could vary for different cell lines;
kinetic parameters of all reactions  have the same values for both cell lines, whereas the value of kd could be regulated by various entities and, moreover, is dependent on the initial concentration of anti-cd <dig>  <cit> ;
type i and type ii cells conform to different reaction chains of caspases- <dig> activation  <cit> .
we constructed the composite model of cd <dig> and nf-κb signaling in the following way.
firstly, according to the assumption , we replaced reaction nr <dig> in the neumann’s model by a chain of three reactions:
bid–casp8→tbid,pro9–tbid→casp <dig> pro3–casp9→casp <dig>
the first reaction of this chain had the same kinetic law as br7*.
the other reactions represented br6* and nr <dig> with kinetics modified using the law of mass action.
such modification of br6* was necessary to get the slow increase of caspase- <dig> concentration experimentally described by neumann et al.
<cit> .
modification of nr <dig> consisted in substitution of ccasp <dig> for ccasp <dig> .
secondly, we supplemented the neumann’s model with reaction of caspase- <dig> inhibition by iap and evaluated parameters in this reaction together with parameters in  using the corresponding data by neumann et al.
and optimization tools of biouml  <cit> .
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- <dig>  - <dig>  cleaved product p43/p <dig>  and caspase- <dig> .
all the changes, except step  <dig>  agreed with assumptions  and .
to reproduce the species dynamics, we had to increase the initial concentration of procaspase- <dig> in the case of hela cells by an order of magnitude.
finally, we combined the modules of caspase- <dig> and nf-κb activation with the modules of caspase- <dig> activation and parp inactivation .
for this purpose, we modified the last two modules based on  and refitted their parameters as well as the concentrations mentioned in table  <dig>  using experimental data by bentele et al.
for better fit we also supplemented the model by degradation reactions of some species .
the resulting composite model  consists of  <dig> species ,  <dig> reactions  and  <dig> kinetic parameters, except for the constant concentration of ikk .
the model showed a good fit between the simulation results and the experimental data .
calculation of the mean aic coefficient for this model revealed that it has the same level of complexity as the reduced models .
we also noted that reduction of both models did not significantly alter the steady-state concentrations .
however, the composite model has different steady-states, which are more stable according to the mean sensitivities  in the cases of skw  <dig>  cells and hela cells stimulated with  <dig> ng/ml of anti-cd <dig> .
the reduced models are equivalent to the precursor models with respect to their biological predictions
analyzing the models constructed by neumann et al.
and bentele et al., we divided predictions made by them into two groups: qualitative and quantitative.
the first group concerned the model network and described mechanisms of the protein-protein interactions.
it is this group to which we assigned the experimentally validated neumann's predictions that processed caspases are not required for nf-κb activation, as well as that pro-apoptotic and nf-κb pathways diverge already at disc.
the corresponding network of reactions was preserved in the reduced model.
thus, the prediction remained valid.
the second group of predictions characterized behavior of the models after changing concentrations of some species, such as cd95l, cd95r, procaspase- <dig>  and inhibitors cflipl, cflips, and iap.
analyzing the predictions of the models by neumann et al.
and bentele et al.
listed in the tables  <dig> and  <dig> respectively, we concluded that the reduction of the first of them did not alter the species dynamics .
for the second reduced model, we detected some minor changes in the predictions.
however, in general, the model behavior remained in a good agreement with the bentele's experiments .
all the predictions marked with an asterisk were experimentally tested by neumann et al.
and confirmed, unless otherwise noted.
the simulation time in predictions 3– <dig> was  <dig> min.
the concentration of anti-cd <dig> considered by the authors in predictions 4– <dig> was  <dig> ng/ml.
the simulation time was  <dig> min  in all predictions.
the apoptotic threshold in the prediction  <dig> is the concentration of anti-cd <dig> after which cparp amount exceeds 10% of the initial parp level.
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 .
in the case of skw  <dig>  cell line, we estimated the value of degradation rate parameter kd based on the experimental observations by bentele et al.
for 1– <dig> ng/ml of anti-cd <dig> .
when performing simulation experiments of the model for this cell line, we used three different values of kd depending on whether the initial concentration of cd95l was within the range 1– <dig> ng/ml, within the range 100– <dig> ng/ml or greater than  <dig> ng/ml .
analyzing the model predictions, we detected differences in the behavior of the bentele’s and the composite models .
the first of them describes a threshold mechanism for cd95-induced apoptosis implying that this process is completely stopped when cd95l concentration is below the critical value.
however, the model predicts a sharp increase of cparp immediately when the level of cd95l exceeds the threshold .
in contrast, the composite model shows a gradual increase of cparp with increase of the ligand level.
to determine the apoptotic threshold in this case, we found concentration of cd95l for which cparp ratio is equal to 10% of the initial parp amount  <cit> .
as in the original model, this concentration was in the range of 1– <dig> ng/ml.
additionally, the models revealed sensitivity of the threshold to the concentration of cflip .
both of them predicted the cell death phenotype upon stimulation by  <dig> ng/ml anti-cd <dig> and the level of cflip decreased by ~50%.
the next prediction, which we observed for skw  <dig>  cells, considers a delay of caspase- <dig> activation caused by low concentrations of anti-cd <dig> .
the maximal delay predicted by the bentele’s model was about  <dig> min, whereas the composite model detected a much lower value .
analysis of dependence of caspase- <dig> concentration on initial values of iap and cd95l  also demonstrated differences in the models’ predictions.
thus, low level of iap  in the original model results in complete cell death.
however, under the same condition the composite model shows a smooth increase of caspase- <dig> levels with cd95l changing from  <dig>  nm to  <dig> nm.
if the amount of the ligand is less than  <dig>  nm, iap blocks apoptosis completely.
in addition, the bentele’s model predicts that high concentration of iap prevents a significant increase of caspase- <dig> even for high levels of the ligand, whereas in the composite model high concentration of cd95l leads to cell death.
limitations of the composite model
we analyzed a set of parametric constraints which allowed us to reduce the precursor models .
following them, we can recover a detailed description of the investigated process.
connection of the composite model with the original model by neumann et al.
was preserved, whereas with the original model by bentele et al.
it was partially lost since the chains of species interactions were changed for the sake of combining the models together.
in particular, the pathway of caspase- <dig> activation in disc was derived from the neumann’s model as well as the reaction chain of procaspase- <dig> cleavage was modified.
however, following the constraints of additional file 1: table 10s, we can complete the composite model with such reactions as bid truncation by caspase- <dig> and parp cleavage triggered by caspase- <dig>  we can also recover the details of caspase- <dig> activation mediated by cytochrome c and smac release from mitochondria but with altered kinetics of release.
for now, the question of admissibility of such modifications as well as extension of the composite model based on other apoptosis models remains open and is a challenge for further research.
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  <cit> .
in particular, the researcher has a clear vision of the most important biochemical reactions taking place in the modeled system as well as better understanding of them.
low-dimensional models are easier to analyze and faster to simulate.
this helps to save time and enhances productivity.
the main limitation of model reduction consists in loss of biological information.
however, it should be noted that even if some information  may be lost, it can result in a more clear understanding of the most important interactions and allows focusing on the decisive processes in the model, predictive ability of which is reasonably preserved.
hence, this limitation can be transformed into an advantage.
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  <dig>  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  <dig>  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 et al.
thus, the question of what model  is more realistic requires further experimental investigation.
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  <cit> .
this attempt revealed multiple conceptual and methodological difficulties.
the main obvious obstacles included:
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 ;
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 , their composition may be significantly easier in comparison to composition of quite different models.
in our work, the models by bentele et al.
and neumann et al.
are not directly related, as it could be expected as these models were constructed in the same research group.
for example, they use different reaction chains of caspase- <dig> activation and different values of kinetic parameters in the overlapping reactions.
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  <cit>  and tnf-α  <cit>  signaling, apoptosome-dependent caspase activation  <cit>  and p <dig> oscillation system  <cit> .
it is noteworthy that the model by laussmann et al.
<cit> , describing trail-induced activation of caspase- <dig>  emerged from the original study by bentele et al.
<cit> .
this fact may help us in our future work.
characterising the models used in our work, we can say that we saved our time working with the model by neumann et al.
<cit>  derived from the biomodels database  <cit> , and spent a lot of time to reproduce the model by bentele et al.
<cit>  that had not existed in a ready-to-use format.
thus, we would like to emphasize that for the common benefit, the systematic application of biological standards in modeling  and depositing the working models in public databases  would significantly facilitate analysis of existing biomathematical models and using them as the base for development of new composite systems.
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  <cit> .
using these methods, we composed two models describing pathways of cd95- and nf-κb-signaling in one without affecting the fit to the experimentally measured dynamics and model predictions.
for the model reduction and composition, we used the biouml software that was extended by the required methods of analysis.
the models by bentele et al.
and neumann et al.
reconstructed in biouml and represented using sbml format and sbgn notation, as well as their reduced and decomposed versions are available at http://ie.biouml.org/bioumlweb/#de=databases/the%20composite%20model%20of%20cd95%20and%20nf-kb%20signaling.
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.
supplementary material
additional file 1: table 1s
notations of parameters used in the paper and in the original models.
table 2s.
table of all reactions  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  <dig>  cells.
table 10s.
parametric constraints of the models by bentele et al.
and neumann et al.
click here for file
acknowledgements
this work was supported by fp <dig> grant  <dig> "net2drug", fp <dig> grant  <dig> "lipidomicnet", eu fp <dig>  and agence nationale de la recherche .
az is a member of the team “systems biology of cancer,” equipe labellisée par la ligue nationale contre le cancer.

 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. several aspects of the same process can be modeled using different levels of abstraction involving different reaction chains, chemical kinetics, and incomparable sets of parameters. such models are difficult to merge. meanwhile, merging is an important approach for creation of complex models. thus, development of efficient methods and software, allowing us to combine models, is the object of intense study in systems biology. in our work, we focused on the fact that, generally, complexity of models is not comparable to the volume of experimental data used to adjust their parameters. due to this fact, we turn to the methods of model reduction allowing us to minimize model’s complexity without affecting the model simulation dynamics. 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 and has already found multiple applications in systems biology, including discrete modeling and modeling of metabolic pathways. the principles of this technique have been applied in computational biology and implemented as a part of widely used pathway simulators such as biouml, biocham, copasi, and ginsim. model reduction led to new insights in mechanisms of translation regulation by micrornas and was applied for analysis of such signaling pathways as jak-stat, nf-κb, and egfr. in our investigation, we used the principles of model reduction to construct reasonably accurate minimal size approximations of two different models describing the cd signaling pathways. the first model explores pro-apoptotic properties of cd after stimulation by its natural ligand cd95l or by agonistic antibodies anti-cd implying formation of the death-inducing signaling complex. disc consists of oligomerized cd death domain-containing adaptor fas-associated molecule, procaspases- and − and two isoforms of cellular flip protein. caspase- leads to activation of effector caspase- directly or via stimulation of cytochrome c release from the mitochondria. the latter step requires formation of the apoptosome complex and activation of caspase- once activated, caspase- cleaves poly polymerase, thereby making the apoptotic process irreversible. the second model describes the state when cd not only activates the pro-apoptotic pathway, but also induces transcription factor nf-κb that is an important regulator of cell survival functioning. this is possible due to cflipl cleavage in the disc complex. the cleaved p43-flip directly interacts with the ikk complex and activates it. the activated ikk performs phosphorylation of the iκb protein and thereby frees nf-κb. 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 implemented in the biouml software, we reduced the models so that they still satisfy these data. this allowed us to simplify the overlapping components of the models and find a way for their composition. 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 s = associated with a set of variables c =, …, cm) depending on time t ∈ , t ∈  r+ and representing their concentrations, and a set of biochemical reactions with rates v =, …, vn) depending on a set of kinetic constants k. reaction rates are modeled by mass-action or michaelis-menten kinetics. a system of ordinary differential equations representing a linear combination of reaction rates is used to describe the model behavior over time: dctdt=n·vc,k,t,c0=c here n is a stoichiometric matrix of n by m. we say that css is a steady state of the system if n·vcss,k,t= limt→∞cit=ciss. model reduction implies transformation of the ode system to another one with smaller number of equations without affecting dynamics of variables c  …, cs, which is fixed by a set of experimental points ciexp at given times tij, j =   …, ri, where ri is number of such points for the concentration ci, i =   …,s. to check dynamics preservation in the course of model reduction, we consider the function of deviations defined as a normalized sum of squared differences: fdistc k=∑i=1s∑j=1riωminωi·citij−ciexptij where ωmin=miniωi, weights ωi=ri−1·∑jciexptij are calculated as mean squared values of experimental concentrations for each variable and the normalizing factor ω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 between model parameters. when these relations satisfy certain rules, we can approximate the detailed model by a simpler one. if the parameters are 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. we define the minimal complexity model as a reduced model with the minimal number of elements so that the deviation function value calculated for the reduced model is different from the original model no more than 20%. such threshold is explained by the fact that in this work we considered experimental data obtained by bentele et al. and neumann et al. using western blot technology with the standard deviation 15-20%. reduction of mathematical model complexity is achievable by different methods. description of the methods directly used in our work is provided below. removal of slow reactions. we say that reaction r is much slower than reaction r vr1t≪vr2t if maxtvr1t<k·maxtvr2t, where k=10- when such reactions do not affect experimental dynamics of species, we neglect them. note, that in some cases, it is sufficient to consider k = 10− or even k =  quasi-steady-state approximation assumes that all variables c of the model are split in two groups: basic variables cs and quasi-stationary variables cf for which we could write: dcstdt=ns·vsc,k,t,dcftdt=1ε·nf·vfc,k,t, where ɛ is small parameter. we consider the approximation ɛ →  resulting in the equation nf · vf =  thus, the chain of reactions s → i → p with a quasi-stationary species i can be reduced to the form s → p. lumping analysis refers to reducing the number of model variables by grouping some species. in particular, this can be illustrated by reactions r1:a+s1→p r2:a+s2→p with kinetic rates vr1t=k1·cat·cs1t and vr2t=k2·cat·cs2t. if k1 = k and cs1t=cs2t, then the system can be replaced by a single reaction 2a+s→p, where cst=cs1t=cs2t, cpt=cp1t=cp2t and k = k1 = k note, that in this example lumped variables cs1t, cs2t and cp1t, cp2t are linearly dependent: however, this is not required in general case. for conditions on lumpability in monomolecular reaction networks, see. removal of approximately linearly dependent variables. if variables a and b are approximately linearly dependent: at≈k·bt,t∈ t,k=const, then we can replace one of them with another in all kinetic laws of the model. simplification of the michaelis-menten kinetics. consider a reaction r of the form s –e → p, where an enzyme e converts a substrate s into product p. reaction rate of such reaction is frequently defined by the michaelis-menten formula vrt=k·cst·cetkm+cst where the enzyme concentration is a dynamic variable allowing to use the same kinetics in different regions of the phase space. if km ≫ cs, then we can reduce this formula to the form: vrt=k·cst·cetkm. on the other hand, if cs ≫ km, then can be replaced by the equation vr = k · ce. simplification of equations based on the law of mass action when one reactant dominates others. consider a reaction r of the form s1 + s2 → p with the kinetic law vrt=k·cs1t·cs2t. when cs1t≫cs2t for t ∈ , the formula can be replaced by the linear equation vrt=k·cs10·cs2t. the time t of validity of such pseudo-lineary approximation is defined as a period during which the relative change of cs1t does not exceed some ε, i.e. k∫0tcs2tdt<ε. one way to perform the model reduction is to apply the foregoing methods in the numerical order. note that the methods mr and mr are of the same type, so it does not matter which one to use first. this way assumes that we calculate the value of the distance function at each reduction step to test whether approximation of the experimental data is still within allowable limit or not. model analysis and comparison model comparison using akaike information criterion. this criterion defines the relative complexity of a model based on the goodness of experimental data fit and the number of model parameters |k|: aic=χ2+2·k. the function χ is defined by the formula with weights 1/ σij, where mean deviations σij of experimental values ciexp are calculated using smoothing spline: σij2=14∑k=−22citi,j+k−ciexpti,j+k for j + k <  and j + k > ri we assumed ci = ciexp =  models can be compared using the akaike criterion, if they approximate the same experimental data. when the aic value of one model is less than of others, we say that this model is simpler in terms of this criterion. model comparison with the mean aic. when we want to compare models approximating different sets of experimental data, we could calculate the mean aic coefficients by the formula aicmean=aicnexp where nexp determines the number of experimental points. steady-state analysis finds values of species concentrations according to the rules. it is desirable to reduce a model so that steady-state concentrations of all experimentally measured chemical species have not changed significantly. sensitivity analysis reveals steady-state concentrations response to parameter perturbations. the local sensitivities ciss, i =   …, m, for all parameters pj, j =   …, |k|, are calculated via finite difference approximations: sij=∂ciss∂pj=cisspj+Δpj−cisspjΔpj it is useful to compute the mean log sensitivity for all fitted species: smean=∑i,jlnsijnp·s,sij≠ i= …,s,j= …,np, where np is the number of all parameters |k| or the number of parameters retained after the model reduction. modular modeling considering the mathematical models of cd signaling pathways, we decomposed them into functional modules. this step allowed us to identify overlapping components of the models and simplify their analysis. we defined a module as a submodel with input, output and contact ports. the first two types of ports characterized variables calculated in one module and passed to another through a directed connection. the contact ports declared common variables of modules via undirected connections.  preliminary analysis, problems and inconsistencies in the precursor models the mathematical model by bentele et al. consists of species including a virtual variable xaa intended for quantification of “apoptotic activity” caused by active caspases- - and - and calculated by the formula dxaadt=1−xaa·k36act·ccasp3+k36act·ccasp6+k7act·ccasp7km367act+1−xaa−kdd·xaa, where k36act, k7act, km367act ∈ k. the variable xaa specifies a process of degradation in the model introduced as an exponential decay function fdegr. this function is defined by the formula kd · xaa + kds with kd,kds ∈ k for all active caspases and complexes containing their cleaved products, and by the formula kd · xaa for all other molecules besides cparp, where fdegr is constant. all species in the model are degraded with the exception of cytochrome c and smac stored in the mitochondria, whose concentrations are constant. in total, the model contains reactions including reactions based on mass action kinetics, reactions taken with kinetics of michaelis-menten, reactions of degradation, the reaction of xaa production specified above and two reactions of cytochrome c and smac release from the mitochondria modeled using a discrete event. the latter implies complete cytochrome c and smac release within minutes as soon as tbid reaches a certain level in comparison to bcl-2/bcl-xl. since the authors did not provide exact form of the release function, we proposed the sigmoid function which has the expected behavior: freleaset=1−11+expkcontr·−t+ttrigger+ ·trelease where ttrigger is the time when tbid concentration reaches a value of bcl-2/bcl-xl, trelease is the start time of release and kcontr is the contraction coefficient. prototypes of parameter notations in the original bentele’s model are listed in additional file 1: table 1s. abbreviations used in the table represent Мa – mass action kinetics and m-m – michaelis-menten kinetics. the column with rates includes rate orders for the fast/reduced activation scenarios. the model comprises species and kinetic parameters, which the authors estimated based on the experimental data obtained by western blot analysis for the human cell line skw. cells were stimulated by μg/ml and ng/ml of anti-cd and dynamics of several proteins 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- - - - in the case of μg/ml of anti-cd and procaspases- - in the case of ng/ml. degradation rates of procaspase- and caspases- 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- and specified kd equal to min- and min- for fast and reduced activation scenarios, respectively. the model by neumann et al. includes species, reactions constructed according to the law of mass action, and kinetic parameters. it reproduces experimental measurements of entities obtained using western blot technology for hela cells stably overexpressing cd95–gfp and treated with three different concentrations of agonistic anti-cd antibodies. this is the reduced model suggested by the authors for description of the cd95-mediated pathways of cell death and nf-κb activation. model reduction was based on intuition about what is important in the biochemical mechanisms. it was performed as long as the parameter estimation procedure was able to provide a good fit of the experimental data. particularly, the details of some complex formations were omitted as well as some non-measured species were removed. the authors submitted the model to the biomodels database. thus, we had no problems reproducing it. reduction of the cd95-signaling model we started the bentele’s model reduction by excluding the direct dependence of the virtual species xaa on caspases- and - for this purpose, we approximated the amount of caspases- by the linear function of caspase- concentration: ccasp7≈a·ccasp where a =  for both activation scenarios. in addition, we considered inequalities holding for the parameters in the formula: k7act·ccasp7≫k36act·ccasp3>k36act·ccasp km367act≫1−xaa, and simplified it: dxaadt=k36act+ ·k7actkm367act·1−xaa·ccasp3−kdd·xaa. after that we decomposed the model into modules corresponding to three biological steps: activation of caspase- cytochrome c-induced activation of caspase- and parp cleavage down-regulated by these caspases. below we provide detailed description of the reduction process of all modules. the process is based on the flow chart represented in the figure  if it is enough for a reader to have a general idea of model reduction without going into full details, we suggest to go directly to the last two paragraphs of this section, where you find the summary of the model reduction procedure. caspase- activation module includes species and reactions, which could be divided into three groups: cleavage of procaspase- at the disc complex, its activation triggered by caspase- and inhibition of complexes containing disc by cflipl and cflips. firstly, we eliminated slow reactions br and br according to the method mr next, we applied the quasi-steady-state analysis to the module and removed quasi-stationary intermediates cd95r:cd95l, disc:pro disc:pro disc:p43/p and disc:cflipl. thus, in particular, we got two main reactions instead of br1-br6: fast formation of the disc complex and slower activation of caspase- in this complex. further, we noticed that the consumptions of cflipl and cflips satisfied the same kinetic laws. therefore, these species could be lumped resulting in the reaction br4*: cflip + 2∙ disc + pro → cflip:disc:pro where cflip indicated two isoforms cflipl and cflips so that ccflip = ccflipl = ccflips. we approximated the concentration of caspase- by the linear function b · ccasp b = , and merged reactions br and br of caspases- activation induced by caspases- and mediated by caspases- into one reaction br3*: pro –casp → casp since the reaction br followed the michaelis-menten kinetics with the constants k and km then the reaction rate of br3* was provided by the kinetic law vbr3*=k38·cpro8·ccasp3/km38+cpro where k = · k and km = km analysis of this kinetic law for low level of cd95l ensured that cpro ≫ km in addition, in the case of the fast activation scenario the value of vbr3* was much lower than the rate of caspase- activation mediated by cd thus, we established vbr3* = k38·ccasp without significant changes in the results of the model simulation. complementing the above reduction steps, we removed the elements that were unnecessary to fit the experimental dynamics provided by bentele et al. these elements were degradation reactions of cd95l, cd95r and flip, and the blocked disc complex. as a result, we introduced the conservation law: ccd95l–ccd95r=l0–r0=const, where l and r denote initial concentrations of the ligand and receptor, respectively. solving the differential equation dccd95rdt=−klr·ccd95r·ccd95l, we obtained the analytical function for the receptor concentration: ccd95r=r0−l01−l0/r0·exp−klr·r0−l0·t the module of caspase- activation consisted of species and reactions including activation of caspase- triggered by cytochrome c, cleavage/activation of bid by caspases- and - iap inhibition of caspase- activity and binding of inhibitors by smac. based on the method mr we eliminated the last two reversible reactions, as well as the slow reaction br of the apoptosome complex dissociation. further, since this complex was a quasi-stationary intermediate and the rate of br was linearly dependent on a difference of br and br rates, we reduced the chain of procaspase- activation, using the methods mr and mr to the form: cytcstored→cytc:apaf− pro9–cytc:apaf−1→casp accordingly to bentele et al., we have the following rule of cytochrome c release from mitochondria: ccytcstored=ccytcstored0·freleaset, dccytcreleaseddt=−dccytcstoreddt here frelease satisfies the formula. taking into account and ignoring degradation of the complex cyt c:apaf- we got dccytc:apaf−1dt=−c·dccytcstoreddt where c =  was the coefficient of linearity. then, we noticed that the experimental measurements of procaspase- concentration were presented by bentele and his colleagues only for ng/ml of anti-cd in this case, degradation of procaspase- was insignificant. thus, neglecting it and considering the kinetic law kapop·cpro9·ccytc:apaf- of the second reaction in, we obtained the differential equation of procaspase- dynamics: dcpro9dt=− ·1−freleaset·ccytcstored0·kapop·cpro solving it, we found: cpro9=cpro90·freleaset·expkcontr·t−1+1expkcontr·t ·ccytcstored0·kapopkcontr. finally, analyzing reactions of bid activation, we removed the slower reaction mediated by caspase- the similar reaction involving caspase- followed the michaelis-menten kinetics with the constant km28bid ≫ cbid. therefore, we redefined the kinetics of this reaction based on the law of mass action: vbr7*=k8bid/km8bid·ccasp8·cbid. the module of parp inactivation contained species and reactions. six reactions were based on the michaelis-menten kinetics, four reactions reproduced the reversible inhibition of caspases- and - based on the law of mass action, and one corresponded to the production of xaa, simplification of which was discussed above. we deleted slow reactions of casp3:iap and casp7:iap complexes dissociation and ignored slow degradation of procaspase- and iap. then analyzing the cleavage of procaspases- and - we eliminated reactions triggered by caspase- which were slower in comparison to the similar reactions induced by caspase- for the same reason we ignored the reaction of parp inactivation by caspase-. considering the remaining reactions with the michaelis-menten kinetics, we used the following inequalities km78 ≫ cpro and km376act ≫ cparp, which, in combination with mr and, led to the following kinetic laws: vbr29=k78km78·ccasp8·cpro vbr32=k3actkm367act·ccasp3·cparp, where k3act = · k7act thus, the cd95-signaling model consisting of species, reactions and kinetic parameters was approximated by a model with species, reactions and parameters, except the constant concentration of cytochrome c. figure  shows the comparison of simulated concentrations of these models with experimental dynamics obtained by bentele et al. for the species mentioned above. since this dynamics was expressed in arbitrary units, we uniquely translated it into precise values for all proteins with non-zero initial values. however, for such molecules as caspase- p43/p tbid and cparp, precise values directly depended on the concentration levels at the end or in the middle of time series. accordingly, we computed these levels during the model reduction and recalculated experimental dots if necessary. species without experimental evidence and reactions br1*-br13* directly affect the simulated dynamics of the experimentally measured concentrations. thus, further reduction of the model by removing these elements is impossible. regarding the degradation process, we agree with bentele et al. that experimental data cannot be matched if this process is fully ignored. actually, if we remove any of the retained reactions of degradation, then excepted concentration dynamics will be impaired. thus, we found the minimal complexity approximation of the original model. “experimental data” in this table refers to the data obtained by bentele et al.. reduction of the cd95-mediated and nf-κb signaling model as was mentioned above, neumann et al. simplified their model in order to reduce the large number of free parameters. the authors noted that further simplification of the model was not possible because the model fit significantly decreased. however, we removed slow reactions using the method mr and removed the reaction of nf-κb degradation, which was not significant for reproducing the experimental data of the original model. note that we retained the slow reactions nr and nr the former of which regulates apoptosis inhibition in the case of high levels of cflipl and procaspase- in accordance with the precursor model prediction, and the latter of which will be required to combine this model with the bentele’s model hereafter. based on the method mr we also took into account the inequality cikk ≫ cp43−flip and simplified the kinetic law vnr19=kp43−flip_ikk·cikk·cp43−flip of the reaction nr to the form vnr19*=kp43−flip_ikk·cikk0·cp43−flip. therefore, we reduced the number of model species from to the number of reactions from to and the number of kinetic parameters from to. these modifications did not change the fit to the experimental data provided by neumann et al.. model composition analysis of the reduced models based on the akaike criterion confirmed that they had lower complexity than the original models. the difference between the mean aic coefficients of the models decreased by 60%, relative to the initial value. in addition, the reduced neumann’s model better approximated experimental data. therefore, we used it as the basis for the model composition. we took into account that the considered experimental data were obtained with skw and hela cells, which were shown to behave as type i and type ii cells, respectively. in this regard, we assumed that: initial species concentrations could vary for different cell lines; kinetic parameters of all reactions have the same values for both cell lines, whereas the value of kd could be regulated by various entities and, moreover, is dependent on the initial concentration of anti-cd; type i and type ii cells conform to different reaction chains of caspases- activation. we constructed the composite model of cd and nf-κb signaling in the following way. firstly, according to the assumption, we replaced reaction nr in the neumann’s model by a chain of three reactions: bid–casp8→tbid,pro9–tbid→casp pro3–casp9→casp the first reaction of this chain had the same kinetic law as br7*. the other reactions represented br6* and nr with kinetics modified using the law of mass action. such modification of br6* was necessary to get the slow increase of caspase- concentration experimentally described by neumann et al.. modification of nr consisted in substitution of ccasp for ccasp. secondly, we supplemented the neumann’s model with reaction of caspase- inhibition by iap and evaluated parameters in this reaction together with parameters in using the corresponding data by neumann et al. and optimization tools of biouml. 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- - cleaved product p43/p and caspase-. all the changes, except step agreed with assumptions and. to reproduce the species dynamics, we had to increase the initial concentration of procaspase- in the case of hela cells by an order of magnitude. finally, we combined the modules of caspase- and nf-κb activation with the modules of caspase- activation and parp inactivation. for this purpose, we modified the last two modules based on and refitted their parameters as well as the concentrations mentioned in table  using experimental data by bentele et al. for better fit we also supplemented the model by degradation reactions of some species. the resulting composite model consists of species, reactions and kinetic parameters, except for the constant concentration of ikk. the model showed a good fit between the simulation results and the experimental data. calculation of the mean aic coefficient for this model revealed that it has the same level of complexity as the reduced models. we also noted that reduction of both models did not significantly alter the steady-state concentrations. however, the composite model has different steady-states, which are more stable according to the mean sensitivities in the cases of skw cells and hela cells stimulated with ng/ml of anti-cd. the reduced models are equivalent to the precursor models with respect to their biological predictions analyzing the models constructed by neumann et al. and bentele et al., we divided predictions made by them into two groups: qualitative and quantitative. the first group concerned the model network and described mechanisms of the protein-protein interactions. it is this group to which we assigned the experimentally validated neumann's predictions that processed caspases are not required for nf-κb activation, as well as that pro-apoptotic and nf-κb pathways diverge already at disc. the corresponding network of reactions was preserved in the reduced model. thus, the prediction remained valid. the second group of predictions characterized behavior of the models after changing concentrations of some species, such as cd95l, cd95r, procaspase- and inhibitors cflipl, cflips, and iap. analyzing the predictions of the models by neumann et al. and bentele et al. listed in the tables  and respectively, we concluded that the reduction of the first of them did not alter the species dynamics. for the second reduced model, we detected some minor changes in the predictions. however, in general, the model behavior remained in a good agreement with the bentele's experiments. all the predictions marked with an asterisk were experimentally tested by neumann et al. and confirmed, unless otherwise noted. the simulation time in predictions 3– was min. the concentration of anti-cd considered by the authors in predictions 4– was ng/ml. the simulation time was min in all predictions. the apoptotic threshold in the prediction is the concentration of anti-cd after which cparp amount exceeds 10% of the initial parp level. 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. in the case of skw cell line, we estimated the value of degradation rate parameter kd based on the experimental observations by bentele et al. for 1– ng/ml of anti-cd. when performing simulation experiments of the model for this cell line, we used three different values of kd depending on whether the initial concentration of cd95l was within the range 1– ng/ml, within the range 100– ng/ml or greater than ng/ml. analyzing the model predictions, we detected differences in the behavior of the bentele’s and the composite models. the first of them describes a threshold mechanism for cd95-induced apoptosis implying that this process is completely stopped when cd95l concentration is below the critical value. however, the model predicts a sharp increase of cparp immediately when the level of cd95l exceeds the threshold. in contrast, the composite model shows a gradual increase of cparp with increase of the ligand level. to determine the apoptotic threshold in this case, we found concentration of cd95l for which cparp ratio is equal to 10% of the initial parp amount. as in the original model, this concentration was in the range of 1– ng/ml. additionally, the models revealed sensitivity of the threshold to the concentration of cflip. both of them predicted the cell death phenotype upon stimulation by ng/ml anti-cd and the level of cflip decreased by ~50%. the next prediction, which we observed for skw cells, considers a delay of caspase- activation caused by low concentrations of anti-cd. the maximal delay predicted by the bentele’s model was about min, whereas the composite model detected a much lower value. analysis of dependence of caspase- concentration on initial values of iap and cd95l also demonstrated differences in the models’ predictions. thus, low level of iap in the original model results in complete cell death. however, under the same condition the composite model shows a smooth increase of caspase- levels with cd95l changing from nm to nm. if the amount of the ligand is less than nm, iap blocks apoptosis completely. in addition, the bentele’s model predicts that high concentration of iap prevents a significant increase of caspase- even for high levels of the ligand, whereas in the composite model high concentration of cd95l leads to cell death. limitations of the composite model we analyzed a set of parametric constraints which allowed us to reduce the precursor models. following them, we can recover a detailed description of the investigated process. connection of the composite model with the original model by neumann et al. was preserved, whereas with the original model by bentele et al. it was partially lost since the chains of species interactions were changed for the sake of combining the models together. in particular, the pathway of caspase- activation in disc was derived from the neumann’s model as well as the reaction chain of procaspase- cleavage was modified. however, following the constraints of additional file 1: table 10s, we can complete the composite model with such reactions as bid truncation by caspase- and parp cleavage triggered by caspase- we can also recover the details of caspase- activation mediated by cytochrome c and smac release from mitochondria but with altered kinetics of release. for now, the question of admissibility of such modifications as well as extension of the composite model based on other apoptosis models remains open and is a challenge for further research. 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. in particular, the researcher has a clear vision of the most important biochemical reactions taking place in the modeled system as well as better understanding of them. low-dimensional models are easier to analyze and faster to simulate. this helps to save time and enhances productivity. the main limitation of model reduction consists in loss of biological information. however, it should be noted that even if some information may be lost, it can result in a more clear understanding of the most important interactions and allows focusing on the decisive processes in the model, predictive ability of which is reasonably preserved. hence, this limitation can be transformed into an advantage. 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 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 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 et al. thus, the question of what model is more realistic requires further experimental investigation. 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. this attempt revealed multiple conceptual and methodological difficulties. the main obvious obstacles included: 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; 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, their composition may be significantly easier in comparison to composition of quite different models. in our work, the models by bentele et al. and neumann et al. are not directly related, as it could be expected as these models were constructed in the same research group. for example, they use different reaction chains of caspase- activation and different values of kinetic parameters in the overlapping reactions. 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 and tnf-α signaling, apoptosome-dependent caspase activation and p oscillation system. it is noteworthy that the model by laussmann et al., describing trail-induced activation of caspase- emerged from the original study by bentele et al.. this fact may help us in our future work. characterising the models used in our work, we can say that we saved our time working with the model by neumann et al. derived from the biomodels database, and spent a lot of time to reproduce the model by bentele et al. that had not existed in a ready-to-use format. thus, we would like to emphasize that for the common benefit, the systematic application of biological standards in modeling and depositing the working models in public databases would significantly facilitate analysis of existing biomathematical models and using them as the base for development of new composite systems. finally, we believe that the composite model itself may be useful for further investigation of apoptosis.  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. using these methods, we composed two models describing pathways of cd95- and nf-κb-signaling in one without affecting the fit to the experimentally measured dynamics and model predictions. for the model reduction and composition, we used the biouml software that was extended by the required methods of analysis. the models by bentele et al. and neumann et al. reconstructed in biouml and represented using sbml format and sbgn notation, as well as their reduced and decomposed versions are available at http://ie.biouml.org/bioumlweb/#de=databases/the%20composite%20model%20of%20cd95%20and%20nf-kb%20signaling. 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. supplementary material additional file 1: table 1s notations of parameters used in the paper and in the original models. table 2s. table of all reactions 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 cells. table 10s. parametric constraints of the models by bentele et al. and neumann et al. click here for file acknowledgements this work was supported by fp grant "net2drug", fp grant "lipidomicnet", eu fp and agence nationale de la recherche. az is a member of the team “systems biology of cancer,” equipe labellisée par la ligue nationale contre le cancer.
systematic application of model reduction methods allows implementing this modeling principle and finding models of minimal complexity compatible with the data. 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. combining such models is much easier than of precursor models and leads to new model properties and predictions. 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. 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 cd signaling pathways. many mathematical models characterizing mechanisms of cell fate decisions have been constructed recently. firstly, we suggest a set of model reduction tools that can be systematically applied to a given 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. 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. 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. secondly, we suggest a notion of a minimal complexity model. we show that the reduced models lead to the same dynamical behavior of observable species and the same predictions as in the precursor models.
 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.  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 cd 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.  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. 