Intervention Pathway Discovery via Context-Dependent Dynamic Sensitivity Analysis

TL;DR

This paper introduces a comprehensive framework for analyzing the sensitivity of biological networks, combining static and dynamic methods to identify key intervention pathways and influential elements in context-dependent scenarios.

Contribution

The paper presents a novel framework that integrates static and dynamic sensitivity analysis for discrete biological network models, enabling context-aware identification of intervention pathways.

Findings

Dynamic analysis captures context-dependent transient states.

Framework effectively identifies influential pathways for interventions.

Application to T cell differentiation demonstrates practical utility.

Abstract

The sensitivity analysis of biological system models can significantly contribute to identifying and explaining influences of internal or external changes on model and its elements. We propose here a comprehensive framework to study sensitivity of intra-cellular networks and to identify key intervention pathways, by performing both static and dynamic sensitivity analysis. While the static sensitivity analysis focuses on the impact of network topology and update functions, the dynamic analysis accounts for context-dependent transient state distributions. To study sensitivity, we use discrete models, where each element is represented as a discrete variable and assigned an update rule, which is a function of element's known direct and indirect regulators. Our sensitivity analysis framework allows for assessing the effect of context on individual element sensitivity, as well as on element criticality in reaching preferred outcomes. The framework also enables discovery of most influential pathways in the model that are essential for satisfying important system properties, and thus, could be used for interventions. We discuss the role of nine different network attributes in identifying key elements and intervention pathways, and evaluate their performance using model checking method. Finally, we apply our methods on the model of naive T cell differentiation, and further demonstrate the importance of context-based sensitivity analysis in identifying most influential elements and pathways.

Figures

Equations

$$\frac{\partial f_j}{\partial x_i}=f_j{|}_{x_i=0}\oplus f_j{|}_{x_i=1}=f_{ji}(x_1,x_2,...,x_{i-1},x_{i+1},...,x_N)$$

$${\alpha }_i^j={\alpha }_i^{f_j}=E(\frac{\partial f_j}{\partial x_i})$$

$$s^j(\mathbf{v}=\mathbf{v})=\sum _{i=1}^{k(j)}\frac{\partial f_j}{\partial x_i}({\mathbf{v}}^i={\mathbf{v}}^i)$$

$$s^j=E(s^j(\mathbf{v}=\mathbf{v}))=\sum _{i=1}^{k(j)}E(\frac{\partial f_j}{\partial x_i}({\mathbf{v}}^i={\mathbf{v}}^i))=\sum _{i=1}^{k(j)}{\alpha }_i^j$$

$${\alpha }_i^j={\alpha }_i^{f_j}=E(\frac{\partial f_j}{\partial x_i})=\frac{1}{2^{k(j)}}\sum _{\mathbf{v}}\frac{\partial f_j(\mathbf{v})}{\partial x_i}$$

$$pe(\mathbf{v})=\frac{n(T_{\mathbf{v}})}{n(T)}$$

$$p(\mathbf{v})=\frac{pe(\mathbf{v})}{{\sum }_{\mathbf{v}}pe(\mathbf{v})}$$

$${\alpha }_i^j={\alpha }_i^{f_j}=E(\frac{\partial f_j}{\partial x_i})=\sum _{\mathbf{v}}\frac{\partial f_j(\mathbf{v})}{\partial x_i}p(\mathbf{v})$$

$$\{x_s,x_{p_1},...,x_{p_m},x_t\}\triangleq \underset{x_s,...,x_t}{\mathrm{argmin}}{w}_{s{p}_1}+\sum _{i=1}^{m-1}{w}_{p_i{p}_{i+1}}+w_{p_{m}t}$$

$$s(path)=exp(-[w_{s{p}_1}+\sum _{i=1}^{m-1}{w}_{p_i{p}_{i+1}}+w_{p_{m}t}])$$

Full text

1 Introduction

Understanding sensitivity is a critical step in studying system’s adaptability to environment, and its robustness against perturbations, both of which are considered indispensable in living organisms (Shmulevich et al.,, 2004). In biological systems, the sensitivity assessment can indicate how susceptible the system (e.g., a cell) and its components (e.g., proteins or genes) are to internal or external changes. In turn, the analysis of sensitivity can help locate critical regulatory pathways in the system, and give insights into the robustness and responsiveness to interventions or treatments. To study the sensitivity in biological systems, researchers rely on carefully designed in vivo or in vitro experiments, or on analysis of computational models that can represent and mimic the behavior of the real systems in silico. The sensitivity analysis approaches, and their applications, vary depending on studied systems and model types. In this work, we focus on the sensitivity of biochemical, intra- and inter-cellular networks, that is, signal transduction, metabolic, and gene regulatory networks (Shmulevich et al.,, 2004).

There are several different approaches to modeling biochemical networks, which depend on the size of the network, and even more, on the information available about network components and their interactions. A common approach uses ordinary differential equations (ODEs) derived from the reaction network (Faeder et al.,, 2009). However, reaction networks grow exponentially with the number of network components (e.g., receptors, ligands, kinases, etc.), and while fast methods are available to numerically solve large sets of ODEs, a common issue with this type of modeling is the lack of knowledge about network details necessary for accurate ODE-based models. It is often the case that we are familiar only with indirect cause-effect relationships for a number of interactions in the network, as there is no information (or existing knowledge) about exact mechanisms and the parameters necessary to create ODEs. Therefore, for such systems, a modeling approach is required which can accurately integrate all known interactions, and efficiently extend models with newly found interactions.

The discrete modeling approach has been suggested in the past to overcome the issue of missing information while still providing important insights into system behavior (Albert et al.,, 2014; Miskov-Zivanov et al.,, 2013; Sun et al.,, 2014; Gan et al.,, 2016). Although based on relatively simple discrete formalism, this approach allows for automated integration of both prior knowledge and data into models. In this work, we focus on sensitivity in discrete models, in general, and logical or Boolean models, in particular.

