Decoding complex biological networks - tracing essential and modulatory parameters in complex and simplified models of the cell cycle

Background One of the most well described cellular processes is the cell cycle, governing cell division. Mathematical models of this gene-protein network are therefore a good test case for assessing to what extent we can dissect the relationship between model parameters and system dynamics. Here we combine two strategies to enable an exploration of parameter space in relation to model output. A simplified, piecewise linear approximation of the original model is combined with a sensitivity analysis of the same system, to obtain and validate analytical expressions describing the dynamical role of different model parameters. Results We considered two different output responses to parameter perturbations. One was qualitative and described whether the system was still working, i.e. whether there were oscillations. We call parameters that correspond to such qualitative change in system response essential. The other response pattern was quantitative and measured changes in cell size, corresponding to perturbations of modulatory parameters. Analytical predictions from the simplified model concerning the impact of different parameters were compared to a sensitivity analysis of the original model, thus evaluating the predictions from the simplified model. The comparison showed that the predictions on essential and modulatory parameters were satisfactory for small perturbations, but more discrepancies were seen for larger perturbations. Furthermore, for this particular cell cycle model, we found that most parameters were either essential or modulatory. Essential parameters required large perturbations for identification, whereas modulatory parameters were more easily identified with small perturbations. Finally, we used the simplified model to make predictions on critical combinations of parameter perturbations. Conclusions The parameter characterizations of the simplified model are in large consistent with the original model and the simplified model can give predictions on critical combinations of parameter perturbations. We believe that the distinction between essential and modulatory perturbation responses will be of use for sensitivity analysis, and in discussions of robustness and during the model simplification process.


