pathPSA: A Dynamical Pathway-Based Parametric Sensitivity Analysis

Article pubs.acs.org/IECR

pathPSA: A Dynamical Pathway-Based Parametric Sensitivity Analysis Thanneer Malai Perumal† and Rudiyanto Gunawan*,‡ †

Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Esch/Alzette 4362, Luxembourg Institute for Chemical and Bioengineering, ETH Zurich, Zurich 8093, Switzerland

‡

S Supporting Information *

ABSTRACT: Normal functioning of biological systems relies on coordinated activities of biomolecules participating in a complex network of biological interactions. As human intuition is often incapable in analyzing how biological functions emerge from such complex interactions, the use of mathematical models and systems analysis tools has become critical in understanding biological system behavior. While biological functions have been associated with the connectivity (structure) and pathways of the network and the kinetics of the processes involved, most model analyses address each of these aspects separately. In contrast, we present a novel sensitivity analysis, called pathway-based parametric sensitivity analysis (pathPSA), which combines the analysis of both network structure and kinetics. Unlike existing sensitivity analyses, pathPSA relies on perturbing the kinetics of pathways in the network, using persistent perturbations or impulse perturbations at varying time. Consequently, the sensitivity coefficients can give insights on the dominant pathways in a network and any transient shift in the rate controlling mechanism. The efficacy of pathPSA is demonstrated through an application to understand competing signaling pathways in programmed cell death network.

1. INTRODUCTION A complex system is described as a system that is (1) complicated or intricate, (2) nondeterministic, (3) nonlinear, ill-posed, or chaotic, and (4) predisposed to unexpected outcome, which is referred to as emergent behavior.1 Complex systems have been used to describe problems in different fields of studies, from science and engineering to sociology.2 Examples of complex systems can readily be found among natural and engineered systems, such as social and biological networks, combustion reactions, ecosystems, worldwide-web, and pandemic infections.2 In biology, complexity has been argued to arise as a product of evolution,3 conferring traits such as cellular robustness to common perturbations.4 This robustness however may come at a cost of fragility to rare changes and could turn into an Achilles heel. For instance, the mechanisms that confer the robustness property in our cells can be hijacked, leading to hard-to-treat diseases such as cancer.5 The understanding of complexity, robustness, and its trade-offs in cellular networks has been proposed as the basis of systems oriented drug discovery.6 Advances in high throughput (omics) technologies have enabled scientists to obtain nearly complete information of cellular components, allowing a closer investigation of complexity in cells. At the same time, mathematical modeling has become an indispensable tool for studying emergent behavior in biological networks.7 Numerous kinds of mathematical models have been employed to describe cellular networks, among which ordinary differential equations (ODEs) have become the most popular modeling paradigm.8 The mathematical models are used in lieu of the real system to study, among other things, biological mechanisms that give rise to the emergent functional behavior in networks. Such studies © 2014 American Chemical Society

typically involve performing systems analysis of the models, for example, using sensitivity analysis. Parametric sensitivity analysis (PSA) maps out the parametric dependence of model behavior using sensitivity coefficients.9 The coefficients reflect how network outputs change when one or more model parameters are perturbed. Parameters with high sensitivity magnitudes thus correspond to rate-controlling processes. The PSA has been regularly used to investigate emergent behavior in biochemical networks, such as those in programmed cell death,10 budding yeast cell cycle control,11 IL-6 signaling pathway,12 circadian rhythm models,13 and MAPK/PI3K signal transduction pathway.14 Aside from the above applications, sensitivity analysis has also been used for model calibration and parameter identifiability analysis,15,16 model validation and reduction,15,17 identification of bottlenecking processes,10b,18 and the analysis of cellular robustness.4,19 We have shown elsewhere that classical PSA could lead to an incorrect analysis of dynamical behavior.10b This caveat arises because sensitivity coefficients in the classical PSA are defined using a (Heaviside) step parameter perturbation, and thus changes in model outputs are integrated over the duration of the parameter perturbation. Consequently, any transient shift in the rate-controlling parameters may not be immediately apparent from the classical sensitivity coefficients. This drawback has motivated the development of two sensitivity Special Issue: Massimo Morbidelli Festschrift Received: Revised: Accepted: Published: 9149

October 2, 2013 January 15, 2014 January 16, 2014 January 16, 2014 dx.doi.org/10.1021/ie403277d | Ind. Eng. Chem. Res. 2014, 53, 9149−9157

Industrial & Engineering Chemistry Research

Article

interest in biology, such as metabolic, signalling, and gene regulatory networks 29. A few examples of such ODE models include the modelling of apoptosis30, cell cycle,31 and circadian rhythms.32

analyses for ODE models, namely, the Green’s function matrix (GFM) analysis10c and the impulse parametric sensitivity analysis (iPSA).10b The crucial difference between the GFM/ iPSA and classical PSA is that, in the former, perturbations are introduced at varying time and implemented as impulses, either on the model state (in the GFM analysis) or on the model parameters (in the iPSA). The corresponding (impulse) sensitivity coefficients can directly reveal transient changes in the molecules or parameters that control system behavior. In addition, the GFM analysis and iPSA can be used to study signal propagation through the network. Most of the existing sensitivity analyses, including those mentioned earlier, focus on the effects of perturbing individual states or parameters. However, cellular functions often rely not on a single molecule or process but rather on a group of molecules and reactions belonging to a pathway.20 By using higher order sensitivities, the effects of perturbing multiple processes simultaneously can be analyzed.21 Nevertheless, such an analysis may be limited to only a handful of reactions. Meanwhile, biological network topology or structure, i.e., the connectivity among biological entities, has been shown to be closely related to network function.22 Numerous methods have been developed and used to understand the relationships between network structure and functionality, for example, CellNetAnalyzer,23 simple centrality measures,24 communicability measures,25 SigFlux method,26 simple path analysis,27 and elementary signaling mode (ESM) analysis.27 Many of these methods are based on graph theory, in which biological networks are represented as (directed) graphs. Such analyses do not usually consider the kinetics of cellular pathways. In this work, we present a new method, called pathway parametric sensitivity analysis (pathPSA), for an integrated analysis of network structure and kinetics. The pathPSA consists of (1) identification of independent pathways in a biological network and (2) sensitivity analysis of pathway kinetics. The merits of pathPSA are demonstrated through applications to FasL-induced programmed cell death model.28

