Simple network model
To illustrate the shortcoming of local PSA in analyzing system dynamics, consider a simple six state model involving three reactions with Michealis-Menten (MM) kinetics, as shown in Figure 1(a) (model parameters, rate kinetics and initial concentrations are given in Additional File 1: Supplementary Table S1). In this network, the activation of x6 followed a switch-like dynamics in response to the stimulus x1, as illustrated in Figure 1(b) (nominal). The model describes two pathways that contribute to x6 activation: a direct x2 pathway and an indirect x2, x3, and x5 pathway.
In this example, in silico knock-out experiments were performed by removing each pathway individually in order to assess the dominance of one pathway over the other in x6 activation. Both full network and knock-out (KO) simulations were performed under a stimulus of x1(t0 = 0) = 1. As illustrated in Figure 1(b), while the initial x6 activation in the indirect pathway KO remained the same as that of the original model, the switch-like activation was much less pronounced. On the other hand, the original switching behaviour was preserved in the direct pathway KO, but the switching time was delayed due to a slower initial activation. Taken together, these KO simulations suggested that the x6 activation is mainly accomplished through the indirect pathway, while the direct pathway contributes mainly to the initial x6 activation.
Parametric sensitivity analysis for dynamical systems: A caveat
Local parametric sensitivity analysis was also used to study the pathway dominance in this simple network. Mathematically, the parametric sensitivity coefficient is given as
where x
i
is the i-th state in an ODE model and p
j
is the j-th kinetic parameter of an ODE model (for a detailed description of sensitivity coefficient derivation, see Methods). These sensitivity coefficients describe the change in system output (state trajectory) at time t with respect to (an infinitesimal) perturbation on the system parameter values at time τ. Here, the PSA was performed for the same stimulus x1(0) = 1 with τ = 0 and the sensitivity coefficients were computed for the time range of 0-15 time units. The term local sensitivity analysis refers to the fact that the results will depend on the nominal parameter values around which the derivatives in Eq. (1) are calculated.
The sensitivities of x6 with respect to all model parameters are ranked in Figure 2 using consolidated sensitivity metrics, i.e. infinite norm [44] (Figure 2(a)), FIM [43] (Figure 2(b)), and time-integral [34] (Figure 2(c)), and using sensitivity magnitudes at switching time (t = 7.12 time units; Figure 2(d)). The conclusion from these rankings was the same: (1) the largest sensitivity was associated with the kinetics of x1 conversion to x2 and (2) the direct pathway (r2) parameters have larger sensitivities than those from the indirect pathway (r4), suggesting larger influence of the direct pathway on the x6 activation. Hence, the conclusion from the PSA is in direct contradiction with the findings from in silico KO experiments.
This discrepancy can be explained by taking a closer look at the way parametric sensitivity coefficients in eqn. (1) are calculated:
where τ = t0 is the usual perturbation time,
is the time derivative of sensitivity coefficient
(see Methods for detail) and H(t) is the Heaviside step function. In this case, since model parameters are static or time-invariant, the parametric perturbations in the PSA consist of step changes in the parameter values, as depicted in Figure 3(a). Hence, the sensitivity coefficients at time t represent an integrated or accumulated change in the states from τ to t due to a persistent parameter change started at time τ (see Figure 3(b)). Indeed, substituting the full equation of
(see Methods) in eqn. (2) gives
Here, one sees two terms in the integrand that contribute to the sensitivity coefficients at time t: (1) the first is related to the (integrated) sensitivities that are carried over from the initial perturbation time τ and (2) the second accounts for the instantaneous rate changes due to the parametric perturbations that still persist at time
Thus, in the PSA, a large sensitivity magnitude of Si,j (t,τ) indicates the importance of the j-th parameter in time window of τ and t, during which the perturbation is applied to the system. Hence, the use of these coefficients to infer the dynamical importance of parameters is inappropriate and can even be misleading.
For the reason above, the PSA of the simple network model gave an incorrect conclusion regarding direct versus indirect pathway activating x6. As seen in the in silico KO experiments, the direct pathway regulates the initial activation of x6, while the actual switching is carried out by the indirect pathway. In the PSA of this model (τ = 0), the early importance of the direct pathway and also the reaction r1 persisted beyond the initial times in the sensitivity coefficients due to the aforementioned integrated effect. In this case, the importance of the indirect pathway was not apparent from the parameter sensitivity rankings in the background of large (integrated) sensitivities with respect to r1 and the direct pathway. Correspondingly, the use of any time-consolidated sensitivity metrics will only worsen this problem.
Impulse parametric sensitivity analysis (iPSA)
As the problem with local PSA mentioned above is rooted from the persistent parameter perturbations, which is also done in global PSA [45], a new sensitivity analysis is formulated here that introduces impulse perturbations to model parameters as shown in Figure 3(c-d). The corresponding impulse sensitivity coefficients iSi,j (t,τ) reflect the change in the i-th state x
i
at time t due to an impulse perturbation on the j-th parameter p
j
at time τ (see Methods for the derivation and definition). Since impulse perturbations on parameters cause an immediate state changes at the perturbation time (see Methods), the inference of dynamical parametric importance can be obtained from impulse sensitivities by varying the time of perturbation. However, as with the local PSA, impulse perturbations are also local in nature and thus the impulse sensitivities will depend on the nominal parameter values. The global equivalent of iPSA can be formulated using pulse perturbations and is currently under development. Finally, from time-varying impulse perturbations, the iPSA can give answer to the following questions about state dynamics: which are the important parameters and when do they become important?
The iPSA was also applied to the simple network model under the same stimulus and for the same time range as above. Since the impulse perturbations are delivered at different times, the impulse sensitivity iSi,j (t,τ) has two-time dependence, with respect to the time of perturbation (τ) and the time of observation (t). Figure 4 shows the iPSA sensitivity coefficients of x6 at the end of simulation time (t = 15 time units), with respect to the four most important parameters at different perturbation times (for complete iPSA sensitivities, see Additional File 1: Supplementary Figure S1). In agreement with the KO simulations and in contrast to the findings from PSA, the impulse sensitivities gave support to the dominance of the indirect pathway. Specifically, the results showed that x6 activation: (1) is initiated by r1 (high initial sensitivity to parameter kf1); (2) is accomplished mainly by the direct pathway prior to the switching time (sensitivity to r
2
is higher than to r4 during these times); and (3) is subsequently carried by the indirect pathway (highest sensitivity to r4 during switching times). A higher resolution analysis using a heat map of the complete iPSA coefficients with t and τ between 0 and 15 time units gave the same conclusion (see Additional file 1: Supplementary Figure S2).
Fas-induced cell death model of human Jurkat T-cells
In the second case study, we consider a more complex biological network taken from the modelling of programmed cell death in Jurkat cancer T-cells (see Figure 5). The model equations and parameters were identified from experimental data [37] (see Additional file 1: Supplementary Table S2 for detailed reaction rates and parameters). In this network, the cell death is decided by the cleaving of procaspase-3 into caspase-3 [46], in response to FasL stimulus. The model simulation showed that caspase-3 is switched ON at around t = 6000 seconds (see Figure 5, inset) and like the simple network above, there exist two activating pathways: the direct caspase-8 (type-I) and the indirect mitochondria-dependent pathway (type-II). As in the previous case study, the classical PSA and iPSA were applied to this network to elucidate the pathway dominance.
The classical and iPSA analyses were calculated under a constant FasL stimulation (FasL = 2 nM) over the time range of 10,000 seconds. The rankings of the important parameters that control caspase-3 according to the PSA are shown in Figure 6, using consolidated metrics: infinite norm [44](Figure 6(a)), FIM [43] (Figure 6(b)), and time integral [34](Figure 6(c)), and using the sensitivity magnitudes at switching time (Figure 6(d)) (see Additional file 1: Supplementary Figure S3 for detailed sensitivity rankings). From these rankings, one could not obtain any definitive conclusion regarding the dominance of one pathway over the other. On the contrary, the iPSA sensitivity of caspase-3 in Figure 7 clearly supported a type-II dependent caspase-3 switching with an early type-I dependent activation, in agreement with two previous analyses of this model using the Green's function matrix [47] and model reduction [48]. In Figure 7, the high caspase-3 sensitivities with respect to impulse perturbations to upstream reactions were expected during the initial times, as these served as the early cell death signal response. The cleaving of procaspase-3 was carried out by caspase-8 directly (r5 of type-I) before t = 4000 seconds, after which the mitochondrial pathway (r14 of type-II) became the main route of activating caspase-3 (see Additional file 1: Supplementary Figures S4 and S5 for more detail).
The discrepancy between the PSA and iPSA results can again be explained in the context of persistent versus impulse perturbations. As seen in the PSA parameter rankings in Figure 6 and following the insights offered by the iPSA in Figure 7, the effect of perturbing early response processes, including type-I pathway, was integrated over time in the PSA. For example, the highest ranked parameters in the PSA were associated with the first three reactions, r1 to r3. Such integration masked the dynamical importance of different parameters in this analysis. The conclusion from the iPSA is in agreement with the simulations of KOs of type-I and type-II pathways in Figure 8. While type-II knock-out could describe the caspase-3 activation at early times, it failed to capture the switching behaviour. On the other hand, the type-I knock-out was able to reproduce the switching of caspase-3, albeit with a short delay.
The two examples above illustrate the problem of using the classical PSA in identifying the controlling mechanisms of a dynamical system. Of course, this does not mean that the PSA of dynamical models is incorrect, but rather the interpretation of the sensitivity coefficients should be carefully managed. In particular, a large sensitivity magnitude with respect to a parameter suggests the importance of this parameter in the time period between the perturbation time τ and the state observation time t. In contrast, the iPSA is developed with dynamics in mind, where the impact of a single perturbation on the system is realized only at the perturbation time and subsequently they are delivered at varying perturbation times. By doing so, the iPSA coefficients can elucidate the way system dynamics x(t) is achieved, by indicating which and when parameters or processes are essential. Because of the persistent nature of perturbations used in the PSA, it is still not possible to reproduce the conclusions of the iPSA by varying the time of perturbation (see Additional file 1: Supplementary Figures S6 and S7).