In order to understand and explain the dynamics of biochemical network models, we use sensitivity analysis not only to study the influence from immediate, directly connected regulators, but also to explore global influences from indirectly connected elements. Previous research on sensitivity analysis in models of biochemical networks explored several directions. First, in (Shmulevich et al.,, 2004; Shmulevich et al.,, 2002; Shmulevich et al.,, 2002; Anon,, 2007), the authors studied ‘‘local’’ average sensitivity of a model element and the influence of its regulators on this element, while (Qian et al.,, 2009; Trinh et al.,, 2016) focused on the ‘‘global’’ long-run sensitivity of how likely a mutation is to change the converging attractor (steady-state) of a biological network. In this work, we start with the former direction to explore sensitivity ‘‘locally’’, and then extend the analysis by assigning a sensitivity value to each interaction in the network, and use these values to study sensitivity ‘‘globally’’. Second, considering in particular discrete models of biological networks, sensitivity analysis has been applied only on two types of models, the probabilistic Boolean networks (PBN), which are inferred from gene expression profiles (Shmulevich et al.,, 2002; Shmulevich et al.,, 2002; Qian et al.,, 2009) and random Boolean networks (RBN), randomized according to biological observations (Shmulevich et al.,, 2004; Trinh et al.,, 2016; Kochi et al.,, 2014). To the best of our knowledge, there has been no related sensitivity study on models assembled from individual element interactions using knowledge sources such as published literature, or information from experts, thus we refer to them as knowledge-based models (KBM). In contrast to PBN and RBN models, KBM models allow for various function types in network nodes (i.e., model elements), not only Boolean functions. Finally, previous research assumed static state distributions in studied models, such as uniform state distribution (Shmulevich et al.,, 2004; Anon,, 2007; Kochi et al.,, 2014) or steady-state distribution (Shmulevich et al.,, 2002; Shmulevich et al.,, 2002; Anon,, 2007), neither of which captures transient context-dependent trajectories, and the dynamics of models between initial and steady state. In this work, we propose a more comprehensive sensitivity analysis framework, in which we define element influence and sensitivity with respect to the transient and scenario-dependent state distributions of the system.

We perform both static and dynamic sensitivity analysis, the former assuming uniform state distribution, and the latter using a dynamic distribution estimated from stochastic simulation (Miskov-Zivanov et al.,, 2013) trajectories under a particular scenario. Static sensitivity analysis can offer insights about the topological structure of the model and canalizing degree (Anon,, 2007) of Boolean functions, but it entails no knowledge about the system beyond individual element update rules. Dynamic sensitivity analysis, on the other hand, incorporates the context-based bias into both individual and joint element state distributions at the intermediate points along trajectories. We propose a method to determine these joint distributions from dynamic simulation traces, and we demonstrate the importance of performing dynamic analysis.

Our sensitivity analysis approach enables the extraction of globally important pathways, which is critical in studying biological systems. Once we identify these pathways, we can design control strategies and interventions to tune system inputs or alter element values during the transient process. Through the construction of a weighted directed graph for the model, we have reduced the problem of extracting important paths to a minimum-cost graph search problem.

We also validate our sensitivity analysis results using the model checking approach (Miskov-Zivanov et al.,, 2016; Wang et al.,, 2016) and demonstrate that element influence (especially in the case of dynamic analysis), outperforming many other network attributes, is a useful measure in pathway extraction. As part of the validation process, we propose a procedure to apply our sensitivity analysis framework to biological networks via the best-first search algorithm.

To demonstrate the advantages of our proposed work, we applied our sensitivity analysis framework on the previously published model of the intra-cellular network that controls naive T cell differentiation (Miskov-Zivanov et al.,, 2013). With element influence and sensitivity measures, we evaluate, under different scenarios, individual elements and interactions, as well as the overall model. Finally, we highlight the significance of sensitivity analysis by comparing and contrasting observations from static network-based analysis and context-dependent dynamic analysis.

2 Background and Related Work

In this section, we provide the background of the modeling and simulation approaches that we use as part of our sensitivity analysis framework.

2.1 Discrete Modeling Approaches

In discrete models, each element is assigned a set of discrete values, and an update rule, which is a function of its regulators. The model elements are, therefore, connected into a complex network via their update rules.

Among discrete modeling approaches, Bayesian Networks and Boolean Networks have been two common structures used in studying biological networks, in particular, gene regulatory networks (Kaderali et al.,, 2008; Mihajlovic et al.,, 2001; Pe’er et al.,, 2001; Davidich et al.,, 2008). A Bayesian Network describes the probabilistic dependence interactions among elements. The model is assembled via structure learning (integrating prior knowledge from literature, experts) and parameter learning (Kaderali et al.,, 2008; Mihajlovic et al.,, 2001). However, it fails to capture cyclic interactions and neglects feedback loops (Friedman et al.,, 2000). Boolean Networks (BN) (Kaderali et al.,, 2008), together with Probabilistic Boolean Networks (PBN) (Shmulevich et al.,, 2002; Shmulevich et al.,, 2002) and Random Boolean Network (RBN)s (Shmulevich et al.,, 2004; Trinh et al.,, 2016), aim to model the biological phenomena with simple Boolean formalism, where each element is assigned a Boolean variable to represent its state, and one or multiple candidate Boolean update function. In the BN-based approaches, state transitions are usually assumed to be synchronous, that is, all elements are updated simultaneously. In contrast to Bayesian Networks, BNs allow modeling of cyclic regulation and intertwined feedback loops. On the other hand, BNs inferred from data are highly dependent on expression profiles and cannot incorporate prior knowledge of interactions between elements. In addition, specific to PBNs, their high computational complexity limits the application to large-scale networks.

The discrete modeling approach that we use in this work allows for feedback loops, integration of both prior knowledge and data, more realistic simulation of state transitions, as well as analysis of large hybrid networks that include protein-protein interactions, gene regulation, and even metabolic pathways and cellular processes.

2.2 Modeling Preliminaries

The construction of a model begins with identifying key system components, and their interactions, usually through literature reading, data analysis or discussions with experts (Sayed et al.,, 2017). The extraction of this information from knowledge sources allows modelers to identify the set of model elements $V=\{x_{1},x_{2},...,x_{N}\}$, where $N$ is the number of elements in the model. For each element $x_{i}$ in $V$, we define influence set $V_{IS}^{i}\subset V$ as the set of other elements that regulate it, as well as the polarity (positive or negative) of these regulations, and we also define the influence set size for element $x_{i}$ as $k(i)=|V_{IS}^{i}|$. The influence sets in a model can be illustrated as an influence graph $G(V,E)$, where nodes are model elements from $V$, and directed edges in $E$ are regulatory interactions between elements.

For a given element $x_{i}$, we can also define the number of discrete values representing different levels of the element’s activity, $n_{i}$, such that $x_{i}\in X_{i}:\{0,1,2,...,n_{i}-1\}$. We will call the value assigned to element $x_{i}$ a state of element $x_{i}$, and by assigning values to all elements in the model, we obtain overall model state. Once we define model elements, their influence sets, and levels of activity, we create an executable model that can be simulated, by assigning update functions to a subset (or all) of model elements. Update functions can be derived according to the information available about element regulations, as described in (Miskov-Zivanov et al.,, 2013). Thus, an influence graph $G(V,E)$, together with additional information about regulatory functions, is used to create a discrete model, $M(V,F)$, with $F=\{f_{1},f_{2},...,f_{N}\}$, where, for each model element $x_{i}$, we define a deterministic discrete element update function, $f_{i}$, by mapping a $k(i)$-dimensional non-negative vector to a non-negative integer in the set $X_{i}$. Boolean models are considered a special case of discrete models where the domain of all elements is $B={0,1}$ and the primitive operators used to create any function in $F$ include AND, OR, and NOT. In Fig. 1(a), we show an example three-element Boolean model, its influence graph, and element update functions.