3. CLASSICAL PARAMETRIC SENSITIVITY ANALYSIS (PSA) Sensitivity analysis is a well-established systems analysis tool in science and engineering.9,33 This analysis is used to investigate the dependence of model outputs on model parameters through the calculation of sensitivity coefficients.9 In the classical PSA, the sensitivity coefficients are defined as Si , j(t , τ ) =

= lim

Δpj → 0

x(t0 , p) = x 0

Δpj (τ )

∂g ∂g d S(t ) = S(t ) + ; dt ∂x ∂p

=

∂xi(t ) ∂pj (τ )

(3)

S(t0) =

dx 0 dp

(4)

while the FDM uses a finite difference approximation of the derivatives in eq 3. Many off-the-shelf software packages can provide an integrated and user-friendly computational platform for simulation and analyses of biological ODE models.36 These and other software for sensitivity analysis have been summarized in review articles.37

(1)

In biological models, the state vector x ∈  typically describes the concentrations or levels of biomolecules, such as mRNAs and proteins. Meanwhile, the vector-valued function g(x,p) describes the rate of appearance and disappearance of biomolecules due to cellular processes (e.g. transcription, translation, phosphorylation, and dephosphorylation). In the pathPSA, the function g(x,p) is further expanded to n

g (x , p) = Nr(x , p)

Δxi(t , p + Δpj ej)

where xi is the ith system state, pj is the jth model parameter, and ej is the jth column of the m × m identity matrix. The magnitude of sensitivity coefficients is thus related with how much individual parametric perturbation in a system at time τ affects the model output behavior at time t(t ≥ τ). Consequently, one of the common uses of PSA in systems biology is to infer the importance (or lack of importance) of cellular processes and to provide mechanistic explanations for biological behavior (e.g., oscillations).12,34 In eq 3, the sensitivity coefficients depend on the model parameter values, and hence the classical PSA above is considered to be a local analysis. In general, the parametric perturbation can be introduced at any time τ(t0 ≤ τ ≤ t),35 but τ is commonly set to be the initial time t0. Therefore, the argument τ can be dropped out of eq 3 and the sensitivity coefficients are computed for varying time t. A variety of approaches have been developed to calculate sensitivity coefficients for ODE models, including finite difference method (FDM), direct differential method (DDM), and the Green’s function method.9 The DDM and Green’s function method rely on taking a derivative of the model with respect to p, giving

2. MATHEMATICAL MODELS One of the most commonly used paradigms for modelling dynamical biological systems depends on ordinary differential equations. An ODE model describes the transient changes in the system states according to dx(t , p) = g (x , p); dt

change in xi at time t change in pj at time τ

4. IMPULSE PARAMETRIC SENSITIVITY ANALYSIS (IPSA) In many studies, the classical PSA is used to rank parameters based on the magnitude of sensitivity coefficients, either taken at a specific time or using consolidated sensitivity metrics, such as time-integral, average, or norm of sensitivity coefficients.10b,38 The parameter ranking is subsequently used to make conclusions about the important mechanisms or the property of the biological system behavior (such as robustness). We have shown previously that the classical PSA cannot reveal transient changes in the controlling parameters.10b Importantly, when these sensitivity coefficients are used to infer the

(2)

where N denotes the stoichiometric or connectivity matrix and the vector-valued function r(x,p) comprises the rate equations. Finally, the rates of cellular processes depend on a set of kinetic parameters p ∈ m. In such model formulation, the network structure can be analyzed based on the matrix N, while the analysis of model dynamics will be based on r(x,p). The ODE formulation above is general enough to model most systems of 9150

dx.doi.org/10.1021/ie403277d | Ind. Eng. Chem. Res. 2014, 53, 9149−9157


Article

Figure 1. A simple network model. (a) A simple network with 7 species and 10 reactions. Straight arrows connect substrates to products. (b) Reformulated network model after introducing specific pathway perturbation parameters.

described in detail, using a signal transduction network in Figure 1 for illustration purposes. Pathway Identification. The topology of biological networks is often drawn as a directed graph (digraph) G(V,E), in which the vertices (or nodes) V represent biological species and the directed edges E represent biological interactions or reactions. When the reactions (edges) involve multiple substrates and products, the networks can be described by a generalized graph, called hypergraph. The pathway identification can be achieved either manually based on biological knowledge or systematically using pathway finding methods, such as elementary mode analysis,39 extreme pathway analysis,40 simple path analysis, and elementary signaling mode analysis.27 An automated pathway finding would be necessary when the size of the network becomes prohibitively large for a manual curation of pathways. Here, we have used the concept of elementary signaling modes (ESMs),27 but depending on the particular application, the pathPSA could be easily modified to use other definitions of pathways. An ESM is defined as the minimal set of molecules (nodes) that can perform signal transduction from the input (stimulus) to the output (response). The minimality signifies that ESMs comprise a nonredundant set of molecules. Moreover, ESMs account for both directionality and synergistic interactions and thus are also applicable to hypergraphs.27 The identification of ESMs involves rewriting a given network graph in an enriched representation and applying a depth-first search algorithm of the enriched representation for given input and output nodes. For the simple network shown in Figure 1a, the ESMs from x1 (input) to x4 (output) consist of {r1,r2,r3} and {r1,r4,r5,r3}. As shown in this example, a reaction can appear in more than one ESM. Because of combinatorial aspects, the computational cost of finding ESMs increases rapidly with the size of the network, and the number of ESMs in a large network is difficult to count. By definition, ESMs do not include any feedback loop in the network. Moreover, any side processes (e.g., degradation reactions pointing to null in Figure 1a) are also not accounted. While these motifs may not be relevant in the context of signal transduction topology, the kinetics of the associated processes could be important in determining the dynamics of the system response. Thus, the pathway identification step in pathPSA also includes the mapping of possible feedback and side processes that are associated with each ESM. In this case, the algorithm for ESM identification is applied in the reverse direction. For every ESM, the depth first search algorithm is used to identify any (feedback) path connecting a downstream node to an upstream node in an ESM. In addition, finding side processes becomes equivalent to finding synergistic nodes that are one node distance away from an ESM. For example, for the network in Figure 1a, the set {r6,r7,r8} and {r9} are the common