Background
In recent years, numerous studies characterizing static topological properties of molecular networks derived from heterogeneous sources such as metabolic, transcriptomic, transcription factor and protein-protein interaction data have been published [1]. Concurrently, several studies have revealed the dynamics of small-scale molecular networks using innovative applications of reverse-engineering techniques, modeling and computer simulation of kinetic equations [2]. However, observations regarding overall parameter robustness of biological circuits may imply that beyond the apparent model complexity there exists, in a mathematical sense, one or several core systems or operational principles driving the system dynamics that are shielded by the more complex original dynamical equations [3][4][5][6]. Such a dynamical core is not necessarily apparent from the connectivity graph or from the original dynamical equations. The existence of such a core has been supported by studies employing an ensemble simulation approach in the analysis of various computational models [7,8]. Essentially, a large number of different parameter configurations exist that produce the "same" dynamical output from the system, as emphasized in a study by Gutenkunst et al [9], whereby, using a large set of biochemical models, they found a "sloppy" parameter sensitivity spectra to be a universal characteristic of these models.
To advance our understanding of molecular networks, we need to develop methods that permit an analysis of such underlying core dynamics. Research on model reduction and reformulation has been at the forefront within the systems biology community over the last few years (e.g. [4,[10][11][12][13][14][15][16]). Several studies have targeted the elaborate original model using different approaches, and have arrived at a more compact system description, which in different ways capture some aspects of the essence of the original model dynamics. Using model reduction tools, the number of degrees of freedom are reduced, while still remaining within the same model formalism as the original model. By reformulating the model, the original model is approximated to another hopefully more transparent modeling formalism, with a structure that better captures the key qualities of the original system. Examples of model reduction [17] include lumping of variables [10,11], separation of timescales [14] (or, for example, the classical Michaelis-Menten equation describing enzyme kinetics), sensitivity analysis based methods, and methods based on identifiability analysis [13]. Examples of transformation of modeling formalism include boolean approximations [4,12], hybrid stochastic approximations [18], or the simplified model used throughout this study, a delayed piecewise linear approximation [15] of an ordinary differential model.
The simplified model has subsequently been used to characterize parameters of the original model, as in Radulescu et al [14], where critical parameters are found through reduction to dominant subsystems, and after this analytically mapped back to the original model. Here, in contrast, we started from the original model, by exploring the dynamics with sensitivity analysis, and next compared this to the predictions from the simplified model. The aim was to investigate to what extent it is feasible, by using an underlying core description [15], to explain dynamical properties of the original model ( Figure 1).
The cell cycle is one of the most extensively studied and essential biological regulatory circuits. The differential computational models pioneered and developed by Novak and Tyson [2,19] describing cell cycle regulation of fission yeast, incorporates and accounts for a large body of experimental data. During the cell cycle, the cell grows, DNA replicates, and the cell divides into two daughter cells. To maintain cell size over several generations, it is essential that there is a controlled balance between cell growth and cell division. The model used in the present study [19] describes the key components of that regulatory mechanism, consisting of proteins and protein complexes that drive the cell cycle in fission yeast. It is detailed in the Methods (equations (3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)) and is referred to throughout as the NT-model. The NT-model is thus our point of departure for the model simplification process and sensitivity analyses. For a detailed description of the NT-model see [19]. For a more general discussion of cell cycle dynamics see [2].
To examine the dynamical properties of the original NT-model, we performed a sensitivity analysis. Two different types of output responses to parameter perturbations were considered. One was qualitative and described whether the system was still working, i.e. whether there were oscillations (cell divisions) corresponding to parameters that here we call "essential". The other response pattern was quantitative and measured changes in cell size (cell growth), which described "modulatory" effects of parameter perturbations. This is primarily a methodological distinction, i.e. response effects related to parameter perturbations. In principle, parameters could be both essential and modulatory, depending on the perturbations in question. In practice, when conducting sensitivity analyses, the key question is which parameters affect the system output, how, and under what conditions. Here, we extend the scope of the sensitivity analysis to evaluate the correspondence between perturbation effects in the original model and analytical results of the simplified model, focusing on two types of possible perturbation effects -essential and modulatory changes -anticipating that the simplified model should be able to account for both of these situations. The analysis includes three key steps: (1) perturbations of the parameters of the original NT-model; (2) a corresponding recalibration of the parameters of the simplified model; and (3) systematic comparison of the perturbation effects on the orginal and simplified models ( Figure 1). The perturbation effects of the original model are investigated by numerical simulations whereas the corresponding effects of the simplified model are retrieved by evaluation of analytical expressions.
In [15] we introduced a scheme by which we could simplify the NT-model using step functions and a time delay. Introducing discrete step-functions resulted in a hybrid model containing both continuous and discrete variables. In the case of the NT-model, this hybrid model turned out to be piecewise linear. This delayed piecewise linear (DPL)-model can be found in the Methods (equations (17)(18)(19)(20)), and is referred to throughout as the hybrid, or the DPL-model. A piecewise linear description facilitates system analysis, as each individual linear system can be analyzed separately using well established linear system analysis tools. The combined results from these separate analyses exposes the behaviour of the full system ( Figure 2).
Sensitivity analysis refers to a broad range of mathematical and statistical analysis methods that serve to evaluate the relative merits of a set of parameters in controlling system dynamics. The aim is to relate variation (uncertainty) in the system output to variation (uncertainty) in the system input. The analysis method chosen depends on the goal. Campolongo et al [20] distinguish between three types of analysis: (1) factor screening; (2) local sensitivity analysis; and (3) global analysis. The first type refers to an analysis targeting the effects of parameter variation one at a time, i.e. variation in system output is related to parameters one at a time. It is a more or less qualitative assessment of the relevant factors in system dynamics. By contrast, local and global sensitivity analyses imply detailed quantitative assessments of relative effects of the isolated factors of concern. Local analysis makes use of differential calculus, and stability analysis of equilibrium points. Global sensitivity analysis targets response effects in the whole parameter space, as well as interaction effects between parameters. It relies to a larger extent on statistical methods. Our purpose in this study was to clarify the relevance of parameters in cell cycle dynamics. As a starting point, we performed a screening study to investigate the relative merits of single parameters in cell cycle regulation.

Results
Our analysis is based on a set, P org , of parameter perturbations of the original NT-model, consisting of 470 single parameter changes. This set was used to analyse the relationship between parameter values and model output (cell size and cycle time) through two different tracks, see Figure 1. During the first track, the sensitivity analysis, numerical simulations of the NT-model were conducted for all perturbations in P org , and scores of sensitivity were constructed for each model parameter. During track two, all perturbations of P org were translated to corresponding perturbations of the hybrid model and, instead of simulations, mathematical conditions were evaluated in order to predict the sensitivity to the perturbations, and thereby construct scores of sensitivity. Two measures of sensitivity were used in both tracks, targeting essential and modulatory effects of perturbations, respectively. Finally, the sensitivity scores from the two tracks were compared in order to i) validate the predictions from the hybrid model, and ii) investigate the role of essential and modulatory parameters and perturbations. The results are divided into the following sections: (I) sensitivity analysis of original model; (II) analytical predictions from the hybrid model; and (III) comparison between the analytical predictions and the sensitivity analysis. Finally, (IV) the hybrid model is used for more complicated predictions.