2.3 Model Simulation

Discrete models have been studied either formally (Miskov-Zivanov et al.,, 2016; Wang et al.,, 2016) or using simulations (Sayed et al.,, 2017). Since we focus in this work on discrete models $M(V,F)$, and also account for the randomness in the occurrence time of biochemical reactions, we use DiSH simulator (Sayed et al.,, 2017). DiSH supports several simulation schemes, which can be divided into two groups: simultaneous (SMLN) and sequential (SQ). In the SMLN update scheme (commonly used for BN-based models), current state values of all elements are used to simultaneously compute their next state values, thus, the SMLN simulation scheme is deterministic. In the SQ update scheme, elements are updated sequentially, one after the other. The order of element updates is random, and the update rate for each element depends on the information available about its regulatory mechanisms. While DiSH simulator supports several SQ-based schemes, in this work we use the RSB-SQ (random step-based sequential) scheme, in which, at a given simulation step, only one randomly selected element is updated according to its update rule, while all other elements are not updated during that step. Fig. 1(b) shows the state transition graph for our example model (Fig. 1(a)) under both the SMLN and RSB-SQ update schemes.

In order to obtain a trajectory for the model, that is, state transitions of all elements between initial state and final state, we need to define a simulation scenario. Each scenario includes information about initial values (states) of all non-input model elements (i.e., nodes in the influence graph that have arrows pointing at them), initial values for all model inputs, and, when available, perturbations that are assumed to happen at a particular model element, at a specified time point. The simulation is then executed following the scenario, from the initial state, until a pre-specified final state, which is indicated with the number of simulation steps. One simulation run provides a trajectory of each model element between an initial and a final state.

3 Sensitivity Analysis

In this section, we discuss the details of our methodology. We outline in Fig. 2 the data flow diagram of the framework. We use as inputs a model $M(V,F)$ and a scenario under which the model will be analyzed. The model alone is sufficient for static sensitivity analysis (described in detail in Section 3.1, 3.2, 3.3), while the scenario definition is required for dynamic sensitivity analysis (see Section 3.4). As mentioned in Section 2.3, to obtain model trajectories for the dynamic sensitivity analysis, we run simulations using RSB-SQ scheme. We then define element influence and element sensitivity and show how to compute them. With the influence and sensitivity computation results, we can extend a discrete model graph $G(V,E)$ to a weighted directed graph $G(V,E,W)$, as described in Section 3.5. An overall algorithm of our sensitivity analysis framework is also given in Section 3.5.

3.1 Element Influence

For a given set of model elements $V=\{x_{1},x_{2},...,x_{N}\}$, we define a vector $\boldsymbol{v}=(x_{1},x_{2},...,x_{N})$ to represents the model’s state. We are interested in computing a sensitivity of element $x_{j}$ to changes in the value of element $x_{i}$, where $i,j\in\{1,2,...,N\}$, and $i\neq j$. To find the sensitivity of function $x_{j}=f_{j}(x_{1},x_{2},...,x_{N})$ to element $x_{i}$, that is, the influence of element $x_{i}$ on element $x_{j}$, we need to calculate the partial derivative of function $f_{j}$ with respect to element $x_{i}$. For logical models with Boolean variables and Boolean functions, the partial derivative is defined as (Shmulevich et al.,, 2004):

[TABLE]

where $f_{j}|_{x_{i}=0}$ and $f_{j}|_{x_{i}=1}$ are the co-factors of $f_{j}$ with respect to $x_{i}$, which evaluate $f_{j}$ assuming $x_{i}=0$ and $x_{i}=1$, respectively.

It is obvious from Equation (1) that, $\frac{\partial f_{j}}{\partial x_{i}}$ does not depend on $x_{i}$, and only depends on other model elements $x_{l}\in V,l\neq i$, as determined from the update function $f_{j}$ of element $x_{j}$. Therefore, to find whether element $x_{j}$ can be influenced by $x_{i}$, we need to identify all possible values of vector $(x_{1},x_{2},...,x_{i-1},x_{i+1},...,x_{N})$, for which partial derivative, $\frac{\partial f_{j}}{\partial x_{i}}$, is true (i.e., equal 1). It can be seen from Equation (1), that the partial derivative will be equal 1 if and only if, for given values of $(x_{1},x_{2},...,x_{i-1},x_{i+1},...,x_{N})$, functions $f_{j}|_{x_{i}=0}$ and $f_{j}|_{x_{i}=1}$ have different values. In such cases, the value of function $f_{j}$ changes when $x_{i}$ changes, that is, element $x_{i}$ can influence element $x_{j}$.

The influence of element $x_{i}$ in function $f_{j}$ is defined as expectation of partial derivative of function $f_{j}$ with respect to element $x_{i}$:

[TABLE]

Since the partial derivative $\frac{\partial f_{j}}{\partial x_{i}}$ itself is a Boolean function with two possible values {0,1}, we can consider the expectation $E(\frac{\partial f_{j}}{\partial x_{i}})$ as a probability $P\{(x_{1},x_{2},...,x_{i-1},x_{i+1},...,x_{N})|\frac{\partial f_{j}}{\partial x_{i}}=1\}$, that is, a ratio of the number of vectors $(x_{1},x_{2},...,x_{i-1},x_{i+1},...,x_{N})$ for which $\frac{\partial f_{j}}{\partial x_{i}}=1$ to the number of all possible vectors $(x_{1},x_{2},...,x_{i-1},x_{i+1},...,x_{N})$. This quantitatively describes the relationship between the regulator $x_{i}$ and regulated element $x_{j}$, as it provides the probability that a change in element $x_{i}$ will change the value of $x_{j}$.

In general, the procedure of calculating the influence of variable $x_{i}$ on variable $x_{j}$ (ELEMENT_INFLUENCE($x_{i},x_{j}$) in the algorithm described in Section 3.5), can be divided into two main steps:

(1) COMPUTE_INFLUENTIAL_VECTORS($f_{j},x_{i}$), which takes as arguments update function $f_{j}$ and one of its regulator variables $x_{i}$, returns the set of vectors for which the change in $x_{i}$ leads to a change in $f_{j}$;

(2) EST_PROB($\boldsymbol{v}$), which computes the probability of occurrence of a vector of element values $\boldsymbol{v}$ from simulation trajectories; as default probability, we assume that values follow uniform distribution.

Example 1: Given a Boolean function $x_{5}=f_{5}=x_{1}x_{2}+x_{1}^{{}^{\prime}}x_{3}+x_{2}x_{3}^{{}^{\prime}}x_{4}$, how can we determine the influence of element $x_{1}$ on $x_{5}$? First, we find all combinations of $x_{2},x_{3},x_{4}$ for which $\frac{\partial f_{5}}{\partial x_{1}}=x_{2}x_{3}^{{}^{\prime}}x_{4}^{{}^{\prime}}+x_{2}^{{}^{\prime}}x_{3}=1$. As can be seen from the truth table in Fig. 3(left), these are the value combinations that have probabilities of occurrence $p_{3},p_{4}$, and $p_{5}$. Next, we add the probabilities of these input vectors and compute the influence $\alpha_{1}^{5}=p_{3}+p_{4}+p_{5}$. Similarly, we compute all other influences, $\alpha_{2}^{5},\alpha_{3}^{5},\alpha_{4}^{5}$ (listed in Fig. 3(right)).