controlling mechanisms of dynamical system behavior (such as oscillations and switching), the parameter rankings from the classical PSA may become misleading. The reason stems from the fact that, in the classical PSA, the sensitivity coefficients are defined using a persistent perturbation on the model parameters. Thus, the sensitivity magnitudes correspond to the integrated impact of parameter perturbations up to time t (see eq 3). Therefore, while the classical PSA can reveal parameter importance over duration of time, the sensitivities do not directly indicate when these perturbations matter. In order to overcome this drawback, we presented the iPSA in a previous publication, which involves perturbing parameters using impulses and at varying times.10b By using impulses, a parameter perturbation affects the system output only at a specific time. Thus, the impulse sensitivity coefficients provide information not only on which model parameters matter but also on when they matter. In contrast to the classical PSA, the impulse sensitivity coefficient iSi,j(t,τ) is defined as the change in the state xi at time t due to an impulse perturbation on the parameter pj at time τ. The iPSA coefficients are calculated as follows:10b n

iSi , j(t , τ ) =

∑ k=1

∂xi(t ) ∂gk (τ ) ∂xk(τ ) ∂pj

(5)

for varying observation and perturbation times (t and τ, respectively) with t ≥ τ. In eq 5, the state sensitivity coefficients Sxi,k(t,τ) = ∂xi(t)/∂xk(τ) are calculated by solving:10c ∂g x d x S (t , τ ) = S (t , τ ); dt ∂x

S x (τ , τ ) = I n × n

(6)

where In×n is an n × n identity matrix.

5. PATHWAY-BASED PARAMETRIC SENSITIVITY ANALYSIS (PATHPSA) The fundamental difference between pathPSA and existing sensitivity analyses is the consideration of pathways. The pathPSA is based on parametric perturbations, but the perturbations are rendered on pathways in the system rather than individual reactions. In this case, a pathway is defined as a set of species and processes (e.g., reactions) that can carry out system response to a stimulus. Analogous to PSA, the sensitivity coefficients in the pathPSA are defined as the change in the state xi at time t due to a perturbation on the jth pathway at time τ, respectively. The pathPSA coefficients thus portray the importance of pathways in regulating the observed output behavior. The computation of the pathPSA sensitivity coefficients consists of three steps: (i) pathway identification, (ii) model reparameterization, and (iii) sensitivity analysis on the reparameterized model. Below, each of these steps is 9151



Article

Here, the pathway parameters αj have nominal values of 1. Analogous to the interpretation of common PSA, a large (small) sensitivity magnitude indicates important (unimportant) pathway. Since iSpath i,j (t,τ) involves the time of perturbation τ and the time of observation t, each coefficient can be visualized as a heat-map, as shown in Figure 2. The heat-map illustrates the

feedback loop and side process, respectively, between the two ESMs. Meanwhile, {r10} is a side process of the ESM {r1,r4,r5,r3}. Finally, the pathways in the pathPSA are formed by generating all possible combinations of each ESM with its feedback loop(s) and side process(es). These pathways can be further curated depending on the purpose of the study. Model Reparameterization. The pathway perturbations are performed by scaling the rate equations of all processes in a pathway simultaneously and in equal proportion. More specifically, a surrogate parameter αk is introduced as a premultiplicative factor to every rate equation in r(x,p) that is involved in the kth pathway. The reparameterization of the ODE model in eq 1 gives d x(t , α ) = Nr *(x , α) = g*(x , α) dt

(7)

where α ∈ npath and npath is the number of pathways. The nominal value of each αk is set to 1, and thus the solution of the reparameterized ODE model stays the same as the original model. As a reaction can participate in more than one pathway (e.g., r3 appears in all pathways in Figure 1a), multiple α’s can appear in a rate equation of r*(x,α). Consider the simple network in Figure 1a. As an illustration, we define two pathways of interest in the network, the first involving {r 1 ,r 2 ,r 3 ,r 6 ,r 7 ,r 8 } and the second involving {r1,r4,r5,r3,r6,r7,r8} (i.e., each ESM in the network combined with the feedback loop). Correspondingly, the network model is reformulated by premultiplying {r1,r2,r3,r6,r7,r8} with α1 and {r1,r4,r5,r3,r6,r7,r8} with α2. Figure 1b shows the reparameterized network in this example. Pathway Perturbation. The third and final step in the pathPSA is to perform parametric sensitivity analysis of the reparameterized ODE model with respect to the pathway parameters α. The sensitivity analysis can be based on persistent or time-varying impulse perturbations, referred to as persistent and impulse pathPSA, respectively. Following eqs 4 and 5, the persistent pathPSA sensitivity coefficients Spath i,j (t) are computed by solving the differential equation

Figure 2. Heat-map of impulse pathPSA sensitivity coefficient. The heat-map depicts the (i,j)th pathPSA coefficient, showing the change in the state xi with respect to a perturbation on the jth pathway parameter αj. The x-axis shows the perturbation time τ at which an impulse perturbation is introduced, while the y-axis represents the observation time t at which the changes of xi are observed.