I Sensitivity analysis of the original model identifying essential and modulatory perturbation effects
The parameters of the original model, k1,..., μ = p org (equation (16)), were perturbed one at a time, starting from a default set of parameter values corresponding to the wild-type cell (equation (16)). A perturbation consists of a change of the parameter value by a factor ps j (the relative Perturbation Size), where ps j ps = (10 -1 , 10 -0.8 , 10 -0.6 , 10 -0.4 , 10 -0.2 , 10 0 , 10 0.  simulation, the output trajectory was examined for i) whether there were consistent oscillations or not, ii) cycle time t CT and iii) final (before cell division) cell mass M end .
i) Out of a total of 470 perturbations (47 × 10, not counting ps j = 1), 437 corresponded to trajectories with consistent oscillations. ii) When inspecting these trajectories, two distinct patterns emerged -either the cycle time t CT was constant within the trajectory, or, for a few of the perturbations, t CT alternated in a repeated pattern, e g. t CT = ..., 111.0, 166.25, 111.0, 166.25, 111.0,... (denoted quantized cell cycles using the formulation of Novak and Tyson [19]). Surprisingly, for perturbations with constant cycle times, most had the same cycle time t CT = 138.6 (in total 423 perturbations), giving a very narrow distribution of cycle times (see Figure 3). One exception to this was the parameter μ, which determines the speed of cell growth, and naturally has a large impact on cycle time. iii) For final cell mass M end , the pattern reversed. The parameter perturbations generated a broader distribution of M end ( Figure 3).
We next defined measures of sensitivity summarizing how sensitive the system was to perturbations of the single parameter p org i . Two measures were defined, one based on the cycle time of the oscillations, denoted S osc i , and the other based on final cell mass, S mass i . Since the cycle time was, in almost all cases, either constant t CT = 138.6 or non-existing (quantized cycle times or no oscillations), we decided to consider only whether we had proper cell cycle oscillations or not in the definition of S osc i . Proper oscillations were defined as trajectories that still oscillated after 4000 min, and had a narrow distribution of cycle times (to remove quantized behaviour). We therefore ended up with one discrete-qualitative measure S osc i and one continuous-quantitative S mass i . In the calculation of S mass i , only simulations with proper oscillations were considered. We denoted perturbations with a large effect on S osc essential, and with a large effect on S mass i modulatory.
The sensitivity analysis was carried out using two different ranges of parameter perturbations, narrow range perturbations close to the default parameter value, and wide range perturbations both close, but also further away. For narrow range perturbations, S osc i and S mass i were calculated using only the two perturbations closest to the default value, i.e. ps narrow = (10 -0.2 , 10 0 , 10 0.2 ), while for wide range, all perturbation sizes were used, ps wide = ps.

Essential parameters
The sensitivity measure, S osc i , corresponds to the number of perturbed simulations of parameter p org i that ended in cell collapse, i.e. a lack of proper oscillation, S osc i = j s ij, where s ij = 0 if perturbation j resulted in proper oscillations, and s ij = 1 otherwise, j 1,..., n pert , where n pert = 11 for wide and n pert = 3 for narrow range perturbations. For narrow range perturbations, all parameters have S osc i = 0 (Table 1 original model), and therefore the system is robust against perturbations of this size when proper oscillations are considered. We here define robustness as persistent system dynamics after perturbation. It is quantified with the sensitivity scores (being the opposite to sensitivity). Even for larger perturbations, oscillation is a robust behaviour for most parameter perturbations (Table 1 original model, wide range perturbations). Still, some perturbations do affect oscillations, e.g. the rate constants k 9 and k 10 of the hypothesized intermediate enzyme IEP, of the Slp1/APC subsystem, which remove the protein complex MPF at the end of mitosis. In the coming analysis, we denote parameters with S osc i > 0 as essential, and with S osc i = 0 as non-essential.

Modulatory parameters
The sensitivity measure, S mass i , corresponds to the average effect on final cell mass M end of the n perturbations of parameter p org i that had proper oscillations (s ij = 0) ( Figure 4). It is defined as S mass is the final cell mass after the j'th perturbation of parameter p org i , and M end def is the final cell mass using the default parameter set (16).
With narrow range perturbations (Table 2 original model) the parameters with the highest sensitivity (S mass i ) are related to either the proteins Cdc25 or Wee1, or to the formation or destruction of the protein Cdc13. This is in agreement with the fact that, in the original NT-model, final cell mass is in large determined by the SNIC bifurcation [19] corresponding to the G2/M checkpoint. This checkpoint is realized by the phosphorylation and dephosphorylation of Cdc25 and Wee1.
With wide range perturbations (Table 2 original model), the top scoring parameters also include k 16 and k 15 , which are related to TF (a transcription factor) activity, and k 4 related to the activity of the protein Ste9. This is most likely a result of the G1/S checkpoint not

