Model composition through model reduction: a combined model of CD95 and NFκB signaling pathways
 Elena Kutumova^{1, 2}Email author,
 Andrei Zinovyev^{3, 4, 5},
 Ruslan Sharipov^{1, 6} and
 Fedor Kolpakov^{1, 2}
DOI: 10.1186/17520509713
© Kutumova et al.; licensee BioMed Central Ltd. 2013
Received: 21 September 2012
Accepted: 5 February 2013
Published: 15 February 2013
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 nonunique choice of parameter values, nonreproducible 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 (APO1/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 [1]. 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 wellestablished technique in many fields of biochemical research and engineering. It has been used for many years in chemical kinetics (for reviews, see [2–4]) and has already found multiple applications in systems biology, including discrete modeling [5] and modeling of metabolic pathways [6, 7]. The principles of this technique have been applied in computational biology [8] and implemented as a part of widely used pathway simulators such as BioUML [9], BIOCHAM [10], COPASI [11], and GINsim [12]. Model reduction led to new insights in mechanisms of translation regulation by microRNAs [13, 14] and was applied for analysis of such signaling pathways as JAKSTAT [15], NFκB [16], and EGFR [17].
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 [18, 19]. The first model explores proapoptotic properties of CD95 after stimulation by its natural ligand CD95L or by agonistic antibodies antiCD95 implying formation of the deathinducing signaling complex (DISC) [18, 20]. DISC consists of oligomerized CD95, death domaincontaining adaptor FASassociated molecule (FADD), procaspases8 and −10, and two isoforms of cellular FLIP (CFLAR) protein (cFLIP long and cFLIP short). Caspase8 leads to activation of effector caspase3 directly (type I cells) or via stimulation of cytochrome C release from the mitochondria (type II cells) [21]. The latter step requires formation of the apoptosome complex and activation of caspase9. Once activated, caspase3 cleaves poly(ADPribose) polymerase (PARP), thereby making the apoptotic process irreversible. The second model describes the state when CD95 not only activates the proapoptotic pathway, but also induces transcription factor NFκB that is an important regulator of cell survival functioning [19]. This is possible due to cFLIPL cleavage in the DISC complex. The cleaved p43FLIP 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 [2] 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
where ${\omega}_{\mathit{min}}=\underset{\mathit{\text{i}}}{\mathit{\text{min}}}{\omega}_{i},$ weights ${\omega}_{i}=\sqrt{{r}_{i}^{1}\xb7{\displaystyle \sum}_{j}{\left({C}_{i}^{\mathit{exp}}\left({t}_{\mathit{ij}}\right)\right)}^{2}}$ 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.
We define the minimal complexity model as a reduced model with the minimal number of elements (species and reactions) so that the deviation function value (3) 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.[18] and Neumann et al.[18] using Western blot technology with the standard deviation 1520%.
Reduction of mathematical model complexity is achievable by different methods [2]. Description of the methods directly used in our work is provided below.
(MR1) Removal of slow reactions. We say that reaction r_{1} is much slower than reaction ${r}_{2},{v}_{{r}_{1}}\left(t\right)\ll {v}_{{r}_{2}}\left(t\right)$ if $\underset{\mathit{\text{t}}}{\mathit{\text{max}}}\left{v}_{{r}_{1}}\left(t\right)\right<k\xb7\underset{\mathit{\text{t}}}{\mathit{\text{max}}}\left{v}_{{r}_{2}}\left(t\right)\right$, where k=10^{2}. 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^{−1} or even k = 1.
where ɛ is small parameter. We consider the approximation ɛ → 0 resulting in the equation N^{ f } · v^{ f }(C, K, t) = 0. Thus, the chain of reactions S → I → P with a quasistationary species I can be reduced to the form S → P.
where ${C}_{S}\left(t\right)={C}_{{S}_{1}}\left(t\right)={C}_{{S}_{2}}\left(t\right)$, ${C}_{P}\left(t\right)={C}_{{P}_{1}}\left(t\right)={C}_{{P}_{2}}\left(t\right)$ and k = k_{1} = k_{2}. Note, that in this example lumped variables ${C}_{{S}_{1}}\left(t\right)$, ${C}_{{S}_{2}}\left(t\right)$ and ${C}_{{P}_{1}}\left(t\right)$, ${C}_{{P}_{2}}\left(t\right)$ are linearly dependent: however, this is not required in general case. For conditions on lumpability in monomolecular reaction networks, see [23].
then we can replace one of them with another in all kinetic laws of the model.
where the enzyme concentration is a dynamic variable allowing to use the same kinetics in different regions of the phase space. If Km ≫ C_{ S }(t), then we can reduce this formula to the form: ${v}_{r}\left(t\right)=\frac{k\xb7{C}_{S}\left(t\right)\xb7{C}_{E}\left(t\right)}{\mathit{Km}}$. On the other hand, if C_{ S }(t) ≫ Km, then (5) can be replaced by the equation v_{ r }(t) = k · C_{ E }(t).
When ${C}_{{S}_{1}}\left(t\right)\gg {C}_{{S}_{2}}\left(t\right)$ for t ∈ [0, T], the formula (6) can be replaced by the linear equation ${v}_{r}\left(t\right)=k\xb7{C}_{{S}_{1}}\left(0\right)\xb7{C}_{{S}_{2}}\left(t\right)$. The time T of validity of such pseudolineary approximation is defined as a period during which the relative change of ${C}_{{S}_{1}}\left(t\right)$ does not exceed some ε, i.e. $k\underset{0}{\overset{T}{{\displaystyle \int}}}{C}_{{S}_{2}}\left(t\right)\mathit{dt}<\mathit{\epsilon}$.
Model analysis and comparison
For j + k < 0 and j + k > r_{ i } we assumed C_{ i }(t_{i,j+k}) = C_{ i }^{ exp }(t_{i,j+k}) = 0.
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.
where n_{ exp } determines the number of experimental points.
Steadystate analysis finds values of species concentrations according to the rules (2). It is desirable to reduce a model so that steadystate concentrations of all experimentally measured chemical species have not changed significantly.
where n_{ p } is the number of all parameters K or the number of parameters retained after the model reduction.
Modular modeling
Considering the mathematical models of CD95 signaling pathways [18, 19], 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 (including several species, reactions and parameters of the model) 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 (for more details of modular modeling, see [26]).
Results
Preliminary analysis, problems and inconsistencies in the precursor models
Summary of reactions from the original model by Bentele et al. and the reduced model
Original model  Reduced model  

№  Reactions (Kinetics)  Rates (nM/min)  №  Reactions (Kinetics) 
Activation of caspase8 induced by CD95  
br1  CD95L + CD95R → СD95R:СD95L (MA)  10^{1}/10^{0}  br1*  CD95L → DISC (k_{ LR } · C_{CD95R} · C_{CD95L}, C_{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 (k_{DISC_pro8} · C_{pro8} · C_{ DISC }) 
br6  DISC:p43/p41 → casp8 + DISC (MA)  10^{0}/10^{1}  
Activation of caspase8 by caspase3  
br7  pro6 casp3 → casp6 (MM)  10^{0}/10^{1}  br3*  pro8 –casp3 → casp8 (k_{38} · C_{casp3}) 
br8  pro8 casp6 → casp8 (MM)  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 (k_{DISC_FLIP} · C_{ FLIP } · C_{ 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 (k_{DFp8} · C_{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 caspase9 triggered by cytochrome C  
br15  Cyt C stored → Cyt C (f_{ release }(t)·C_{ CytCstored })  10^{2}  br6*  pro9 → cleavage (the formula (16)) 
br16  Apaf1 + Cyt C → Cyt C:Apaf1 (MA)  10^{2}  
br17  Cyt C:Apaf1 → Apaf1 + Cyt C (MA)  10^{1}  
br18  pro9 + Cyt C:Apaf → Apop (MA)  10^{0}  
br19  Apop casp3 → casp9 + Cyt C:Apaf1  10^{0}  
(MM)  
br20  Apop → casp9 + Cyt C:Apaf1 (MA)  10^{2}  
Activation of Bid  
br21  pro2 –casp3 → casp2 (MM)  10^{0}  br7*  Bid –casp8 → tBid ($\frac{{k}_{8\mathit{Bid}}}{K{m}_{8\mathit{Bid}}}\xb7{C}_{\mathit{casp}8}\xb7{C}_{\mathit{Bid}}$) 
br22  Bid –casp8 → tBid (MM)  10^{0}  
br23  Bid –casp2 → tBid (MM)  10^{1}/10^{0}  br8*  pro2 –casp3 → cleavage ($\frac{{k}_{32}\xb7{C}_{\mathit{casp}3}\xb7{C}_{\mathit{pro}2}}{{C}_{\mathit{pro}2}+K{m}_{32}}$) 
Blocking of IAP by Smac  
br24  Smac stored → Smac (f_{ release }(t)·C_{ Smac stored })  10^{2}  –  
br25  Smac + IAP → IAP:Smac (MA)  10^{5}  
br26  IAP:Smac → Smac + IAP (MA)  10^{5}  
Activation of caspases3 and −7  
br27  pro3 –casp8 → casp3 (MM)  10^{1}/10^{0}  br9*  pro3 –casp8 → casp3 (MM) 
br28  pro3 –casp9 → casp3 (MM)  10^{3}/10^{0}  
br29  pro7 –casp8 → casp7 (MM)  10^{0}/10^{1}  br10*  pro7 –casp8 → cleavage ($\frac{{k}_{78}}{K{m}_{78}}\xb7{C}_{\mathit{casp}8}\xb7{C}_{\mathit{pro}7}$) 
br30  pro7 –casp9 → casp7 (MM)  10^{5}/10^{3}  
PARP inactivation  
br31  PARP casp3 → cPARP (ММ)  10^{1}  br11*  PARP casp3 → cPARP ($\frac{{k}_{3\mathit{act}}}{K{m}_{367\mathit{act}}}\xb7{C}_{\mathit{casp}3}\xb7{C}_{\mathit{PARP}}$) 
br32  PARP –casp7 → cPARP (ММ)  10^{0}/10^{1}  
Inhibition of caspases3, 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 x _{ aa }  
br39  → x_{ aa } (1 − x_{ aa }/Km_{367act} + (1 − x_{ aa }) · (k_{36act} · C_{casp3} + k_{36act} · C_{casp6} + k_{7act} · C_{casp7}))  10^{0}  br13*  → x_{ aa } ($\frac{{k}_{36\mathit{act}}+0.18\xb7{k}_{7\mathit{act}}}{K{m}_{367\mathit{act}}}\xb7\left(1{x}_{\mathit{aa}}\right)\xb7{C}_{\mathit{casp}3}$) 
The model comprises 43 species (including x_{ 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 antiCD95 (fast and reduced activation scenarios, respectively) and dynamics of several proteins (Bid, tBid, PARP, cPARP, procaspases2, 3, 7, 8, 9, cleaved product of procaspase8 p43/p41, and caspases8) 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 procaspases2, 3, 7, 8 in the case of 5 μg/ml of antiCD95 and procaspases2, 9 in the case of 200 ng/ml. Degradation rates of procaspase8 and caspases8 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 k_{ 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.
Reduction of the CD95signaling model
After that we decomposed the model into modules (Figure 3B) corresponding to three biological steps: activation of caspase8, cytochrome Cinduced activation of caspase9 and PARP cleavage downregulated 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 2. 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.
Caspase8 activation module includes 20 species and 14 reactions (besides the reactions of degradation), which could be divided into three groups: cleavage of procaspase8 at the DISC complex, its activation triggered by caspase3, and inhibition of complexes containing DISC by cFLIPL and cFLIPS (Table 1, Figure 3A, reactions br1br14).
Firstly, we eliminated slow reactions br12 and br13 according to the method MR1. Next, we applied the quasisteadystate analysis (MR2) to the module and removed quasistationary intermediates CD95R:CD95L, DISC:pro8, DISC:pro8_{2}, DISC:p43/p41 and DISC:cFLIPL. Thus, in particular, we got two main reactions instead of br1br6: fast formation of the DISC complex (Table 1, Figure 3B, br1*) and slower activation of caspase8 in this complex (br2*).
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 C_{ cFLIP } = C_{ cFLIPL } = C_{ cFLIPS }.
We approximated the concentration of caspase6 by the linear function b · C_{casp3}, b = 0.145 (MR4), and merged reactions br7 and br8 of caspases8 activation induced by caspases3 and mediated by caspases6 into one reaction br3*:
pro8 –casp3 → casp8.
where k_{38} = 0.145 · k_{68} and Km_{38} = Km_{68}. Analysis of this kinetic law for low level of CD95L ensured that C_{pro8} ≫ Km_{38}. In addition, in the case of the fast activation scenario the value of v_{ br3* } was much lower than the rate of caspase8 activation mediated by CD95. Thus, we established v_{ br3* } = k_{ 38 }·C_{ casp3 } without significant changes in the results of the model simulation (MR1, MR5).
where c = 0.59 was the coefficient of linearity.
The module of PARP inactivation contained 13 species (including the virtual species x_{ aa }) and 11 reactions (excepting reactions of degradation). Six reactions (activation of caspases3 and 7 by caspases8 and 9, and cleavage/inactivation of PARP) were based on the MichaelisMenten kinetics (br27br32), four reactions reproduced the reversible inhibition of caspases3 and 7 based on the law of mass action (br33, br34, br36, br37), and one corresponded to the production of x_{ 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 procaspase7 and IAP. Then analyzing the cleavage of procaspases3 and 7, we eliminated reactions triggered by caspase9, which were slower in comparison to the similar reactions induced by caspase8. For the same reason we ignored the reaction of PARP inactivation by caspase3 (MR1).
where k_{ 3act } = 0.18 · k_{ 7act }
The role of species and reactions in the apoptotic process described by the reduced model
Species  Reactions  The role in the apoptosis process 

CD95L  br1*  Triggering the apoptosis process. 
DISC  br2*  Activation of caspase8 in the DISC complex. In the case of a high ligand concentration, DISC is a slow species, so reactions br1* and br2* cannot be combined. 
caspase3  br3*  Activation of caspase8 via the negativefeedback loop. When the concentration of CD95L is reduced, the rate of br3* is an order of magnitude higher than the rate of br2*. Thus, activation of caspase8 is mainly caused by caspase3. 
br8*  Cleavage/activation of procaspase2 is confirmed by the experimental data.  
br9*  Caspase3 activation plays a crucial role in the cleavage of procaspases2, 8 and PARP.  
br11*  PARP inactivation is confirmed by the experimental measurements of PARP and cPARP concentrations.  
cFLIP  br4*  Inhibition of the DISC complex. cFLIP blocks the activity of DISC preventing a significant increase of caspase8 concentration. This conclusion is consistent with the observations by Bentele et al. asserting that downregulation of cFLIP resulted in cell death occurring already upon low concentration of antiCD95 (1 ng/ml). 
DISC:cFLIP:pro8  br5*  Production of blocked p43/p41, whose concentration is detected by the experimental data. 
br6*, br7*, br10*  Activation of Bid and cleavage of procaspases7 and −9. These reactions determine dynamics of the mentioned species in agreement with the experimental data.  
IAP  br12*  Inhibition of caspase3. In the case of the reduced activation scenario, IAP prevents caspase3 from reaching a significant level and, therefore, blocks significant increase of caspase8 due to br3*. 
x _{ aa }  br13*  The virtual variable x_{ aa } specifies the process of degradation. 
Reduction of the CD95mediated 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 (nr4, nr11nr13, Figure 4A) using the method MR1, and removed the reaction of NFκB degradation (nr23), which was not significant for reproducing the experimental data of the original model. Note that we retained the slow reactions nr6 and nr7, the former of which regulates apoptosis inhibition in the case of high levels of cFLIPL and procaspase8 in accordance with the precursor model prediction (see the analysis of the predictions below), and the latter of which will be required to combine this model with the Bentele’s model hereafter.
Model composition
Comparison of the investigated models according to Akaike’s information criterion ( AIC )
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 
 (A)
initial species concentrations could vary for different cell lines;
 (B)
kinetic parameters of all reactions (besides the degradation rate modeled as the function k _{ d } · x _{ aa } ^{2}) have the same values for both cell lines, whereas the value of k _{ d } could be regulated by various entities and, moreover, is dependent on the initial concentration of antiCD95 [18];
 (C)
type I and type II cells conform to different reaction chains of caspases3 activation [21].
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 caspase3 concentration experimentally described by Neumann et al.[19]. Modification of nr15 consisted in substitution of C_{casp9} for C_{casp8} .
Secondly, we supplemented the Neumann’s model with reaction of caspase3 inhibition by IAP and evaluated parameters in this reaction together with parameters in (17) using the corresponding data by Neumann et al. and optimization tools of BioUML [29].
Modifications of the reduced Neumann’s model to reproduce the experimental data by Bentele et al
Step  Changed parameters and concentrations  Initial values, modifications  Reason  Ranges and initial values of fitting 

1  CD95R:FADD  91.266 nM, increase  Case: antiCD95 = 5 μg/ml Procaspase8 is cleaved in about 30 minutes. This is possible only if DISC concentration reaches a level sufficient to accelerate the reaction nr2.  [10^{2}, 10^{3}], 442.821, Bentele et al. 
2  procaspase8  64.477 nM, increase  Case: antiCD95 = 5 μg/ml The concentration of caspase8 should reach its maxim value in about 20 minutes.  [10^{2}, 10^{3}], 442.821, Bentele et al. 
3  FLIPS  5.084 nM, increase  Case: antiCD95 = 200 ng/ml As mentioned in Table 2, activation of caspase8 is mainly caused by caspase3 and, therefore, is delayed approximately for 30 minutes. Since reactions nr2 and nr5 have the same rates, we cannot observe the delay, but when we reduce the rate of nr5, we can. The latter is achieved by upregulation of the DISC:pro8 complex inhibition.  [10^{1}, 10^{2}], 65.021, Bentele et al. 
4  k10  0.121 nM^{1}min^{1}, decrease by an order of magnitude  Case: antiCD95 = 200 ng/ml Preventing a rapid growth of caspase8 concentration.  0.012 
5  procaspase3  1.443 nM, refitting  Case: antiCD95 = 5 μg/ml Improving approximation of the experimental data by Bentele et al.  [10^{0}, 10^{1}], 1.443, Neumann et al. 
6  Bid, procaspase9  5.003 nM, 2.909 nM, replacement by Bentele’s values  Case: antiCD95 = 200 ng/ml Reproducing the experimental dynamics of procaspase9.  231.760 (Bid), 245.101 (pro9), Bentele et al. 
Calculation of the mean AIC coefficient for this model revealed that it has the same level of complexity as the reduced models (Table 3). We also noted that reduction of both models did not significantly alter the steadystate concentrations (Additional file 1: Table 5S, Additional file 1: Table6S). However, the composite model has different steadystates, which are more stable according to the mean sensitivities (9) in the cases of SKW 6.4 cells and HeLa cells stimulated with 250 ng/ml of antiCD95 (Additional file 1: Table 7S).
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 proteinprotein 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 proapoptotic and NFκB pathways diverge already at DISC. The corresponding network of reactions was preserved in the reduced model. Thus, the prediction remained valid.
Analysis of predictions regarding apoptosis in HeLa cells as formulated by Neumann et al
№  The composite model behavior  Predictions by Neumann et al. 

1* 
 The concentration of antiCD95 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 12fold. 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κBsignaling pathways. To test this prediction, Neumann et al. considered levels of caspases8 cleavage and IκBα degradation for the original (solid lines) amount of CD95 and the amount decreased by about 12fold (dashed lines). The concentration of antiCD95 was 500 ng/ml. 
3* 
 Along with increasing the concentration of antiCD95 from 500 ng/ml to 1500 ng/ml, p43/p41 peaks earlier, while there is almost no difference for p43FLIP. 
4 
 Increased concentrations of cFLIPS inhibit both apoptotic and NFκB pathways, although p43FLIP generation is inhibited at a lower threshold than p43/p41 generation. 
5* 
 Increasing the concentration of cFLIPL leads to a steep increase in p43FLIP generation until it reaches a maximum, after which the curve drops. Lowered levels of cFLIPL result in very little p43FLIP but almost unchanged levels of p43/p41. 
At very high concentrations of cFLIPL no p43FLIP is generated. This dropoff was not observed experimentally by the authors.  
6* 
 Only an intermediate level of cFLIPL promotes NFkB activation. Decreased levels of procaspase8 lead to a significantly lower amount of p43FLIP and, subsequently, NFκB. The figures show of logarithmic dependence of the maximal NFκB concentration on the initial values of procaspase8 and cFLIPL 
7* 
 High cFLIPL or low procaspase8 concentrations cause suppression of apoptosis. The figures show the same dependence as considered in the previous prediction, but with caspases3 instead of NFκB. 
Predictions of the models for SKW 6.4 cells
№  The models behavior  Experimental observations by Bentele et al.and predictions of the models 

1 
 Experiments by Bentele et al.: 
− for 1 ng/ml of antiCD95, PARP cleavage was not observed;  
− the measured death rate for 10 ng/ml of antiCD95 was 2030%.  
Original model (red):  
− the apoptotic threshold is 1.9 ng/ml;  
− cPARP concentration rises dramatically within an extremely narrow interval of antiCD95 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 antiCD95 level.  
2 
 Experiments by Bentele et al.: downregulation of cFLIP in SKW 6.4 cells by addition of cyclohexamide resulted in cell death (40% for 1 day) already upon 1 ng/ml of antiCD95. The level of cFLIP was decreased to 70%. 
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 antiCD95.  
3 
 Experiments by Bentele et al.: in the 10 ng/ml activation scenario, a significant increase of caspase8 was observed after more than 4 hours. 
Original (red) and composite (blue) models: antiCD95 concentrations which are slightly above the apoptotic threshold result in caspase8 activation after a delay of many hours.  
The figure shows peak times of caspase8 concentration exceeding 0.1% of the initial procaspase8 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 caspase3 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 caspases3 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 5).
In the case of SKW 6.4 cell line, we estimated the value of degradation rate parameter k_{ d } based on the experimental observations by Bentele et al. for 1–10 ng/ml of antiCD95 (Table 6). When performing simulation experiments of the model for this cell line, we used three different values of k_{ 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 1: Table 3S).
Analyzing the model predictions, we detected differences in the behavior of the Bentele’s and the composite models (Table 6). The first of them describes a threshold mechanism for CD95induced 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 (Table 6, № 1). 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 [30]. As in the original model, this concentration was in the range of 1–10 ng/ml. Additionally, the models revealed sensitivity of the threshold to the concentration of cFLIP (Table 6, № 2). Both of them predicted the cell death phenotype upon stimulation by 1 ng/ml antiCD95 and the level of cFLIP decreased by ~50%.
The next prediction, which we observed for SKW 6.4 cells, considers a delay of caspase8 activation caused by low concentrations of antiCD95 (Table 6, № 3). The maximal delay predicted by the Bentele’s model was about 2500 min, whereas the composite model detected a much lower value (700–800 min). Analysis of dependence of caspase3 concentration on initial values of IAP and CD95L (Table 6, № 4) also demonstrated differences in the models’ predictions. Thus, low level of IAP (less than 1 nM) in the original model results in complete cell death. However, under the same condition the composite model shows a smooth increase of caspase3 levels with CD95L changing from 0.3 nM to 1000 nM. If the amount of the ligand is less than 0.3 nM, IAP blocks apoptosis completely. In addition, the Bentele’s model predicts that high concentration of IAP prevents a significant increase of caspase3 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 (Additional file 1: Table 10S). 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 caspase8 activation in DISC was derived from the Neumann’s model as well as the reaction chain of procaspase9 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 caspase2 and PARP cleavage triggered by caspase7. We can also recover the details of caspase9 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 lowdimensional models [8]. 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. Lowdimensional 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 (for example, about very slow biochemical reactions) 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 p43FLIP 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 antiapoptotic 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 et al. Thus, the question of what model (original or composite) is more realistic requires further experimental investigation.

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. singlecell 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 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 caspase8 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 [31–33] and TNFα [34, 35] signaling, apoptosomedependent caspase activation [36] and p53 oscillation system [37]. It is noteworthy that the model by Laussmann et al.[33], describing TRAILinduced activation of caspase8, emerged from the original study by Bentele et al.[18]. 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.[19] derived from the BioModels database [27], and spent a lot of time to reproduce the model by Bentele et al.[18] that had not existed in a readytouse format. Thus, we would like to emphasize that for the common benefit, the systematic application of biological standards in modeling (e.g. SBML [38] and SBGN [39]) and depositing the working models in public databases (e.g. BioModels [27]) 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 [26]. Using these methods, we composed two models describing pathways of CD95 and NFκBsignaling 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%20NFkB%20signaling.
Declarations
Acknowledgements
This work was supported by FP6 grant 037590 "Net2Drug", FP7 grant 090107 "LipidomicNet", EU FP7 (project APOSYS) and Agence Nationale de la Recherche (project ANR08SYSC003 CALAMAR). AZ is a member of the team “Systems Biology of Cancer,” Equipe labellisée par la Ligue Nationale Contre le Cancer.
Authors’ Affiliations
References
 Hlavacek WS: How to deal with large models?. Mol Syst Biol. 2009, 5: 240PubMedPubMed CentralView ArticleGoogle Scholar
 Gorban AN, Radulescu O, Zinovyev AY: Asymptotology of Chemical Reaction Networks. Chemical Engineering Science. 2009, 65: 23102324.View ArticleGoogle Scholar
 Gorban AN, Karlin IV, Zinovyev AY: Constructive methods of invariant manifolds for kinetic problems. Physics Reports. 2004, 396: 197403.View ArticleGoogle Scholar
 Klonowski W: Simplifying principles for chemical and enzyme reaction kinetics. Biophys Chem. 1983, 18 (2): 7387.PubMedView ArticleGoogle Scholar
 Calzone L, Tournier L, Fourquet S, Thieffry D, Zhivotovsky B, Barillot E, Zinovyev A: Mathematical modelling of cellfate decision in response to death receptor engagement. PLoS Comput Biol. 2010, 6 (3): e1000702PubMedPubMed CentralView ArticleGoogle Scholar
 Liebermeister W, Baur U, Klipp E: Biochemical network models simplified by balanced truncation. FEBS J. 2005, 272 (16): 40344043.PubMedView ArticleGoogle Scholar
 Bruggeman FJ, Westerhoff HV, Hoek JB, Kholodenko BN: Modular response analysis of cellular regulatory networks. J Theor Biol. 2002, 218 (4): 507520.PubMedView ArticleGoogle Scholar
 Radulescu O, Gorban AN, Zinovyev A, Noel V: Reduction of dynamical biochemical reactions networks in computational biology. Frontiers in genetics. 2012, 3: 131PubMedPubMed CentralView ArticleGoogle Scholar
 Kolpakov FA, Valeev IT, Tolstyh NI, Kisele I, Kondrakhin YV, Kutumova EO, Yevshin IS: BIoUML  open source Integrated platform for building virtual cell and virtual physiological 1 human. Virtual Biology. in press, DOl: 10.127041vb1e12
 Calzone L, Fages F, Soliman S: BIOCHAM: an environment for modeling biological systems and formalizing experimental knowledge. Bioinformatics. 2006, 22 (14): 18051807.PubMedView ArticleGoogle Scholar
 Hoops S, Sahle S, Gauges R, Lee C, Pahle J, Simus N, Singhal M, Xu L, Mendes P, Kummer U: COPASI — a COmplex PAthway SImulator. Bioinformatics. 2006, 22: 30673074.PubMedView ArticleGoogle Scholar
 Naldi A, Berenguier D, Fauré A, Lopez F, Thieffry D, Chaouiya C: Logical modelling of regulatory networks with GINsim 2.3. Biosystems. 2009, 97 (2): 13413.PubMedView ArticleGoogle Scholar
 Zinovyev A, Morozova N, Nonne N, Barillot E, HarelBellan A, Gorban AN: Dynamical modeling of microRNA action on the protein translation process. BMC Syst Biol. 2010, 4: 13PubMedPubMed CentralView ArticleGoogle Scholar
 Morozova N, Zinovyev A, Nonne N, Pritchard LL, Gorban AN, HarelBellan A: Kinetic signatures of microRNA modes of action. RNA. 2012, 18 (9): 16351655.PubMedPubMed CentralView ArticleGoogle Scholar
 Quaiser T, Dittrich A, Schaper F, Mönnigmann M: A simple workflow for biologically inspired model reduction – application to early JAKSTAT signaling. BMC Syst Biol. 2011, 5: 30PubMedPubMed CentralView ArticleGoogle Scholar
 Radulescu O, Gorban AN, Zinovyev A, Lilienbaum A: Robust simplifications of multiscale biochemical networks. BMC Syst Biol. 2008, 2: 86PubMedPubMed CentralView ArticleGoogle Scholar
 Conzelmann H, SaezRodriguez J, Sauter T, Bullinger E, Allgöwer F, Gilles ED: Reduction of mathematical models of signal transduction networks: simulationbased approach applied to EGF receptor signalling. Syst Biol. 2004, 1 (1): 159169.View ArticleGoogle Scholar
 Bentele M, Lavrik I, Ulrich M, Stößer S, Heermann DW, Kalthoff H, Krammer PH, Eils R: Mathematical modeling reveals threshold mechanism in CD95induced apoptosis. J Cell Biol. 2004, 166 (9): 839851.PubMedPubMed CentralView ArticleGoogle Scholar
 Neumann L, Pforr C, Beaudouin J, Pappa A, Fricker N, Krammer PH, Lavrik IN, Eils R: Dynamics within the CD95 deathinducing signaling complex decide life and death of cells. Mol Syst Biol. 2010, 6: 352PubMedPubMed CentralView ArticleGoogle Scholar
 Krammer PH, Rüdiger A, Lavrik IN: Life and death in peripheral T cells. Nat Rev Immunol. 2007, 7: 532542.PubMedView ArticleGoogle Scholar
 Scaffidi C, Fulda S, Srinivasan A, Friesen C, Li F, Tomaselli KJ, Debatin KM, Krammer PH, Peter ME: Two CD95 (APO1/Fas) signaling pathways. EMBO J. 1998, 17 (6): 16751687.PubMedPubMed CentralView ArticleGoogle Scholar
 Choi J, Yang KW, Lee TY, Lee SY: New timescale criteria for model simplification of bioreaction systems. BMC Bioinformatics. 2008, 9: 338PubMedPubMed CentralView ArticleGoogle Scholar
 Wei J, Kuo JC: A lumping analysis in monomolecular reaction systems: Analysis of the exactly lumpable system. Ind Eng Chem Fundam. 1969, 8: 114123.View ArticleGoogle Scholar
 Akaike H: A new look at the statistical model identification. IEEE Transactions on Automatic Control. 1974, 19 (6): 716723.View ArticleGoogle Scholar
 Schilling M, Maiwald T, Hengl S, Winter D, Kreutz C, Kolch W, Lehmann WD, Timmer J, Klingmüller U: Theoretical and experimental analysis links isoformspecific ERK signaling to cell fate decisions. Mol Syst Biol. 2009, 5: 334PubMedPubMed CentralView ArticleGoogle Scholar
 Kutumova EO, Kiselev IN, Sharipov RN, Lavrik IN, Kolpakov FA: A modular model of the apoptosis machinery. Adv Exp Med Biol. 2012, 736 (2): 235245.PubMedView ArticleGoogle Scholar
 Le Novère N, Bornstein B, Broicher A, Courtot M, Donizelli M, Dharuri H, Li L, Sauro H, Schilstra M, Shapiro B, Snoep JL, Hucka M: BioModels Database: a free, centralized database of curated, published, quantitative kinetic models of biochemical and cellular systems. Nucleic Acids Res. 2006, 34: D689D691.PubMedPubMed CentralView ArticleGoogle Scholar
 Engels IH, Stepczynska A, Stroh C, Lauber K, Berg C, Schwenzer R, Wajant H, Jänicke RU, Porter AG, Belka C, Gregor M, SchulzeOstho K, Wesselborg K: Caspase8/FLICE functions as an executioner caspase in anticancer druginduced apoptosis. Oncogene. 2000, 19: 45634573.PubMedView ArticleGoogle Scholar
 Kutumova EO, Ryabova AS, Valeev TF, Kolpakov FA: BI0UML plugin for nonlinear parameter estimation using multiple experimental data. Virtual Biology. 1: in press. DOl: 10.127041vb1e10
 Gaudet S, Spencer SL, Chen WW, Sorger PK: Exploring the contextual sensitivity of factors that determine celltocell variability in receptormediated apoptosis. PLoS Comput Biol. 2012, 8 (4): e1002482PubMedPubMed CentralView ArticleGoogle Scholar
 Albeck JG, Burke JM, Aldridge BB, Zhang M, Lauffenburger DA, Sorger PK: Quantitative analysis of pathways controlling extrinsic apoptosis in single cells. Mol Cell. 2008, 30 (1): 1125.PubMedPubMed CentralView ArticleGoogle Scholar
 Albeck JG, Burke JM, Spencer SL, Lauffenburger DA, Sorger PK: Modeling a snapaction, variabledelay switch controlling extrinsic cell Death. PLoS Biol. 2008, 6 (12): e299PubMed CentralView ArticleGoogle Scholar
 Laussmann MA, Passante E, Hellwig CT, Tomiczek B, Flanagan L, Prehn JH, Huber HJ, Rehm M: Proteasome inhibition can impair caspase8 activation upon submaximal stimulation of apoptotic tumor necrosis factorrelated apoptosis inducing ligand (TRAIL) signaling. J Biol Chem. 2012, 287 (18): 1440214411.PubMedPubMed CentralView ArticleGoogle Scholar
 Rangamani P, Sirovich L: Survival and apoptotic pathways initiated by TNFalpha: modeling and predictions. Biotechnol Bioeng. 2007, 97 (5): 12161229.PubMedView ArticleGoogle Scholar
 Cho KH, Shin SY, Lee HW, Wolkenhauer O: Investigations into the analysis and modeling of the TNFαmediated NFκBsignaling pathway. Genome Res. 2003, 13: 24132422.PubMedPubMed CentralView ArticleGoogle Scholar
 Rehm M, Huber HJ, Dussmann H, Prehn JH: Systems analysis of effector caspase activation and its control by Xlinked inhibitor of apoptosis protein. EMBO J. 2006, 25 (18): 43384349.PubMedPubMed CentralView ArticleGoogle Scholar
 Hamada H, Tashima Y, Kisaka Y, Iwamoto K, Hanai T, Eguchi Y, Okamoto M: Sophisticated framework between cell cycle arrest and apoptosis induction based on p53 dynamics. PLoS One. 2008, 4 (3): e4795View ArticleGoogle Scholar
 Hucka M, Finney A, Sauro HM, Bolouri H, Doyle JC, Kitano H, Arkin AP, Bornstein BJ, Bray D, CornishBowden A, Cuellar AA, Dronov S, Gilles ED, Ginkel M, Gor V, Goryanin II, Hedley WJ, Hodgman TC, Hofmeyr JH, Hunter PJ, Juty NS, Kasberger JL, Kremling A, Kummer U, Le Novère N, Loew LM, Lucio D, Mendes P, Minch E, Mjolsness ED: The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics. 2003, 19 (4): 524531.PubMedView ArticleGoogle Scholar
 Le Novère N, Hucka M, Mi H, Moodie S, Schreiber F, Sorokin A, Demir E, Wegner K, Aladjem MI, Wimalaratne SM, Bergman FT, Gauges R, Ghazal P, Kawaji H, Li L, Matsuoka Y, Villéger A, Boyd SE, Calzone L, Courtot M, Dogrusoz U, Freeman TC, Funahashi A, Ghosh S, Jouraku A, Kim S, Kolpakov F, Luna A, Sahle S, Schmidt E: The Systems Biology Graphical Notation. Nat Biotechnol. 2009, 27 (8): 735741.PubMedView ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.