perturbation−output relationship, where a positive (negative) iSpath i,j (t,τ) implies that a positive impulse perturbation to the jth pathway at time τ will cause an increase (decrease) in xi at time t. Therefore, the dynamical aspect of the system output regulation by pathways can be deduced from the peak/trough of the corresponding sensitivities. The lower right half of the plot is null as the system is assumed to be causal, i.e., the perturbation will only affect the states at time t ≥ τ.

∂g* path ∂g* d d x(t , α ) d S (t ) + ; = Spath (t ) = dα dt dt ∂x ∂α Spath (t0) = 0n × n

6. CASE STUDY: APPLICATION TO FASL-INDUCED APOPTOTIC CELL DEATH The efficacy of the pathPSA is demonstrated by an application to a model of cell death in human Jurkat T-cell lines.30 As shown in Figure 3, the cell death network consists of 28 molecular species (nodes) and 32 reactions (edges), where the model parameters have been obtained through a parameter fitting of experimental data. The output of interest is caspase-3, which is a protease that cleaves many protein substrates.41 The activation of caspase-3 resembles a bitoggle switch behavior in response to Fas ligand, as shown in Figure 3 (see inset), with a switching time of t = 6000 s. The network includes two major routes of caspase-3 activations: mitochondria-independent (type-I) and mitochondria-dependent route (type-II). In this case study, the pathPSA sensitivities are used to elucidate the key pathways in the activation of caspase-3 by FasL. In the first step of the pathPSA (pathway identification), two ESMs (type-I and type-II) and one feedback (involving caspase-6) and four side processes (i.e., DISC inhibition by FLIP, caspase-3 inhibition by XIAP, mitochondria inhibition by Bcl-2, and the apoptosome inhibition by XIAP) were identified, as shown in Figure 3. In addition to these processes, one

(8)

and the coefficients of impulse pathPSA iSi,jpath(t,τ) are evaluated using n

iSipath , j (t , τ ) =

∑ k=1

∂xi(t ) ∂gk* (τ ) ∂xk(τ ) ∂αj

(9)

where 0n×n denotes the n × n zero matrix. Note that the calculation of sensitivities in eq 9 requires integrating the ODE in eq 6. The sensitivities in eq 8 could also be obtained from the solution of eq 6 through the (adjoint) Green’s function method.9 For the ranking of pathways, the sensitivity coefficients should be normalized as follows: path Si̅ path , j (t ) = Si , j (t )

αj xi(t )

path iSi̅ path , j (t , τ ) = iSi , j (t , τ )

(10)

αj(τ ) xi(t )

(11) 9152



Article

Figure 3. Model of FasL-induced programmed cell death in Jurkat T cell lines. The network consists of two ESMs: mitochondria-independent type-I (T1) and mitochondria dependent type-II (T2), a feedback involving caspase-6 (FB), and five side processes: Flip inhibition of DISC (DI), caspase3 inhibition by XIAP (C3I), mitochondrial inhibition by Bcl-2 (MI), apoptosome inhibition by XIAP (AI), and XIAP inhibition by Smac (XI). The rate equations and parameters are available in the original publication by Hua et al. Caspase-3 activation follows a switch-like behavior in response to a constant FasL stimulus of 2 nM (see inset).

pathways differed only during early times. The consolidated sensitivity metrics with respect to type-I clusters using 2-norm were higher than those with respect to type-II clusters, indicating a type-I dominant activation of caspase-3 in agreement with the conclusion from a previous analysis using path the classical PSA.10c However, the time profiles of SCaspase ̅ ‐ 3, j(t ) in Figure 5 indicated that the activation of caspase-3 depended on the type-I pathway in early times (t < 4000 s) and on the type-II pathway later on, including the switching period. The shift from type-I to type-II pathway above agreed with previous analyses using the iPSA and GFM analysis.10b,c Thus, the persistent pathPSA could not provide a definitive conclusion regarding the relative importance of type-I and -II pathways. The impulse pathPSA revealed a different insight than past analyses of the model and the persistent pathPSA above. For path iSCaspase ̅ ‐ 3, j(t , τ ), we obtained only three clusters (see Figure 4b): two clusters involving type-I pathway (clusters B1 and B2) as before, but only one cluster involving type-II pathway (cluster B3). The clusters involving type-I pathway were path identical to those for SCaspase ̅ ‐ 3, j(t ) in Figure 4a. The heat-maps of the average sensitivities of each cluster are shown in Figure 6. The high sensitivities associated with cluster B1 near t = τ in Figure 6a indicated that type-I pathway was triggered early in the response to Fas ligand, leading to an immediate effect (production) in active caspase-3. The sensitivities with respect to pathway perturbations of cluster B2 also peaked at early perturbation time, but the effects on caspase-3 were more lasting than that of cluster B1 (as indicated by high sensitivity magnitudes in the region t > τ in Figure 6b). Taken together, the sensitivities above suggested that XIAP could effectively inhibit the initial caspase-3 activation, by sequestering caspase-3 produced through type-I pathway (see Figure 3).