IIA Analytical predictions from the hybrid model on essential and modulatory perturbation effects
In [15] we introduced a scheme by which we could simplify the NT-model. Dynamical switching modules were identified and replaced by time lagged step-functions, resulting in the hybrid DPL-model. This model approximation process is summarized in the Methods. The hybrid model (equations (17)(18)(19)) consists of four linear systems, defined by the matrices A kl , k, l {1, 2}, where the linear system at each time-point t is determined by, so called, switching rules (equation (18)). Figure 2 displays a numerical simulation of the DPL-system, indicating the different linear systems. A switch between systems is initiated when the output y(t) = [MPF](t), corresponding to the concentration of the protein complex MPF, crosses a threshold θ 25/wee or θ slp/ste . Note that when θ 25/wee is passed, the switch is immediate, whereas, when θ slp/ste is passed, the switch is delayed with τ minutes. The switching rule can therefore be considered as equipped with a memory, and we have denoted the model delayed piecewise linear. An analysis of the DPL-model [15] showed that each separate linear system, defined by the system matrix A kl , was stable, having an asymptotically stable fixed point (x kl = −A −1 kl Bu ext , with real negative eigenvalues). Combining the four linear systems with the switching rules revealed, however, that some of the fixed points could never be reached without passing a switching threshold (and thereby a change of linear systems). Thus, these  fixed points were no longer stable. The relationship between the location of the fixed points and the switching thresholds, therefore, determines the dynamical behaviour of the system. This relationship between fixed point and switching thresholds corresponds to analytical expressions including the parameters of the hybrid model, and allows us to present two types of important parameter/output relationships. The first type refers to necessary constraints (NC), which must be satisfied in order for the model to generate proper oscillations. Therefore NC identifies essential parameter changes. The second type is a function describing the size of the cell, or Size Function (SF). SF describes how cell size is influenced by hybrid model parameters, and therefore SF identifies modulatory parameter changes and their effects on cell size. When using these expressions in order to find out whether a parameter change is modulatory, essential or both, we always have the default parameter set as the point of departure for the change.
With another set of default parameter values, the result could of course change. Note that NC is a qualitative test -either a parameter change violates the conditions (and the change is considered to be essential), or it does not; whereas SF is quantitative, giving numerical predictions on cell size.
Analytical results identifying essential parameter changesnecessary constraints NC Two constraints, that are necessary in order to have cell cycle oscillations, were retrieved in the analysis of the hybrid model [15], and are further extended in the Methods section (equations (24,25)). These correspond to inequalities, including hybrid model parameters, that define regions of parameter values for which there cannot be proper oscillations (note that the necessary conditions only state that if not satisfied, there will not be oscillations, but not the opposite, i.e. they do not assure oscillations, which agrees with the definition of "essential" in the sensitivity analysis, i.e. anhilation of oscillation due to parameter perturbations). A change of hybrid parameter τ...μ = p hyb (equation (20)), that violates these constraints is considered as essential.
Analytical results identifying modulatory parameter changes -the size function SF The second type of parameter/output relationship, the function SF (equations (21,22), Methods), describes how the final size of the cell (i.e. just before cell division), is influenced by hybrid model parameters. SF corresponds to the formation of a concentration threshold. When the concentration of the protein complex MPF passes this threshold (under normal circumstances), cell division is initiated. In cell cycle physiology, it corresponds to the G2/M checkpoint [2].

IIB Translating the analytical results from the hybrid model to the original model
To be able to compare the hybrid and original models (i.e. compare their respective parameter characterizations) we needed to translate the analytical result of the hybrid model to the original model. For this translation, we first defined a mapping from original model parameters to hybrid model parameters. We then used this mapping to backtrack the analytical results of the hybrid model to the original model parameters.