3.2 Element Sensitivity

Besides evaluating influence of individual regulators, we also investigate overall sensitivity of elements with respect to the combined influence of all their regulators (i.e., their influence sets). For a given element $x_{j}$, and its influence set $\{x_{j_{1}},x_{j_{2}},...,x_{j_{k(j)}}\}$, where $k(j)$ is the size of the influence set of $x_{j}$, we can decrease the size of vector $\boldsymbol{v}$ to only include $x_{j}$’s influence set elements, such that $\boldsymbol{v}=(x_{j_{1}},x_{j_{2}},...,x_{j_{k(j)}})$. Obviously, there are $k(j)$ ways to flip only one value in vector $\boldsymbol{v}$. Each value flip could potentially change the value of $f_{j}$. In some cases, $f_{j}$ can be very sensitive to its input vectors (e.g., the sensitivity of function $f=x_{j_{1}}\cdot x_{j_{2}}\cdot...\cdot x_{j_{k(j)}}$ to input vector (1,1,…,1)). To take this into consideration, for a given element $x_{j}$, we define its element sensitivity to a certain value $\mathbf{v}$ of its influence vector $\boldsymbol{v},s^{j}(\boldsymbol{v}=\mathbf{v})$ as the number of different ways individual elements of the influence set can be flipped within $\boldsymbol{v}$ to change the value of $x_{j}$. If we define $\boldsymbol{v}^{i}$ as a vector of size $k(j)-1$, same as vector $\boldsymbol{v}$ but missing element $x_{i}$, then we can compute $s^{j}(\boldsymbol{v}=\mathbf{v})$ as:

[TABLE]

We then define an overall element sensitivity of $x_{j}$, to all possible influence set value vectors, with respect to a given value distribution, as:

[TABLE]

As shown in Equation (4), $s^{j}$ has been proven equivalent to the sum of the influences of all its regulators in function $f_{j}$ (Shmulevich et al.,, 2004), where individual regulator influences can be computed according to Equation (2). We later call this procedure of computing element sensitivity for $x$ as ELEMENT_SENSITIVITY($x$).

We could see that the upper bound of $s^{j}$ is $k(j)$ as element influence $\alpha_{i}^{j}$ is bound by 1, and that the sensitivity of element $x_{j}$ is also dependent on the size of its influence set $k(j)$. For the same network as in Example 1, the element sensitivity of $x_{5}$ to input vector $(x_{1},x_{2},x_{3},x_{4})$ is given by $\alpha_{1}^{5}+\alpha_{2}^{5}+\alpha_{3}^{5}+\alpha_{4}^{5}$.

As shown in (Shmulevich et al.,, 2004) and (Anon,, 2007), if an element’s sensitivity is greater than 1, this would lead to instability of this element (i.e., it enables this element’s change to propagate out of control). Additionally, (Shmulevich et al.,, 2004) has shown that, even with the same element sensitivity, unbalanced influence distribution of its regulators can make elements behave more stable and robust than elements with balanced regulator influence distribution.

3.3 Static Analysis

In this section, we describe our approach to analyzing static sensitivity of model elements, with respect to the distribution of all possible model states. If we assume a uniform distribution of model states, the influence of element $x_{i}$ in the regulation of $x_{j}$ can be expressed as

[TABLE]

In general, static analysis provides a method to estimate the distribution of $\boldsymbol{v}$. In our algorithm shown in Section 3.5, we denote this calculation as EST_PROB($\boldsymbol{v},static$), where $\boldsymbol{v}$ is any vector and $static$ is an indicator of the method. As can be seen from Equation (4), the sensitivity of model elements to changes in values of other elements is determined by the fixed set of element update rules. In other words, the static sensitivity analysis approach relies solely on $M(V,F)$, therefore, it provides information only about the connectivity and logic rules of the model. Thus, for the model in Example 1, if we assume that $p_{1}=p_{2}=...=p_{8}=1/8$, $\alpha_{1}^{5}$ will be equal 3/8, which is solely determined by $\frac{\partial f_{5}}{\partial x_{1}}=x_{2}x_{3}^{{}^{\prime}}x_{4}^{{}^{\prime}}+x_{2}^{{}^{\prime}}x_{3}=1$. Under the same assumption, we can easily compute the other three influences for Example 1 as $\alpha_{2}^{5}=5/8,\alpha_{3}^{5}=3/8$, and $\alpha_{4}^{5}=1/8$.

However, the assumption that the states are uniformly distributed fails to capture the dynamics of changes in model states, and consequently, the details about transient behavior along trajectories between the initial and final model states.

3.4 Dynamic Analysis

In a real biological system, many combinations of element states may never occur under particular conditions, or, given the structure of the system, they are not even possible. Additionally, scientists are often focusing on the system’s transient response to interventions (e.g., drugs or treatments), since early response states could determine long-term outcomes (Miskov-Zivanov et al.,, 2013). Therefore, in conducting sensitivity analysis, we account for (i) uncertainty in the information available about the system, (ii) influence of the system’s initial state on its response, and (iii) stochasticity in system’s response to stimulations and interventions. In other words, we assume that (1) the available information may not be sufficient to derive the exact state distribution, and that (2) the distribution of states may vary, depending on initial states, interventions and perturbations.

To tackle the challenge (1) above, we estimate the distribution of states through simulations. We can simulate the model for a pre-determined number of steps, and we use element trajectories that we obtain through simulations to derive the state preference (i.e., which states are more likely to be reached). To tackle the challenge (2) above, we conduct simulations for all the initial states that are of interest when studying a given system.

As described in Section 2.3, DiSH simulator supports several different simulation schemes for discrete models. Here, we are interested in the random order update scheme (RSB-SQ), where simulation is run multiple times, on the same model and under the same scenario, and the trajectories obtained vary between different runs. Additionally, due to the randomness of trajectories, a particular state can be reached multiple times within a single run. Thus, we define the preference, $pe(\boldsymbol{v})$, of a certain state as a ratio of the number of trajectories $n(T_{\boldsymbol{v}})$ that reached this state at least once, to the total number of trajectories $n(T)$. Finally, we estimate the probability $p(\boldsymbol{v})$ of a certain state $\boldsymbol{v}$ via normalization of $pe(\boldsymbol{v})$, that is:

[TABLE]

[TABLE]