reaction (Smac inhibition of XIAP) was added as a side process for completeness. Combining each of the two ESMs with the appropriate feedback and side processes led to a total of 72 different pathways, with 8 involving type-I and 64 involving type-II. The detailed list of pathways is provided in Supporting Information Section A.1. After reparameterisation of the ODE model, the pathPSA was carried out under a constant FasL stimulation (FasL = 2nM) over a time range of 10 000 s, during which caspase-3 activation was completed. For the purpose of this case study, path the pathPSA sensitivity coefficients of caspase-3 (i.e., SCaspase ̅ ‐ 3, j path and iSCaspase ̅ ‐ 3, j ) with respect to perturbations on each of the 72 pathways were computed (see Supporting Information Figures A.1−5). As many of these pathways had overlapping processes, some sensitivity coefficients showed strong correlation with each other. Therefore, a clustering of the pathPSA sensitivities was performed based on their cosine distances. Figure 4 shows the resulting binary classification tree of the pathPSA path path path coefficients SCaspase ̅ ̅ ̅ ‐ 3, j (t) and iSCaspase ‐ 3, j (t), ‐ 3, j (t,τ). For SCaspase we identified four clusters: two clusters for type-I (clusters A1 and A2) and two clusters for type-II pathway (clusters A3 and A4), as shown in Figure 4a. The main difference between the two clusters of type-I as well as of type-II pathway was the inclusion of caspase-3 inhibition by XIAP (in clusters A2 and A4). The addition of the feedback involving caspase-6 and other side processes to the two ESMs did not appreciably change the pathPSA sensitivity coefficients, in agreement with previous analyses of the model using the iPSA and GFM analysis. As shown in Figure 5 (see Supporting Information Figure A.2 for the complete sensitivity profiles), the average path SCaspase ̅ ‐ 3, j(t ) of the two clusters within type-I and type-II

9153



Article

outcome could also be obtained by using infinite-norm and integrated sensitivities as the consolidation metric (see Supporting Information Figure A.6). Thus, the impulse pathPSA suggested that type-II pathway persistently dominated over type-I in activating caspase-3, which was in contrary to the outcomes of other sensitivity analyses of the model. The difference between the conclusions from the impulse pathPSA and the other sensitivity analyses, including iPSA,10b GFM,10c and persistent pathPSA above, demonstrated the importance of combining pathway and time-varying impulse perturbations in analyzing the dynamics of biological networks. The iPSA and GFM analyses, which were based on perturbing individual reactions and species at different times, could not reveal the relative importance of pathways in a network (at least not directly). Meanwhile, the persistent pathPSA sensitivities could not be used to discriminate direct and indirect effects of the perturbations as the effects of the perturbations were integrated. For example, as shown in Figure 5, the sensitivities path SCaspase ̅ ‐ 3, j(t ) with respect to type-I pathway perturbation increased quickly because of the fast response of type-I pathway to FasL stimulation in the model. However, the path pathPSA sensitivities iSCaspase ̅ ‐ 3, j(t , τ ) showed that impulse perturbations to type-II pathway caused larger changes in caspase-3 activation than type-I pathway, regardless of the time of the perturbations (see Figure 7 and Supporting Information Figure A.6). In this case, the effects of type-II pathway perturbations during the early time period were delayed (see Figure 6c) and therefore were not immediately accounted for in path the sensitivities SCaspase ̅ ‐ 3, j(t ) (see Figure 5). In agreement with

Figure 4. Dendrogram of pathPSA sensitivity coefficients. Binary classification tree based on cosine distances among normalized pathPSA coefficients of active caspase-3 (a) persistent pathPSA and (b) impulse pathPSA. T1: type-I, T2: type-II, FB: Caspase-6 feedback, DI: Flip inhibition of DISC, C3I: Caspase-3 inhibition by XIAP, MI: Mitochondrial inhibition by Bcl-2, AI: Apoptosome inhibition by XIAP and XI: XIAP inhibition by Smac. The clustering was performed using a cosine distance threshold of 10−3.

path the results from iSCaspase ̅ ‐ 3, j(t , τ ), a simpler model of the cell death signaling network incorporating only type-II pathway could reproduce the switch-like caspase-3 response to FasL.10b

7. DISCUSSION Understanding nonintuitive behavior of biological networks has been done using a variety of systems analysis based on a mathematical description of these networks. In existing analyses, the structure or topology of the network is often studied separately from the dynamics of the network. However, both structure and dynamics have been tied to network function and emergent behavior. In this work, we presented a novel sensitivity analysis, called pathPSA, which considers both structural and kinetic information of biological networks. The pathPSA involves two main components: (1) delineation of relevant pathways in a network graph using the concept of elementary signaling mode (ESM) and (2) perturbation of pathways through a reparameterization of the kinetic ODE model of the network. Specifically, the pathPSA directly gives insights about the pathways that are critical in regulating system behavior, and by employing time-varying impulse perturbation, the specific timing when these pathways matter can also be obtained. The information gained from the pathPSA sensitivities can be used in various applications, such as robustness analysis, model identification and model reduction, drug discovery research for understanding the mechanism of drug action, and identifying the target pathways. The choice of using the pathPSA over other sensitivity analysis (i.e., classical PSA, GFM, and iPSA10b,c) depends on the end applications. In particular, the pathPSA is especially advantageous in the analysis of a large scale network for which pathways of interest

Figure 5. Persistent pathPSA sensitivities of active caspase-3. The normalized sensitivity coefficients of caspase-3 were averaged over the pathways in each cluster: A1 (*), A2 ( × ), A3 (□), and A4 (●). The consolidated sensitivity metrics using 2-norm (over time t) for cluster A1, A2, A3, and A4 were 21.44, 20.55, 18.16, and 17.67, respectively.

Nevertheless, caspase-3 activation was most sensitive to (impulse) perturbations of type-II pathway, as shown in Figure 6c. Moreover, the consolidated sensitivities of type-II pathway using 2-norm (over the observation time t) were higher than those of type-I pathway for all perturbation times. The same 9154



Article

Figure 6. Impulse pathPSA sensitivities of active caspase-3. The heat-maps show the average of normalized pathPSA sensitivity coefficients over the pathways within cluster: (a) B1, (b) B2, and (c) B3. The x-axis shows the time τ at which impulse perturbations are applied, while the y-axis indicates the observation time t. Each plot is scaled to have values between −1 and +1 by the maximum magnitude of the sensitivities coefficients, shown in the lower right corner of each heat-map.