Mapping original model parameters onto hybrid model parameters
There is no simple algebraic relationship that could be used to define mappings between original and hybrid model parameters parameters, p org i ∈ p org , i = 1...n org and p hyb m ∈ p hyb , m = 1...n hyb , i.e. there is no analytical function p hyb m = f (p org ) for all m that can be used to relate the original and hybrid models. Instead, we constructed mappings numerically for different specific values of p org , namely those values corresponding to the perturbations used in the sensitivity analysis. First, each original parameter p org i was perturbed n pert = 11 times according to the relative perturbation sizes, ps, described earlier, i. e. the parameter was given the value p org i = p org ij , j = 1... n pert (note here that single index p org i denotes a parameter, whereas a double index p org ij denotes a specific parameter value). Next, the effect from the perturbation was translated to the hybrid parameters p hyb . The translation was done using a recalibration process, adapting the hybrid model to the original perturbed one through re-performing the hybrid approximation process and measuring the effect of the perturbation p

Backtracking the analytical result from the hybrid model to the original model
After recalibration of the hybrid model, we translated the analytical result in terms of modulatory and essential effects, by backtracking the analytical conditions (NC and SF) to the original model. This means that for each perturbation of an original model parameter p org ij , we constructed a new hybrid model and examined the analytical expressions determining essential and modulatory effects. From this, we finally characterized the original parameter in terms of essential and/or modulatory. To sum up, the essential and modulatory effects of the original parameter perturbations are measured in two main ways: with respect to the original model (track 1), and with respect to analytical constraints and implications after recalibration of the hybrid model (track 2). Bactracking the essential conditions To backtrack the necessary constraints (NC, corresponding to essential parameter changes), to original parameter characterizations, we counted the total number (out of n pert ) of perturbations of p org i which have a mapped effect such that the necessary constraints (equations (24,25)) are violated. From this we calculated a sensitivity score similar to the one from the sensitivity analysis, but based on the hybrid model, and therefore denoted hyb S osc i . Thus There are two hybrid parameters that are not considered in the correspondence analysis. The first one is μ, which is related to the growth of the cell, i.e. the increase in cell mass M. Since the hybrid model is based on the assumption that the increase in M is very slow compared to other variables, perturbations of μ are not considered. The other parameter is τ, which has not been included in the analysis for practical reasons (Methods). Not all original parameter perturbations have a mapping to a hybrid parameter (data not shown). This is partly a reflection of the fact that parts of the original system are inactive during the experimental context that the hybrid model seeks to reproduce, and has therefore not been included in the hybrid model.

III Correspondence analysis of hybrid and original parameter characterizations
Once we had conducted the sensitivity analysis of the original model (I), and retrieved the corresponding predictions from the hybrid model (II), we were now in a position to compare the results (Figure 1). The parameters were divided into essential or non-essential (Table 1) and modulatory or non-modulatory (Table 2), with respect to the model in question: original or hybrid. We considered both narrow range and wide range perturbations. A parameter was characterized as modulatory if S mass i > 0.1 (or hyb S mass i > 0.1) and as essential if S osc i > 0 (or hyb S osc i > 0), and the parameter characterizations of the original and hybrid models were compared.
No parameter was detected as being essential by narrow range perturbations, neither from the original nor hybrid model (Table 1). Wide range perturbations were needed to get an essential effect. For modulatory characterizations, the correspondence was better for the models at narrow range perturbations than at long range perturbations ( Table 2). The result is summarized by a binary (essential/non-essential or modulatory/non-modulatory) classification test (Table 3). We view the parameter characterizations from the hybrid model as a predictor of the class of the original model parameters, and calculate the sensitivity and specificity of the predictions for each class (Table 3). For the narrow range modulatory perturbations, the prediction is close to perfect, with a sensitivity of 100% and a specificity of 95% (Table 3). For narrow range essential predictions the specificity is 100% and no hybrid or original model parameters are found to be essential. When wide range perturbations are considered, there are more misclassifications. For modulatory parameter predictions the sensitivity decreased to 60%, while the specificity was 96% and for essential predictions, the sensitivity was 67% and the specificity remained 100% (Table 3).
In the hybrid model, modulatory parameters are defined by the size function SF, which corresponds, as earlier described, to the G2/M checkpoint of the cell cycle. The proteins behind this checkpoint are those that are most important to regulate the size of the cell in this context. If we did not know this and were to identify the mechanism regulating cell size, from this model comparison, we can see that we would be better of by the use of narrow perturbations rather than wide ones. If wide range perturbations are used, other parameters, which are not involved in this checkpoint, are also incorrectly identified, and therefore confusing the picture. We can conclude that the correspondence between models is better near the default values of the model parameters (narrow range), and becomes worse further from this point. At the same time, it seems impossible to identify essential parameters without using perturbations further away from the default parameter setting. Finally, if we use narrow range perturbations to identify modulatory parameters and wide range perturbations to identify essential parameters, then we note that most parameters are either essential or modulatory for this set of perturbations.