The dynamic trajectories obtained from simulation are often highly dependent on the initial state, and thus, the state preference results obtained from the sample trajectories, as well as our final estimated state distribution will be different for different initial states. In the algorithm (Fig. 5), we name the procedure of estimating state probabilities in the dynamic-based analysis as EST_PROB($\boldsymbol{v}$, $dynamic$, TRAJECT). This procedure takes as arguments an influence vector $\boldsymbol{v}=(x_{j_{1}},x_{j_{2}},...,x_{j_{k(j)}})$, indicator variable $dynamic$, and a trajectory TRAJECT, obtained from simulations. Furthermore, this procedure could be divided into two steps, COUNT($\boldsymbol{v}$, TRAJECT) which returns the count of trajectories that reach state $\boldsymbol{v}$ (at any simulation step, not only the final step), and NORMALIZE(COUNT), which normalizes the counts according to the trajectory set size and scales them to sum to 1.

In the dynamic sensitivity analysis approach, the influence of element $x_{i}$ in function $f_{j}$ is defined by considering the state distribution $p(\boldsymbol{v})$.

[TABLE]

For the model in Example 1, we denote the influence vector of $x_{5}$ as $\boldsymbol{v}=(x_{1},x_{2},x_{3},x_{4})$. For instance, the initial state of the network we are interested in is $(x_{1},x_{2},x_{3},x_{4})=(0,0,0,0)$, and RSB-SQ scheme will be applied. We run $n(T)=1500$ times to obtain 1500 trajectories, of which $n(T_{\boldsymbol{v}})=1050$ trajectories have reached state $(x_{1},x_{2},x_{3},x_{4})=(0,0,1,1)$. Then $pe(0,0,1,1)=1050/1500=0.7$. If, for example, we compute $\sum_{\boldsymbol{v}}pe(\boldsymbol{v})=7$, $p(0,0,1,1)$ could be obtained via normalization as $0.7/7=1/10$, no longer $1/8$ as static analysis. Element influences $\alpha_{1}^{5},\alpha_{2}^{5},\alpha_{3}^{5},\alpha_{4}^{5}$ will then be calculated based on these new joint probabilities.

As can be seen from Equation (8), $\alpha_{i}^{j}$ depends on both state dynamics along the trajectories and the topological structure of the model. Given these dependencies, it is possible that some influence vector values that have significant influence under static analysis, never occur under dynamic analysis, thus leading to $\alpha_{i}^{j}=0$ in the dynamic case.

3.5 Sensitivity Analysis with Weights

With our sensitivity analysis framework, we are now able to extend the model influence graph $G(V,E)$ to a weighted directed graph $G(V,E,W)$ by adding weights $w_{ij}=\alpha_{i}^{j}$ to the directed edges pointing to $x_{j}$ from $x_{i}$. Here $\alpha_{i}^{j}$ can be obtained either from static analysis, as shown in Equation (5) or from dynamic analysis, as shown in Equation (8). In Section 4.2, we will also discuss several different graph attributes, and define weights based on those attributes. While for a pre-defined discrete model, the interaction map $G(V,E)$ is fixed, the weighted directed graph $G(V,E,W)$ can vary with respect to $p(\boldsymbol{v})$ of state $\boldsymbol{v}$, which is determined by a simulation scenario. We denote the generation of this graph as GENERATION_WDG($M$, ELEMENT_INFLUENCE($x_{i},x_{j}$)), where $M$ is the model and ELEMENT_INFLUENCE($x_{i},x_{j}$)) provides the weights $w_{ij}=\alpha_{i}^{j}$. To illustrate this, we show in Fig. 4 an example graph $G(V,E,W)$ of the model from (Miskov-Zivanov et al.,, 2013), using RSB-SQ simulation scheme and high antigen dose scenario.

In Fig. 5, we provide the algorithm of our sensitivity analysis method, with standardized procedure of generating a weighted directed graph for a given model with Boolean update functions.

There are four main parts of the algorithm outlined in Fig. 5. Lines 2-6 include model initialization, which returns an influence graph $G(V,E)$ and an executable model $M(V,F)$. Lines 7-8 represent the simulation part: given any scenario of interest, $SCE_{m}$, and any available scheme, $SCH_{n}$, it returns a trajectory set, TRAJECT. Lines 9-22 serve as the main part of the algorithm to compute influence and sensitivity. Lines 23-24 generate a weighted directed graph.

4 Intervention Pathway Discovery

In this section, we discuss application and validation of our sensitivity analysis approach in the discovery of intervention pathways using graph search. Specifically, we propose a method for finding intervention pathways, while emphasizing the importance of element influence and sensitivity attributes in the discovery process.

When studying biological systems, finding the most influential pathways from one component (source) of a system to another (target) is critical in understanding the system, and in developing effective interventions or treatments. Once we identify these pathways in a model, we can easily control and guide the model by tuning its inputs, or by deciding which elements we should change during the transient process. In Section 4.1, we first reduce the problem of extracting pathways to a minimum cost search problem, based on a given definition of importance. Here, we selected to use the best-first search (Heuristics,, 1984) since it cooperates well with weighted graph and avoids unnecessary expansion via priority queue and hash table. In Section 4.2, we list nine candidate attributes that we will use in studying the model graph, and we define the weights that correspond to these attributes. To validate our proposed approach, in Section 4.4, we modify a well-studied biological model (Miskov-Zivanov et al.,, 2013) by removing the discovered important pathways, and then test the modified model to determine whether the model properties have been negatively affected. By testing whether the reduced model satisfies ground truths, we show that model properties are not preserved when we remove the pathways that our algorithm found.

4.1 Best-First Search

Given a weighted directed graph, a source node and a target node, our goal is to find the most ‘‘important’’ path(s) from source to target. The importance of a path could be interpreted in different ways. First, we could assume that the path having the highest cumulative element sensitivity is the most important, since the changes of regulators on this pathway will usually result in significant changes of regulated elements. Second, elements that have large influence sets usually serve as hubs in the network, and often form cycles (Albert,, 2005), therefore, we could also define the most important pathway as the one passing through elements with large influence sets. In other words, given a particular definition of pathway importance, the extraction of important pathways could be reduced to problems of maximizing/minimizing some index along pathways in a graph search.

That is, given graph $G(V,E,W)$, source node $x_{s}$ and target node $x_{t}$, we aim to find a path that has the minimum weight summation (also known as cost) along it, that is, we define optimal path as:

[TABLE]

where $w_{sp_{1}}$ is the weight on the edge between $x_{s}$ and $x_{p_{1}}$, $w_{p_{m}t}$ is the weight on the edge between $x_{p_{m}}$ and $x_{t}$, and $w_{p_{i}p_{i+1}}$ is the weight associated with the edge pointing to node with index $p_{i+1}$ from node with index $p_{i}$. In the following section, we will discuss several approaches for computing $w_{ij}$. For instance, the definition in Section 3.5 is one of our definitions of weight.