for the ESMs. Subsequently, ESMs that have “small” sensitivities need not be further considered in the full pathPSA. In the case study, the most useful pathPSA sensitivities were based on time-varying impulse perturbations of pathways. Like the iPSA and GFM analysis, the sensitivities possessed two time-axes, perturbation and observation times, thereby allowing one to discriminate between direct and indirect effects of the perturbation. While such sensitivities could be presented as heat-maps over the two time-axes, for ease of comparison and ranking of pathways, a consolidating metric was computed as the 2-norm of sensitivity magnitudes over the observation time t (see Figure 7). Depending on the purpose of the analysis, other consolidating metrics could be used, such as the maximum or integral of sensitivity magnitudes over the observation times or the magnitudes at a specific observation time. By comparing the sensitivity metrics as a function of perturbation time τ, one can deduce not only the relative importance of pathways in the network (i.e., pathway ranking) but also any transient changes in the pathway importance. Nevertheless, some information is lost during the consolidation; for example, we will not be able to determine if the effects of a perturbation are immediate or delayed. The pathPSA can be easily modified to accommodate other definitions of pathways, such as elementary modes and extreme pathways,39 and even arbitrary sets of connected reactions. The pathway identification algorithm in the first step of the pathPSA can in fact be substituted with any pathway finding algorithm or with a manual curation of pathways which are relevant in a specific application. The reparameterization of the ODE model is not dependent on the manner by which biological pathways are defined, as this step requires only the indexes of rate equations involved in each pathway. For the cell death signaling case study above, the concept of ESMs was considered appropriate in defining the relevant signaling paths. The major computational cost in the pathPSA is due to the calculation of the sensitivity coefficients in eq 6. The sensitivities in eq 8 could also be solved indirectly using the Green’s function method, which requires solving eq 6.9 By taking advantage of the semigroup property of eq 6, these coefficients can be constructed from the solutions of ∂xi(τk + Δτ)/∂xj(τ); i.e., one only needs to solve the state sensitivity coefficients for one time step Δτ from each τk. The state sensitivities over the whole time duration can be calculated using simple matrix multiplications as described elsewhere.10c

Figure 7. Consolidated impulse pathPSA sensitivities of active caspase3. The consolidation sensitivity metric of clusters B1 (*), B2 ( × ), and B3 (○) was computed as the 2-norm of the impulse pathPSA coefficients over the observation time t.

are known a priori. In this case, the pathPSA would involve analyzing a small set of pathways, instead of the full set of model parameters and states. Moreover, biological species and/ or processes involved in a pathway are often controlled by a cisregulatory element in DNA, i.e., components of pathways are coregulated by a common set of transcription factors. Diseases, such as cancer, often involve deregulation of pathways, for example, through mutations in a set of transcription factors.42 In such a case, the pathPSA better suits the need for studying the impact of perturbation on a pathway level. The pathway identification algorithm considered in this work borrows the concept of ESMs, and extends it to include feedback and side processes. However, in the algorithm described above, we have ignored any side processes that are two-node distance away from an ESM. While the algorithm can be easily modified to identify such interactions, this may lead to a large increase in the number of pathways, as pathways are generated by combining ESMs with the associated feedback loops and side processes. The issue of combinatorial increase in the number of pathways can be avoided by performing a prescreening of ESMs, which involve carrying out pathPSA only 9155



Article

Notes

The application of the pathPSA in the case study took less than 10 min on a dual-core CPU workstation (Intel 6300 @ 1.86 GHz) and 3GB RAM, and roughly 90% of the computational time was involved in solving eq 6. The computational time of pathPSA is mainly associated with the pathway identification and the calculation of sensitivity coefficients, whose complexity directly scales with the size of the network (i.e., the number of nodes and edges). In particular, the computational cost for identifying pathways in pathPSA depends on several network features, including network size and diameter, the number of inputs and outputs, and the number of feedback loops and side reactions. As mentioned above, most of the computing time of pathPSA sensitivities is associated with solving eq 6, and the number of the associated ODEs increases quadratically with the dimension of the states (i.e., the number of nodes). Aside from the computational cost, when applying (impulse) pathPSA to a large network, one may need to analyze and compare a large set of sensitivity heat-maps, which would become difficult and tedious. In the case study, we have overcome such an issue by clustering and consolidating the sensitivities before ranking of pathways.

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS T.M.P. was financially supported by Singapore Millennium Foundation for his Ph.D. The research was also supported by funding from ETH Zurich.

■

ABBREVIATIONS PSA Parametric sensitivity analysis GFM Green’s function matrix iPSA Impulse parametric sensitivity analysis pathPSA Pathway parametric sensitivity analysis

■

8. SUMMARY Biological complexity is associated with both network topology and nonlinear interactions among biomolecules. A novel sensitivity analysis, called pathway-based parametric sensitivity analysis (pathPSA), is presented in this work, for the analysis of kinetic ODE model of biological networks. In contrast to existing sensitivity analyses, the pathPSA relies on perturbations of pathways in the network. The analysis involves a three step process: (1) identifying relevant pathways in the network, (2) reparameterizing the ODE model, and (3) performing parametric sensitivity analysis using either persistent or timevarying impulse perturbations. In the case study of Fas-induced programmed cell death model, the pathPSA sensitivities suggested that while mitochondria-independent (type-I) pathway was triggered early in the cell death response, the activation of caspase-3 depended on mitochondria-dependent pathway (type-II). The effects of perturbing type-II pathway during the early time period, however, appeared after a delay. The presence of caspase-6 feedback mechanism was found to have no effect in the caspase-3 activation. While the pathPSA is created for an application to signal transduction networks, the analysis can be easily modified for analyzing other kinds of biological pathways.

■

ASSOCIATED CONTENT

S Supporting Information *

List of reactions in pathways, the complete persistent and impulse pathPSA sensitivity coefficients, and consolidated impulse pathPSA sensitivities using the integral and maximum sensitivity magnitude as the consolidating metric. This material is available free of charge via the Internet at http://pubs.acs.org.

■

REFERENCES