Mis-predictions
There are two main types of mis-predictions at wide range perturbations. First, parts of the original model, for which there is no counterpart in the hybrid model (parameters indicated with * in Table 1 and Table 2), and that are not active at narrow range perturbations, have in fact an essential or modulatory effect at wide range perturbations. This indicates that the hybrid model performs quite well for predictions locally around the point in parameter space where the hybrid model was constructed. Further away from this point, the predictions are still quite good for those parts of the original system that were included in the hybrid model. The hybrid model cannot, however, capture what will happen to parts that were not included. Second, there are some mis-predictions that are most likely due to the approximation and discrepancy between the hybrid and original models. One such prediction was the parameters k 9 and k 10 , which are found to be more essential in the original model than in the hybrid model. k 9 and k 10 are important actors to obtain quantized output behaviour. Such behaviour cannot be mimicked by the hybrid model; rather, the hybrid model predicts that the system will have proper oscillations for some of the perturbations, which gives quantized oscillations. Also, the fact that we have not included the hybrid model parameter τ in the analysis might affect predictions.
Finally, we have a problem with mappings to the hybrid model parameter θ 25/wee . This parameter is different compared to the other hybrid parameters since it corresponds to the switching threshold of two step functions, s 25 and s wee (Methods), when their respective switching thresholds are at the same level i.e. θ 25/wee = θ 25 = θ wee when θ 25 = θ wee . To be able to map a perturbation of the original model to θ 25/wee p hyb that keeps this relationship (θ 25 = θ wee ), and thus the structure of the hybrid model, at least two original parameters have to be perturbed simultaneously. However, the single parameter perturbations p org ij used in this study, only correspond to a change in either θ 25 or θ wee . To obtain a numerical value for each mapping p org ij → θ 25/wee despite this, we assume that θ 25/wee = min(θ 25 ,θ wee ), when the threshold is passed from below, and that θ 25/wee = max (θ 25 , θ wee ), when the threshold is passed from above. Presuming that as soon as either θ 25 or θ wee is passed, the other threshold will also soon be passed.

IV Combinatorial predictions from the hybrid model
After investigating the correspondance between the original and hybrid model concerning single parameter perturbations of the original model, we used the hybrid model to identify combinations of parameters aimed at a specific critical effect.

Essential combinations of parameters
We used the necessary constraints (24,25) to find hybrid model parameters with an essential effect on model output. To evaluate the constraints, we used, as earlier, a set of perturbations, but this time investigating hybrid model parameters directly instead of translating the result to original model parameters. From these perturbations (Table 4) we could see that (at least) hybrid model parameters θ 25/wee , h wee and h slp/ste can have an essential effect.
The hybrid parameter θ 25/wee was shown to have the largest effect (Table 4), using these perturbations. Simulataneously, none of the original parameters mapped to θ 25 or θ wee were found to be essential (data not shown). This emphasizes that there are important (essential) features in the original model that are robust towards single parameter perturbations. We described earlier that θ 25/ wee p hyb can only be mapped to from combinations of original parameters if the structure of the hybrid model is to be retained (i.e. θ 25 = θ wee ). We therefore expect a  critical effect on the original system by changing such a combination. A pairwise mapping to θ 25/wee p hyb can be achieved by simultaneous changes in V i25 p org and V awee p org or V a25 p org and V iwee p org . By adding the parameter s p org such that V i25 sV i25 , V awee sV awee , and the parameter t p org such that V a25 tV a25 and V iwee tV iwee , then a perturbation of s p org or t p org in the original model corresponds to changing θ 25/wee p hyb . Using the same perturbation scheme as before (wide range) we obtained the following sensitivity scores for the combined perturbations S osc s = 3 and S osc t = 3, which can be compared to the score of S osc = 0 for V i25 , V awee , V a25 and V iwee , respectively. This shows that by combining perturbations of two, by themselves, robust parameters, we can obtain a critical effect.