We apply the best-first search via maintaining a tree of paths originating at the source node and extending those paths one edge at a time until the target is reached. At each node, the search algorithm needs to determine which of its paths to extend based on the cost of the partial path. Specifically, best-first search selects the path that minimizes $g(x)$, where $x$ is the next node on the path, and $g(x)$ is defined as the cost of the path from the source node to node $x$. Our algorithm of best-first search uses a priority queue to perform the repeated selection of nodes with minimum $g$ value. This priority queue is known as the open set. At each step, the node with the lowest $g$ value is removed from the queue, the $g$ values of its successors are updated accordingly, and these successors are then added to the queue. The algorithm continues until the target node is popped up from the queue. The trace of elements along which we reach the target is the path with minimum cost, which is also the most important path based on our interpretation of importance.

We also introduce the following in our best-first search algorithm: (i) the flexibility to allow or prohibit cycles as paths with cycles within a biological system form feedback or feed-forward controls; (ii) the flexibility to return one or multiple minimum cost paths as two paths might have the same minimum cost, and in addition, multiple paths are useful for later evaluation which we will elaborate in Section 4.4.

4.2 Network Attributes

Different interpretations of importance will assign different weights to edges. The influence graph $G(V,E)$ and the discrete model $M(V,F)$ together carry a number of structural and dynamical features that can serve as useful attributes for detecting important pathways. We propose nine different attributes and define the edge weight according to the minimum cost graph search.

IN-DEGREE: $a_{1}(x)$, is the size of the $x$’s influence set. In the context of logical functions (and logic circuits), an element with large influence set tends to have a larger fan-in (the part of the circuit feeding into it), and dominates the circuit behavior. Accordingly, we aim to obtain path with more nodes with high in-degree (equivalently, low value in the reciprocal value of in-degree). Under such interpretation, we define the weight from node $x_{i}$ to node $x_{j}$ as $w_{ij}=\frac{1}{a_{1}(x_{j})}$.

OUT-DEGREE: $a_{2}(x)$, counts how many times an element $x$ occurs in other elements’ influence sets. Note that a self-loop regulation also counts. Generally, elements with higher out-degree also have a larger fan-out (the part of the circuit that this elements feeds into), and propagate the influence more broadly through the network. Thus, such elements often play hub roles in signal transduction. Similarly, we define $w_{ij}=\frac{1}{a_{2}(x_{j})}$.

SHORTEST_LINK($x$): $a_{3}(x)$, is defined as the shortest length path that goes through $x$ and connects given source and target elements. It is obvious that an element $x$ showing in a short path linking the source and target elements has a more significant impact in the relationship between source and target. Therefore, we could define $w_{ij}=a_{3}(x_{j})$.

LOOP_COUNT($x$): $a_{4}(x)$, counts how many loops go through $x$. In biology, feedback loops are common and often have crucial influence on behavior. A minor change in any element within a feedback loop could be amplified by multiple cycles through the loop. Thus, we expect that paths with elements that belong to multiple loops will be more influential, and therefore, we compute the weight as $w_{ij}=\frac{1}{a_{4}(x_{j})}$.

NON-BIAS($x$): $a_{5}(x)$, defined as $1-2\cdot|Pr\{x=1\}-0.5|$ for a Boolean variable node. It is obvious that $a_{5}(x)$ ranges from 0 to 1. If the state of an element is biased towards 0 or 1 (i.e., non-bias degree approaches 0), this element is robust against perturbations and in most cases, it prevents further signal propagation (Shmulevich et al.,, 2004, Anon,, 2007). We define $w_{ij}=-log(a_{5}(x_{j}))$ with the purpose of accumulating non-bias degree along the path.

ELEMENT_INFLUENCE$(x_{i},x_{j})$: $a_{6}(x_{i},x_{j})$, a result of Equations (5) and (8), this attribute is further categorized as $a_{6-static}$ and $a_{6-dynamic}$, which return the static/dynamic analysis of influence of $x_{i}$ in the regulation of $x_{j}$, respectively. We define $w_{ij}=-log(a_{6}(x_{i},x_{j}))$.

ELEMENT_SENSITIVITY($x$): $a_{7}(x)$, a result of Equation (4). Similar to $a_{6}$, this attribute could also be categorized into $a_{7-static}$ and $a_{7-dynamic}$, which return the static and dynamic analysis of element sensitivity of $x$, respectively. We define $w_{ij}=-log(a_{7}(x_{j}))$.

The attributes listed above could be divided into two types: structure-related and dynamics-related. The former category highlights the role of an element in network topological structure. Attributes $a_{1},a_{2},a_{3},a_{4},a_{6-static},a_{7-static}$ fall into this category. The latter category takes network dynamics into account, and it includes attributes $a_{5},a_{6-dynamic},a_{7-dynamic}$. We apply the best-first search based on the above definitions of weights, and investigate in Section 5.2 which of these weight types lead to a better extraction of important pathways.

4.3 Pathway Score

With the algorithm explained in Section 4.1 and weights definition in Section 4.2, we are now able to extract intervention pathways via graph search. Assuming that there is a regulatory pathway between a source node $x_{s}$ and a target node $x_{t}$, we assign a pathway score to path:$\{x_{s},x_{p_{1}},...,x_{p_{m}},x_{t}\}$ such that

[TABLE]

The pathway score gives the influence of source node on the target node following that pathway. The optimal path (i.e., with minimum cost) has the highest path score. It is also clear that, with weight type $a_{6}(x_{i},x_{j})$ (i.e.,$w_{ij}=-log(a_{6}(x_{i},x_{j}))$), computing $s(path)$ will return the element influence product. If we assume the independence between regulatory edges on a pathway, this product will indicate the propagated probability effect along the pathway, thus providing a measure of control from source node to target node.

4.4 Pathway Extraction Evaluation

The underlying idea of our pathway extraction validation method is that the biological properties of the original model will not be preserved when an important pathway in the model (i.e., all of the interactions within this pathway) has been removed, while the removal of a pathway with low influence should not change the properties significantly. We use the biological property mismatch between the partial model after removal and the original baseline model, to verify our pathway extraction approach.

To conduct validation, we follow the algorithm described in Section 4.1 and the weight definitions in Section 4.2. For each source-target pair of interest, and for each weight type, we define the pathway extraction threshold, that is, the number of top paths that we would like to extract and explore, assuming paths are ranked according to their impact. During the extraction, in case we encounter cycles, we allow at most one repeat for some pre-defined (e.g., according to prior knowledge) list of elements, and prohibit any other repeated element occurrences on the paths. After obtaining top paths for each weight definition, we use the following procedure to validate the results of pathway extraction. First, we define a set of properties that are true for the modeled system. For example, a property can state a gene expression level at a particular time during observation interval, or relative changes in proteins in time. Next, we remove all edges (regulations) in the obtained top pathways from the original model to form a partial model. Finally, we use the statistical model checking approach (Miskov-Zivanov et al.,, 2016, Wang et al.,, 2016), combined with DiSH simulations, to test all partial models on the defined set of properties, which returns the probabilities for each model satisfying each property.