(1) Foote, R. Mathematics and complex systems. Science 2007, 318 (5849), 410−412. (2) Strogatz, S. H. Exploring complex networks. Nature 2001, 410 (6825), 268−276. (3) Adami, C.; Ofria, C.; Collier, T. C. Evolution of biological complexity. Proc. Natl. Acad. Sci. U. S. A. 2000, 97 (9), 4463−4468. (4) Stelling, J.; Sauer, U.; Szallasi, Z.; Doyle, F. J.; Doyle, J. Robustness of cellular functions. Cell 2004, 118 (6), 675−685. (5) Kitano, H. Cancer robustness: tumour tactics. Nature 2003, 426 (6963), 125. (6) Kitano, H. A robustness-based approach to systems-oriented drug design. Nat. Rev. Drug Discovery 2007, 6 (3), 202−210. (7) Szallasi, Z.; Stelling, J.; Periwal, V. System Modeling in Cell Biology From Concepts to Nuts and Bolts; The MIT Press: Cambridge, MA, 2006. (8) Ideker, T.; Lauffenburger, D. Building with a scaffold: emerging strategies for high- to low-level cellular modeling. Trends Biotechnol. 2003, 21 (6), 255−262. (9) Varma, A.; Morbidelli, M.; Wu, H. Parametric Sensitivity in Chemical Systems; Cambridge University Press: Cambridge, U.K.: 1999. (10) (a) Eissing, T.; Allgower, F.; Bullinger, E. Robustness properties of apoptosis models with respect to parameter variations and intrinsic noise. Systems Biology (Stevenage) 2005, 152 (4), 221−228. (b) Perumal, T. M.; Gunawan, R. Understanding dynamics using sensitivity analysis: caveat and solution. BMC Syst. Biol. 2011, 5 (1), 41. (c) Perumal, T. M.; Wu, Y.; Gunawan, R. Dynamical analysis of cellular networks based on the Green’s function matrix. J. Theor. Biol. 2009, 261 (2), 248−259. (11) Lovrics, A.; Zsély, I. G.; Csikász-Nagy, A.; Zádor, J.; Turányi, T.; Novák, B. Analysis of a budding yeast cell cycle model using the shapes of local sensitivity functions. Int. J. Chem. Kinet. 2008, 40 (11), 710− 720. (12) Chu, Y.; Jayaraman, A.; Hahn, J. Parameter sensitivity analysis of IL-6 signalling pathways. IET Syst. Biol. 2007, 1 (6), 342−352. (13) (a) Jin, Y.; Peng, X.; Liang, Y.; Ma, J. Uniform design-based sensitivity analysis of circadian rhythm model in Neurospora. Comput. Chem. Eng. 2008, 32 (8), 1956−1962. (b) Gunawan, R.; Doyle, F. J. Isochron-based phase response analysis of circadian rhythms. Biophys. J. 2006, 91 (6), 2131−2141. (c) Gunawan, R.; Doyle, F. J. Phase sensitivity analysis of circadian rhythm entrainment. J. Biol. Rhythms 2007, 22 (2), 180−194. (14) Hu, D.; Yuan, J. M. Time-dependent sensitivity analysis of biological networks: coupled MAPK and PI3K signal transduction pathways. J. Phys. Chem. A 2006, 110 (16), 5361−5370. (15) Gadkar, K. G.; Gunawan, R.; Doyle, F. J., III Iterative approach to model identification of biological networks. BMC Bioinformatics 2005, 6, 155. (16) Srinath, S.; Gunawan, R. Parameter identifiability of power-law biochemical system models. J. Biotechnol. 2010, 149 (3), 132−140. (17) (a) Petzold, L.; Zhu, W. Model reduction for chemical kinetics: An optimization approach. AlChE J. 1999, 45 (4), 869−886.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Author Contributions

T.M.P. designed and performed the experiments. R.G. designed and oversaw the experiments. Both authors wrote and proofread this manuscript. 9156



Article

Members. Ann. N.Y. Acad. Sci. 2009, 1158 (1), 1−13. (d) Iwamoto, K.; Tashima, Y.; Hamada, H.; Eguchi, Y.; Okamoto, M. Mathematical modeling and sensitivity analysis of G1/S phase in the cell cycle including the DNA-damage signal transduction pathway. Biosystems 2008, 94 (1−2), 109−117. (35) Turányi, T. Sensitivity analysis of complex kinetic systems. Tools and applications. J. Math. Chem. 1990, 5 (3), 203−248. (36) (a) Rodriguez-Fernandez, M.; Banga, J. R. SensSB: a software toolbox for the development and sensitivity analysis of systems biology models. Bioinformatics 2010, 26 (13), 1675−1676. (b) Schmidt, H.; Jirstrand, M. Systems Biology Toolbox for MATLAB: a computational platform for research in systems biology. Bioinformatics 2006, 22 (4), 514−515. (c) XPPAUT. http://www.math.pitt.edu/∼bard/xpp/xpp. html (accessed August 2013). (37) (a) Alves, R.; Antunes, F.; Salvador, A. Tools for kinetic modeling of biochemical networks. Nat. Biotechnol. 2006, 24 (6), 667− 672. (b) Klipp, E.; Liebermeister, W.; Helbig, A.; Kowald, A.; Schaber, J. Systems biology standards−the community speaks. Nat. Biotechnol. 2007, 25 (4), 390−391. (38) (a) Yue, H.; Brown, M.; He, F.; Jia, J.; Kell, D. B. Sensitivity analysis and robust experimental design of a signal transduction pathway system. Int. J. Chem. Kinet. 2008, 40 (11), 730−741. (b) Yue, H.; Brown, M.; Knowles, J.; Wang, H.; Broomhead, D. S.; Kell, D. B. Insights into the behaviour of systems biology models from dynamic sensitivity and identifiability analysis: a case study of an NF-kappaB signalling pathway. Mol. BioSyst. 2006, 2 (12), 640−649. (39) Stelling, J.; Klamt, S.; Bettenbrock, K.; Schuster, S.; Gilles, E. D. Metabolic network structure determines key aspects of functionality and regulation. Nature 2002, 420 (6912), 190−193. (40) Price, N. D.; Reed, J. L.; Papin, J. A.; Wiback, S. J.; Palsson, B. O. Network-based analysis of metabolic regulation in the human red blood cell. J. Theor. Biol. 2003, 225 (2), 185−194. (41) Reed, J. C.; Doctor, K. S.; Godzik, A. The domains of apoptosis: a genomics perspective. Science Signaling: STKE 2004, 2004 (239), re9. (42) Goodarzi, H.; Elemento, O.; Tavazoie, S. Revealing global regulatory perturbations across human cancers. Mol. Cell 2009, 36 (5), 900−911.