Modulatory combinations of parameters
From the hybrid model we retrieved an expression for the critical cell mass, M C , above which mitosis is initiated [15] M C = θ 25/wee (k 2 + l slp/ste )(k 2 + l 25 + h wee + l slp/ste ) where all constants are hybrid model parameters. Mathematically, this corresponds to a bifurcation point at which the system moves from a stable fixed point to a limit cycle. Under circumstances corresponding to a normal cell, this expression can be used to approximate final cell mass before cell division.
If we assume that l slp/ste << k 2 , l 25 <<h wee , as for the normal cell (equation (20)), then equation (1) can be approximated with Using this expression, we can compute quantitative predictions of cell size, e.g. if k 1 is reduced by 50%, the cell size is expected to be twice as large. This was tested in the original model by setting k 1 = 0.015 and running a numerical simulation. Final cell mass increased from 2.0 to 4.0. If, together with halving k 1 we also double θ 25/wee and h wee , cell size should increase by about 2 × 2 × 2 = 8 times. Testing this in the original model by setting k 1 = 0.015, V i25 = V awee = 0.5 (corresponding to doubling the hybrid model parameter θ 25/wee ), k wee = 2.6 (doubling h wee ), final cell size increased from 2.0 to 14.4. We can therefore predict multiplicative effects by using the hybrid model.

Summary of results
In this case study of the cell cycle, by comparing the sensitivity analysis of the original model with analytical predictions from the hybrid model, we found that i) predictions for essential and modulatory parameters from the hybrid model were satisfactory for small perturbations, but more discrepancies were seen for larger perturbations; ii) small perturbations were of more use when trying to identify important modulatory mechanisms compared to larger ones whereas, iii) larger perturbations were needed in order to identify essential parameters. Furthermore, iv) most parameters were either essential or modulatory using this set of perturbations. Finally, v) the most sensitive feature of the system (the threshold corresponding to θ 25/wee ) was robust towards single parameter perturbations.