In particular, when we remove the regulation of element $x_{j}$ by element $x_{i}$, we modify $x_{j}$’s update function such that $x_{i}$ is removed while keeping other elements and their original effects in the update function unchanged, which is, in some cases, achieved by setting $x_{i}$ to 0 or 1 permanently in the update function. In Example 1, knocking down the regulation of element $x_{1}$ in $x_{5}=f_{5}=x_{1}x_{2}+x_{1}^{{}^{\prime}}x_{3}+x_{2}x_{3}^{{}^{\prime}}x_{4}$ returns $x_{5}=f_{5}=x_{2}+x_{3}+x_{2}x_{3}^{{}^{\prime}}x_{4}$. Note that this modification does not indicate the deletion of $x_{1}$ from our research scope, nor does it change the effect of $x_{1}$ on other elements than $x_{5}$.

5 Results

In Section 5.1, we provide a brief background on T cell differentiation and the baseline model. In Section 5.2, we confirm, using model checking and the T cell differentiation case study system, that our sensitivity analysis framework is successful in intervention pathway discovery. Finally, we study various features of the T cell differentiation circuitry using our proposed sensitivity analysis methods.

5.1 Case Study: T-cell Differentiation

T cells, one of two primary types of lymphocytes, play a central role in cell-mediated immunity (Gutcher et al.,, 2007). There are several subsets of T cells, and each one has a distinct function in the T-cell mediated immunity. In this work, we investigate T cell differentiation as our case study. Specifically, we use the model of the circuitry that controls differentiation of naive T cells into two types, regulatory cells (Treg), which suppress immune response and reduce the damage caused by autoimmune response, and helper cells (Th), which help promote immune response. Previous research (Gutcher et al.,, 2007) has shown that these two types of T-cells are distinguished by expressions of several key elements. For example, in the Treg type, the transcription factor forkhead box P3 (Foxp3) is expressed, and Interleukin-2 (IL-2) is inhibited, while in the Th type, Foxp3 is inhibited and IL-2 is activated.

The model in (Miskov-Zivanov et al.,, 2013) has overall 42 elements, including receptors TCR, CD28, TGF-$\beta$R, and IL-2R, genes for IL-2R$\alpha$ (i.e., CD25), IL-2, and Foxp3, and cytoplasmic members of PI3K/AKT/mTOR, MAPK, NF$\kappa$B, NFAT pathways. The authors in (Miskov-Zivanov et al.,, 2013) used logical modeling approach to study naive T cell differentiation. In that work, several model elements are implemented as three-level discrete variables with values {0,1,2}, to denote absence, low activity, and high activity of the element, respectively. To accommodate three discrete levels within a logical model, two Boolean variables are associated with these elements with three levels. For example, the T-cell receptor (TCR) is modeled with two variables, TCR_LOW and TCR_HIGH, such that: TCR=0: (TCR_LOW=0, TCR_HIGH=0); TCR=1: (TCR_LOW=1, TCR_HIGH=0); TCR=2: (TCR_LOW=0, TCR_HIGH=1). The authors in (Miskov-Zivanov et al.,, 2013) also use extra, dummy variables, to incorporate delays on several pathways, and therefore, the model in (Miskov-Zivanov et al.,, 2013) has overall 55 logic variables. Here, we will focus on exploring actual element sensitivities, and will thus use element names instead of variable names. Further details of the variable implementation can be found in (Miskov-Zivanov et al.,, 2013).

Besides the list of elements, and their update functions, as described in Section 3.4, we also need to define the scenarios for conducting dynamic sensitivity analysis on the T cell model. Similar to the studies in (Miskov-Zivanov et al.,, 2013), we explore three dynamic scenarios: Scenario 1: high antigen dose (TCR=2), Scenario 2: low antigen dose scenario (TCR=1), and Scenario 3: a toggle scenario (TCR=2 initially, and changed to TCR=0 at a defined time step).

5.2 Pathway Discovery

To demonstrate the use of our pathway discovery method, we start with the model from (Miskov-Zivanov et al.,, 2013) and select two source-target pairs, (TCR, Foxp3) and (CD28, IL-2). Following the procedure described in Section 4.4, for each of the nine network attributes defined in Section 4.2, and for each source target pair, we obtain three paths with minimum costs, computed according to Equation (9). We then create new models by removing these paths from the original model, and test their performance using model checking. While, in general, more complicated properties can be used, for our analysis of the T cell model, we test whether steady-state values of IL-2, Foxp3, AKT, and PTEN are reached. Specifically, we test whether elements (IL-2, Foxp3, AKT, PTEN) reach values (1,0,1,0) in Scenario 1, and (0,1,0,1) in Scenarios 2 and 3. We call the probability of satisfying a property, obtained from model checking, property match probability.

For each model, we combine the property match probabilities into average match probability by computing an average of these four probabilities. In Fig. 6(a), we show the results for average match probability for several model versions under the three scenarios defined in Section 5.1. The models that we compare are the original model (Miskov-Zivanov et al.,, 2013) and partial models after removing pathways based on using different network attributes. As can be seen from Fig. 6(a), in the original model (solid line with triangle marker), properties are highly satisfied (i.e. average value of these four probabilities is close to 1). Note that we illustrate the results for the $a_{1}$-based and $a_{2}$-based removal using a single line in Fig. 6(a), since their results are the same; similarly, we use a single line for $a_{7-static}$-based and $a_{5}$-based removal. It can be observed that all removals affect the property match probability, and moreover, while removing pathways based on in-degree or out-degree distribution ($a_{1}$ and $a_{2}$), minimum length ($a_{3}$), and loop count ($a_{4}$) affects the properties only to some extent under Scenarios 1 and 3, removing pathways based on element influence ($a_{6-static}$, $a_{6-dynamic}$) significantly breaks down the model (match probability approaches 0) under Scenarios 1 and 2.

In addition, the dynamic sensitivity analysis of the $a_{6-dynamic}$-based pathway extraction is conducted with respect to the distribution estimated under Scenario 2. Therefore, we explored further this particular scenario, and illustrated in Fig. 6(b) these four detailed probabilities under Scenario 2. It is clear from the figure that removing pathways based on the dynamic influence value $a_{6-dynamic}$ calculated under Scenario 2, has an even more destructive effect on the model behavior under that scenario. By removing the pathways based on $a_{6-dynamic}$, three of the four key elements have a significant mismatch with their corresponding desired properties (solid line with triangle marker in the figure), thus revealing the critical role of the removed pathways. To further show the effect, we plot in Fig. 7 the change in time (trajectory averaged across 50 runs, each run with 1500 steps) of element IL-2 in the original model (left) and in the partial model (right) after we removed pathways using the $a_{6-dynamic}$ attribute.

5.3 Element-level Analysis: Element Sensitivity