(b) Perumal, T. M.; Madgula Krishna, S.; Tallam, S. S.; Gunawan, R. Reduction of kinetic models using dynamic sensitivities. Comput. Chem. Eng. 2013, 56 (0), 37−45. (18) (a) Ingalls, B. Sensitivity analysis: from model parameters to system behaviour. Essays Biochem. 2008, 45, 177−193. (b) Rabitz, H.; Kramer, M.; Dacol, D. Sensitivity analysis in chemical kinetics. Annu. Rev. Phys. Chem. 1983, 34, 419−461. (19) Perumal, T. M.; Yan, W.; Gunawan, R. Robustness Analysis of Cellular Systems for In Silico Drug Discovery. Proceedings of the 17th World Congress - The International Federation of Automatic Control (IFAC), Seoul, Korea, 2008; IFAC: 2008; pp 12607−12612. (20) Barabasi, A.-L.; Oltvai, Z. N. Network biology: understanding the cell’s functional organization. Nat. Rev. Genet. 2004, 5 (2), 101− 113. (21) Saltelli, A.; Tarantola, S.; Campolongo, F.; Ratto, M. Sensitivity Analysis in Practice: A Guide to Assessing Scientific Models; John Wiley & Sons, Ltd.: 2004. (22) Zhu, X.; Gerstein, M.; Snyder, M. Getting connected: analysis and principles of biological networks. Genes Dev. 2007, 21 (9), 1010− 1024. (23) Klamt, S.; Saez-Rodriguez, J.; Gilles, E. D. Structural and functional analysis of cellular networks with CellNetAnalyzer. BMC Syst. Biol. 2007, 1, 2. (24) (a) Wasserman, S.; Faust, K. Social network analysis: methods and applications; Cambridge University Press: Cambridge; New York, 1994. (b) Perra, N.; Fortunato, S. Spectral centrality measures in complex networks. Phys. Rev. E 2008, 78 (3), 036107. (25) Crofts, J. J.; Higham, D. J. A weighted communicability measure applied to complex brain networks. J. R. Soc., Interface 2009, 6 (33), 411−414. (26) Liu, W.; Li, D.; Zhang, J.; Zhu, Y.; He, F. SigFlux: a novel network feature to evaluate the importance of proteins in signal transduction networks. BMC Bioinformatics 2006, 7, 515. (27) Wang, R. S.; Albert, R. Elementary signaling modes predict the essentiality of signal transduction network components. BMC Syst. Biol. 2011, 5, 44. (28) Hua, F.; Hautaniemi, S.; Yokoo, R.; Lauffenburger, D. A. Integrated mechanistic and data-driven modelling for multivariate analysis of signalling pathways. J. R. Soc., Interface 2006, 3 (9), 515− 526. (29) (a) Fell, D. A. Metabolic control analysis: a survey of its theoretical and experimental development. Biochem. J. 1992, 286 (Part2), 313−330. (b) Conrad, E. D.; Tyson, J. J. Modeling molecular interaction networks with nonlinear ordinary differential equations. In Systems Modeling in Cellular Biology: From Concepts to Nuts and Bolts; Szallasi, Z., Stelling, J., Periwal, V., Eds.; MIT Press: Cambridge, MA, 2006; pp 97−123. (30) Hua, F.; Hautaniemi, S.; Yokoo, R.; Lauffenburger, D. A. Integrated mechanistic and data-driven modelling for multivariate analysis of signalling pathways. J. R. Soc., Interface 2006, 3 (9), 515− 526. (31) Novák, B.; Tyson, J. J. A model for restriction point control of the mammalian cell cycle. J. Theor. Biol. 2004, 230 (4), 563−579. (32) Forger, D. B.; Peskin, C. S. A detailed predictive model of the mammalian circadian clock. Proc. Natl. Acad. Sci. U. S. A. 2003, 100 (25), 14806−14811. (33) (a) Saltelli, A.; Chan, K.; Scott, E. M. Sensitivity Analysis: Gauging the Worth of Scientific Models; John Wiley & Sons, Ltd.: 2000; (b) Saltelli, A.; Ratto, M.; Tarantola, S.; Campolongo, F. Sensitivity analysis for chemical models. Chem. Rev. 2005, 105 (7), 2811−2827. (34) (a) Ihekwaba, A. E.; Broomhead, D. S.; Grimley, R. L.; Benson, N.; Kell, D. B. Sensitivity analysis of parameters controlling oscillatory signalling in the NF-kappaB pathway: the roles of IKK and IkappaBalpha. Systems Biology (Stevenage) 2004, 1 (1), 93−103. (b) Ihekwaba, A. E.; Broomhead, D. S.; Grimley, R.; Benson, N.; White, M. R.; Kell, D. B. Synergistic control of oscillations in the NFkappaB signalling pathway. Systems Biology (Stevenage) 2005, 152 (3), 153−160. (c) Adler, P.; Peterson, H.; Agius, P.; Reimand, J.; Vilo, J. Ranking Genes by Their Co-expression to Subsets of Pathway 9157


pathPSA: A Dynamical Pathway-Based Parametric Sensitivity Analysis

Recommend Documents