∗To whom correspondence should be addressed. Email<EMAIL_ADDRESS>
# Intervention Pathway Discovery via Context-Dependent Dynamic Sensitivity
Analysis

Gaoxiang Zhou ${}^{\text{\sfb 1,}*}$ Kai-Wen Liang ${}^{\text{\sfb 2}}$ and
Natasa Miskov-Zivanov ${}^{\text{\sfb 1}}$ ${}^{\text{\sf 1}}$Department of
Electrical and Computer Engineering, University of Pittsburgh, Pittsburgh,
15260, United States and
${}^{\text{\sf 2}}$Department of Electrical and Computer Engineering, Carnegie
Mellon University, Pittsburgh, 15213, United States

###### 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.

##  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.

Figure 1: Discrete modeling approach and simulation: (a) A toy model $M(V,F)$
with model elements $V:\\{x_{1},x_{2},x_{3}\\}$ and their update functions; ‘‘
$\cdot$ ’’ denotes AND operator, where both regulators are required to be
active (equal 1) to activate regulated element; ‘‘+’’ denotes OR operator,
where either regulator being active is sufficient to activate regulated
element; ‘‘ ′ ’’ denotes NOT operator, where inactivation of regulator
activates regulated element. (b) The state transition graph of the toy model
under SMLN (red arrows) and RSB-SQ schemes (black arrows, labeled with updated
element name).

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.

Figure 2: The flow diagram of our sensitivity analysis framework.

###  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):

$\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})$ (1)

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}$:

$\alpha_{i}^{j}=\alpha_{i}^{f_{j}}=E(\frac{\partial f_{j}}{\partial x_{i}})$
(2)

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)).

Figure 3: Element influence calculation of function
$x_{5}=f_{5}=x_{1}x_{2}+x_{1}^{{}^{\prime}}x_{3}+x_{2}x_{3}^{{}^{\prime}}x_{4}$.

###  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:

$s^{j}(\boldsymbol{v}=\mathbf{v})=\sum_{i=1}^{k(j)}\frac{\partial
f_{j}}{\partial x_{i}}(\boldsymbol{v}^{i}=\mathbf{v}^{i})$ (3)

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:

$s^{j}=E(s^{j}(\boldsymbol{v}=\mathbf{v}))=\sum_{i=1}^{k(j)}E(\frac{\partial
f_{j}}{\partial
x_{i}}(\boldsymbol{v}^{i}=\mathbf{v}^{i}))=\sum_{i=1}^{k(j)}\alpha_{i}^{j}$
(4)

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

$\alpha_{i}^{j}=\alpha_{i}^{f_{j}}=E(\frac{\partial f_{j}}{\partial
x_{i}})=\frac{1}{2^{k(j)}}\sum_{\boldsymbol{v}}\frac{\partial
f_{j}(\boldsymbol{v})}{\partial x_{i}}$ (5)

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:

$pe(\boldsymbol{v})=\frac{n(T_{\boldsymbol{v}})}{n(T)}$ (6)
$p(\boldsymbol{v})=\frac{pe(\boldsymbol{v})}{\sum_{\boldsymbol{v}}pe(\boldsymbol{v})}$
(7)

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})$.

$\alpha_{i}^{j}=\alpha_{i}^{f_{j}}=E(\frac{\partial f_{j}}{\partial
x_{i}})=\sum_{\boldsymbol{v}}\frac{\partial f_{j}(\boldsymbol{v})}{\partial
x_{i}}p(\boldsymbol{v})$ (8)

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.

Figure 4: A weighted directed graph of a biological model (Miskov-Zivanov et
al.,, 2013) under certain 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.

Figure 5: An overall algorithm of our sensitivity analysis framework

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:

$\\{x_{s},x_{p_{1}},...,x_{p_{m}},x_{t}\\}\triangleq\operatorname*{arg\,min}_{x_{s},...,x_{t}}w_{sp_{1}}+\sum_{i=1}^{m-1}w_{p_{i}p_{i+1}}+w_{p_{m}t}$
(9)

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