Discussion and Conclusions
Any biological process is essentially an interaction between different elements in time and space. Given experimental data and prior knowledge we may be able to formulate such interactions in mathematical terms. As a rule, such models suffer from uncertainty in the structure, both from a lack of representation of essential elements and interactions, as well as from erroneous model representation of the phenomena under investigation [21,22]. Moreover, parameters are uncertain in the sense that several different combinations are consistent with a given set of data. When dealing with a specific model, such as the NT-model for the cell-cycle [19], we are in a unique position because a large body of experimental work has already been encoded in mathematical models by the pioneering work of Novak and Tyson [23]. Furthermore, this specific NT-model has been analyzed using both bifurcation [19] and mathematical piecewise linear approximation techniques [15]. Bifurcation analysis indicates the dynamical mechanisms behind the transitions within the cell-cycle, while the mathematical simplification analysis provides proof of the existence of a limit cycle, yields explicit analytical expressions describing how model parameters affect cell size, and also provides conditions necessary for oscillations. However, it is generally very difficult to succeed with a thorough mathematical analysis of a given computational model, and this becomes even harder for larger models. Numerical bifurcation analysis is also a technique that is most efficient for smaller models. For larger models, we are left performing numerical simulations to explore the sensitivity of the model, and thereby increase our understanding of what parts of the model may be important for different aspects of the system's behaviour [24]. For these reasons, we have here interpolated between the results obtained from a sensitivity analysis of the original model, and the outcome from a previous model simplification [15], with the rational being that such an analysis may be instructive for assessing what we can and cannot discover using a sensitivity analysis. Such understanding may prove useful for those cases where an explicit mathematical simplification is not feasible.
Our analysis of the simplified model resulted in the formulation of two types of mathematical constraintsone identifying parameter perturbations with a qualitative effect on model output (NC), and the other describing perturbations with a quantitative effect (SF). From the corresponding sensitivity analysis of the original NT-model, we could also here identify perturbations that had a qualitative effect (removing oscillations), or quantitative effect (changing cell size). Using the same set of parameter perturbations and sensitivity measures in both the hybrid and original models, we could compare our analytical predictions to the sensitivity analysis. Based on our definition on modulatory and essential parameters, both methods characterized the parameters more or less equally for narrow range perturbations; for wide range perturbations there was a discrepancy.
This combined analysis demonstrated that the parameters were, in most cases, characterized as either essential or modulatory using this set of perturbations, where essential parameters were only identified using wide perturbations. This means that there are parameters that do not have any effect at small perturbations, not even modulatory, but which do have a critical (essential) effect at large perturbations. They are therefore very important to the system even though this would not be obvious from a local sensitivity analysis. Modulatory parameters, as defined by the size function (SF) of the hybrid model, were best identified by small perturbations. When larger perturbations were applied, parameters from other parts of the system also showed up, confusing the picture. SF describes the parameter relationship of the G2/M checkpoint -the proteins behind this checkpoint are those that are most important in cell size regulation. If we did not know this and were to try to identify the mechanism regulating cell size, we would be better to use small perturbations. Essential parameters, on the other hand, could not be found using small perturbations alone; rather larger perturbations had to be used. This is not surprising considering that the biological system has to be robust in the default parameter setting, and the corresponding point in parameter space must therefore be located far away from the region where there are, for example, no oscillations. Using the hybrid model for more complex predictions and comparing these to the sensitivity analysis, we found that the threshold corresponding to the hybrid parameter θ 25/wee , which has an essential effect (and modulatory) in the hybrid model, were robust (using S osc ) towards single parameter perturbations; however it had an essential effect with a targeted combined perturbation. We also demonstrated the possibility to find combinations of parameters with a multiplicative effect on cell size.
In the case of single parameter perturbations, it is possible to numerically simulate the behaviour of a large model to all possible parameter perturbations. With perturbations of more than one parameter at a time, the combinatorial explosion has the effect that all parameter combinations cannot be tested. A simplified model therefore has the potential to inform us about which parameter combinations that are important to test. Our study pin-points the difference between parameters that have a quantitative, modulatory impact on model output, where the perturbation size is correlated with the change in model output, and parameters that are essential to the system in a qualitative manner.
Our investigation has several similarities to a recent study of Radulescu et al [14] where they suggest the use of a hierarchy of reduced models to describe complex biological networks. Radulescu also map critical parameters identified from the reduced model, back to the original model. However, there are also marked differences. In Radulescu et al, model reductions are by and large faithful in both structure and components to the original model, letting kinetics guide structural decomposition and parameter elimination, keeping the most important components and not introducing new ones. In the process, critical monomials and thereby critical parameters are identified. This is not the approach taken here. There are no simple and direct analytical relationships between the original and simplified models. Instead we look at different perturbations in the original model and recalibrate the simplified model from those in order to backtrack findings in the simplified model. By using a new model formalism, in this case, a piecewise linear with delay, we can point to other important features of the original model. In this case that the system seems to switch between different linear systems and that the location of the fixed points and the thresholds are most important for the dynamics. We also find that a very sensitive feature of the original model, the threshold described earlier, is robust towards single parameter perturbations.
A simplified model is an approximation of a more complicated system, and there will be instances when it does not work. Its predictions must therefore be confirmed in the full model context. However, a discrepancy between the simplified and original model is of interest in itself, since it can point out an approximation that has failed. Since approximations are a way to learn about the essentials of a system, this gives information on neglected essential features. To combine model simplification with sensitivity analysis is a way to investigate the structure of complex biological networks. We believe that this approach is applicable to systems other than the cell cycle, even though the precise set up used in this study may not be suitable. Sensitivity analysis has to be carefully performed, since it is almost impossible to investigate the full parameter space for any model with more than a few parameters. Here we performed a first step screening study, modifying one parameter at a time. The next logical step would be to consider a more global sensitivity analysis by investigating combinations of parameters, which was beyond the scope of this study. In this comparison between the hybrid model and the original model, we have considered all simulations with a quantized behaviour as not working properly, and therefore disqualified them from the analysis. This is partly because a normal cell does not display this type of behaviour, but also because the hybrid model does not mimic the behaviour of the original model in a few of these quantized instances.

Methods
The original NT-model [TF] = G(k 15 The default parmater values corresponding to the wild-type cell are The hybrid, DPL-model The hybrid model approximation process is described in [15] and summarized further down. Let x(t) = (x Cdc13 T (t) x PreMPF (t)) represent the state of the cell cycle system ( x Cdc13 T corresponds to NT-model variable [Cdc13 T ] and x PreMPF to NT-variable [preMPF]). Further, let u ext (t) = M (t) (cell mass) be the external input, and y(t) the output from the system (DPL-model variable y corresponds to NT-variable [MPF] and we sometimes, a bit sloppy, write y = [MPF]). Then, the DPL-model can be written as slp/ste}, except for θ 25 and θ wee . During the model simplification process we noted that θ SM 25 = θ SM wee (when using default parameters) and this was used as a further simplification step by merging θ 25 and θ wee into one parameter θ 25/wee . Therefore, only perturbations of the original model which keep the θ SM 25 = θ SM wee property are consistent with the hybrid model. These are all pairwise perturbations.