By applying model checking and testing important system properties, we have shown in Section 5.2 that element influence and element sensitivity are good indicators of pathway importance. This is due to the fact that, the calculations of element influence and sensitivity in Equations (4, 5, 8) take both network structure and network dynamics into account. In other words, element influence and sensitivity can be used to define weights and incorporate context-dependent information when determining important pathways, thus enabling extension from local to global analysis.

For the T cell model from (Miskov-Zivanov et al.,, 2013), the interaction map $G(V,E)$ is fixed, while the weighted directed graph $G(V,E,W)$ varies with respect to state distribution, which is determined by simulation scenario. Under different scenarios, elements have different sensitivities. Fig. 8 shows the sensitivity of model elements using static analysis defined by Equation (5) in Section 3.3 (blue bars), as well as dynamic analysis defined by Equation (8) in Section 3.4 under the three different scenarios (red, orange and purple bars for Scenario 1,2, and 3, respectively). Note that in Fig. 8 we omit results for elements that do not have update functions (CD28, TCR, TGF$\beta$, TGF$\beta$R, CD122, CD132), and therefore, cannot be sensitive to changes in any other model elements.

As can be seen from Fig. 8, static analysis results in less variance, ranging from 0 to 1.5. Most elements have sensitivities less than 1, indicating that this network is likely to follow a stable and ordered structured behavior (Shmulevich et al.,, 2004). However, dynamic analysis (in all three scenarios) shows greater variance, ranging from 0 to 2.94. Some elements under dynamic analysis have sensitivities much greater than 1 (e.g. element PKC$\theta$ in low-dose scenario). These results indicate that PKC$\theta$ is almost 100% sensitive to any one-bit change of its regulators (i.e., TCR, CD28, mTORC2) and it will always propagate the change to its immediate downstream connections (i.e., NF$\kappa$B, TAK1). This kind of local ‘‘instability’’ can suggest potential intervention points on the pathways leading to phenotype markers Foxp3 and IL-2. Other elements which also show similar behavior include AKT, IL-2R in high-dose scenario, IL-2R, mTORC2, PIP3 in Scenario 2 (low-dose), and CD25, IL-2-exogenous (IL2-EX) in Scenario 3 (toggle).

Additionally, it can be observed from the results in Fig. 8 that the top elements behaving differently across the three scenarios are AKT and PKC$\theta$. For example, AKT shows robustness against changes in its regulators PDK1 and mTORC2 in Scenario 3 (toggle), while exhibits sensitivity to these regulators in Scenario 1 (high-dose). In other words, high-dose stimulation makes element AKT inclined towards the states that are highly sensitive and flexible against environment changes. This feature is quite useful when designing interventions to efficiently induce certain system response given particular initial state, or for identifying elements with critical role in the context-dependent analysis. Finally, this also emphasizes the importance of the context-dependent analysis itself.

5.4 Interaction-level Analysis: Element Influence

When designing models of systems on the basis of incomplete information, and moreover, when defining intervention strategies for such systems, it is beneficial to also conduct analysis at an element-to-element (interaction) level. This type of analysis can summarize the details of indirect regulations into a single measure, element-to-element influence. Using this measure, we can study influence of an arbitrary source element, $x_{s}$, on an arbitrary target element, $x_{t}$. First, we find all paths from $x_{s}$ to $x_{t}$, and then we summarize the score of all the paths as element-to-element influence. As defined in Equation (10), a single pathway score only implies the influence following that pathway. If we add the scores of all possible paths from $x_{s}$ to $x_{t}$, this summation shows the overall influence that the source element $x_{s}$ has on the target element $x_{t}$. In Fig. 9, we show the element-to-element influence matrices under four cases. Rows in the heatmaps correspond to source elements and columns correspond to target elements.

In general, element-to-element influences in three dynamic scenarios are greater than element-to-element influences under static analysis. There are two main reasons for this. First, under static analysis, the long-run element-to-element influences are quite sensitive to the length of pathways since the direct element influence $\alpha$ under static analysis is relatively small. Second, dynamic scenarios are obtained from real biological observations and thus show stronger homogeneity.

Additionally, it is worthwhile to note that, compared to the static analysis, the influence matrices of the three dynamic scenarios shows significant imbalance (i.e., with wide variance range and high deviation), especially the matrix of toggle case. It has been shown before that systems with unbalanced influence distributions are stable and robust (Anon,, 2007). Thus, the results in Fig. 9 suggest that we may toggle node TCR in the initial state, or toggle nodes like PKC$\theta$ and PI3K during the transient process in order to drive the system to certain states or to maintain system stability.

Intuitively, a high element-to-element influence (darker blue blocks) indicates strong interaction. It is easy to observe examples of such strong interactions from Fig. 9, such as CD28 to JUN, MKK7 in low-dose scenario, and JAK3, STAT5 to CD25, IL-2R in toggle scenario.

5.5 Static vs. Dynamic Analysis

Scenario specifics (i.e., context) are often critical for guiding intervention decisions. For example, detailed information about the effects of a scenario on the model helps evaluate the trade-off between accuracy and effort when designing context-specific vs. generalized interventions.

We have observed that some pathways that appear very influential in the static analysis are actually inactive in dynamic scenarios. These observations suggest that adding interventions following only static pathway and element function analysis may not be effective.

To demonstrate the significance of change in the number of active pathways when switching between static and dynamic analysis, quite common for many source/target combinations, Table 1 gives the comparison in number of active pathways under static analysis and dynamic analysis with high-dose scenario and low-dose scenario.

As can be seen from Table 1, accounting for the real scenarios that can occur in the cell (scenarios within our dynamic analysis) leads to reduction in the number of active pathways. More interestingly, we find that there is only one active pathway from IL-2 EX to IL-2 in high-dose/low-dose scenarios, which is IL-2 EX, JAK3, STAT5, Foxp3, IL-2. Also, in dynamic case, the number of active pathways of regulation for Foxp3 and IL-2 does not change with the change of a scenario, and this holds for all source nodes.

6 Conclusion

In this work, we model and investigate intra-cellular networks via discrete modeling approach, and we propose a framework to study sensitivity in these models. Previous sensitivity analysis assumes uniform state distribution, which is not always true in biology. We perform both static and dynamic sensitivity analysis, the former assuming uniform state distribution, and the latter using a distribution estimated from stochastic simulation trajectories for a given scenario. Within our sensitivity analysis framework, we first compute element-to-element influences, then we extend the element update functions to include weights according to these computed influences. Adding weights to these interaction rules helps to identify key elements in the model, as well as dominant signaling pathways that determine the behavior of the overall model. In the end, we apply our sensitivity analysis framework on the intervention pathway extraction and evaluation in the intra-cellular network that controls T cells differentiation. Our results emphasize the importance of incorporating context in sensitivity analysis and in the selection of intervention pathways.

Funding

This work is supported in part by DARPA award W911NF-17-1-0135 and W911NF-18-1-0017.