$s(path)=exp(-[w_{sp_{1}}+\sum_{i=1}^{m-1}w_{p_{i}p_{i+1}}+w_{p_{m}t}])$ (10)

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.

Figure 6: Probabilities of model property match for the original model and
several partial models: (a) under three different scenarios; (b) under
Scenario 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.

Figure 7: The activity of element IL-2: original model (left), and partial
model after pathway removal based on $a_{6-dynamic}$ (right).

###  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.

Figure 8: The sensitivity distribution of all elements under four cases of
analysis.

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.

Figure 9: The element-to-element influence matrix under four cases: static
analysis (top left), and the three dynamic analysis scenarios (top right and
bottom left and right)

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.

Table 1: Comparison in number of active pathways in different cases

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.

## References

  * Shmulevich et al., (2004) Shmulevich, I., & Kauffman, S. A. (2004). Activities and Sensitivities in Boolean Network Models. Physical Review Letters, 93(4). doi:10.1103/ physrevlett .93.048701
  * Faeder et al., (2009) Faeder, J. R., Blinov, M. L., & Hlavacek, W. S. (2009). Rule-Based Modeling of Biochemical Systems with BioNetGen. Methods in Molecular Biology Systems Biology, 113-167. doi:10.1007/978-1-59745-525-1_5
  * Albert et al., (2014) Albert, R., & Thakar, J. (2014). Boolean modeling: A logic-based dynamic approach for understanding signaling and regulatory networks and for making useful predictions. Wiley Interdisciplinary Reviews: Systems Biology and Medicine, 6(5), 353-369. doi:10.1002/wsbm.1273
  * Miskov-Zivanov et al., (2013) Miskov-Zivanov, N., Marculescu, D., & Faeder, J. R. (2013). Dynamic behavior of cell signaling networks. Proceedings of the 50th Annual Design Automation Conference on DAC 13. doi: 10.1145/ 2463209. 2488743
  * Sun et al., (2014) Sun, Z., Jin, X., Albert, R., & Assmann, S. M. (2014). Multi-level modeling of light-induced stomatal opening offers new insights into its regulation by drought. PLoS computational biology, 10(11), e1003930
  * Gan et al., (2016) Gan, X., & Albert, R. (2016). Analysis of a dynamic model of guard cell signaling reveals the stability of signal propagation. BMC Systems Biology, 10(1). doi:10.1186/s12918-016-0327-7
  * Shmulevich et al., (2002) Shmulevich, I., Dougherty, E. R., Kim, S., & Zhang, W. (2002). Probabilistic Boolean networks: A rule-based uncertainty model for gene regulatory networks. Bioinformatics, 18(2), 261-274. doi:10.1093/bioinformatics/18.2.261
  * Shmulevich et al., (2002) Shmulevich, I., Dougherty, E., & Zhang, W. (2002). From Boolean to probabilistic Boolean networks as models of genetic regulatory networks. Proceedings of the IEEE, 90(11), 1778-1792. doi:10.1109/jproc.2002.804686
  * Anon, (2007) Anon. (2007). Analysis of random Boolean networks using the average sensitivity. (April 2007). Retrieved June 15, 2018 from https://archive.org/details/arxiv-0704.0197
  * Qian et al., (2009) Qian, X., & Dougherty, E. R. (2009). On the long-run sensitivity of probabilistic Boolean networks. Journal of Theoretical Biology, 257(4), 560-577. doi:10.1016/j.jtbi.2008.12.023
  * Trinh et al., (2016) Trinh, H., & Kwon, Y. (2016). Edge-based sensitivity analysis of signaling networks by using Boolean dynamics. Bioinformatics, 32(17), I763-I771. doi:10.1093/bioinformatics/btw464
  * Kochi et al., (2014) Kochi, N., Helikar, T., Allen, L., Rogers, J. A., Wang, Z., & Matache, M. T. (2014). Sensitivity analysis of biological Boolean networks using information fusion based on nonadditive set functions. BMC Systems Biology, 8(1). doi:10.1186/s12918-014-0092-4
  * Miskov-Zivanov et al., (2016) Miskov-Zivanov, N., Zuliani, P., Wang, Q., Clarke, E. M., & Faeder, J. R. (2016). High-level modeling and verification of cellular signaling. 2016 IEEE International High Level Design Validation and Test Workshop (HLDVT). doi:10.1109/hldvt.2016.7748271
  * Wang et al., (2016) Wang, Q., Miskov-Zivanov, N., Liu, B., Faeder, J. R., Lotze, M., & Clarke, E. M. (2016). Formal Modeling and Analysis of Pancreatic Cancer Microenvironment. Computational Methods in Systems Biology Lecture Notes in Computer Science, 289-305. doi:10.1007/978-3-319-45177-0_18
  * Miskov-Zivanov et al., (2013) Miskov-Zivanov, N., Turner, M. S., Kane, L. P., Morel, P. A., & Faeder, J. R. (2013). The Duration of T Cell Stimulation Is a Critical Determinant of Cell Fate and Plasticity. Science Signaling, 6(300). doi:10.1126/scisignal.2004217
  * Kaderali et al., (2008) Kaderali, L., & Radde, N. (2008). Inferring gene regulatory networks from expression data. In Computational Intelligence in Bioinformatics (pp. 33-74). Springer, Berlin, Heidelberg
  * Mihajlovic et al., (2001) Mihajlovic, V., & Petkovic, M. (2001). Dynamic bayesian networks: A state of the art. University of Twente Document Repository
  * Pe’er et al., (2001) Pe’er, D., Regev, A., Elidan, G., & Friedman, N. (2001). Inferring subnetworks from perturbed expression profiles. Bioinformatics, 17(suppl_1), S215-S224
  * Davidich et al., (2008) Davidich, M. I., & Bornholdt, S. (2008). Boolean network model predicts cell cycle sequence of fission yeast. PloS one, 3(2), e1672
  * Friedman et al., (2000) Friedman, N., Linial, M., Nachman, I., & Pe’er, D. (2000). Using Bayesian networks to analyze expression data. Journal of computational biology, 7(3-4), 601-620
  * Sayed et al., (2017) Sayed, K., Telmer, C. A., Butchy, A. A., & Miskov-Zivanov, N. (2017, September). Recipes for Translating Big Data Machine Reading to Executable Cellular Signaling Models. In International Workshop on Machine Learning, Optimization, and Big Data (pp. 1-15). Springer, Cham
  * Wang et al., (2016) Wang, Q., Miskov-Zivanov, N., Liu, B., Faeder, J. R., Lotze, M., & Clarke, E. M. (2016, September). Formal modeling and analysis of pancreatic cancer microenvironment. In International Conference on Computational Methods in Systems Biology (pp. 289-305). Springer, Cham
  * Sayed et al., (2017) Sayed, K., Kuo, Y. H., Kulkarni, A., & Miskov-Zivanov, N. (2017, December). DiSH simulator: Capturing dynamics of cellular signaling with heterogeneous knowledge. In Simulation Conference (WSC), 2017 Winter (pp. 896-907). IEEE
  * Heuristics, (1984) Heuristics, J. P. (1984). Intelligent Search Strategies for Computer Problem Solving Addison.
  * Albert, (2005) Albert, R. (2005). Scale-free networks in cell biology. Journal of cell science, 118(21), 4947-4957.
  * Gutcher et al., (2007) Gutcher, I., & Becher, B. (2007). APC-derived cytokines and T cell polarization in autoimmune inflammation. The Journal of clinical investigation, 117(5), 1119-1127.