Combinatorial Algorithms for Control of Biological Regulatory Networks

TL;DR

This paper addresses the challenge of controlling biological regulatory networks by identifying minimal control nodes, proposing algorithms for specific network structures, and validating methods on real gene networks.

Contribution

It introduces new algorithms for control node selection in Boolean biological networks, exploiting structural properties for computational efficiency.

Findings

Control node set can be computed efficiently for certain network topologies.

Polynomial-time algorithms are developed for networks with nested canalyzing dynamics.

Validation on real gene networks demonstrates practical effectiveness.

Abstract

Biological processes, including cell differentiation, organism development, and disease progression, can be interpreted as attractors (fixed points or limit cycles) of an underlying networked dynamical system. In this paper, we study the problem of computing a minimum-size subset of control nodes that can be used to steer a given biological network towards a desired attractor, when the networked system has Boolean dynamics. We first prove that this problem cannot be approximated to any nontrivial factor unless P=NP. We then formulate a sufficient condition and prove that the sufficient condition is equivalent to a target set selection problem, which can be solved using integer linear programming. Furthermore, we show that structural properties of biological networks can be exploited to reduce the computational complexity. We prove that when the network nodes have threshold dynamics and certain topological structures, such as block cactus topology and hierarchical organization, the input selection problem can be solved or approximated in polynomial time. For networks with nested canalyzing dynamics, we propose polynomial-time algorithms that are within a polylogarithmic bound of the global optimum. We validate our approach through numerical study on real-world gene regulatory networks.

Tables

Table 1. TABLE I: Number of Inputs for Control of Biological Networks

Network	Number of Nodes (Edges)	Number of Inputs
Apoptosis	39 (70)	10
Bordatella Bronchiseptica	33 (79)	2
Breast Cell Development Network	21 (70)	7
Mammalian Cell Cycle	19 (48)	3
T-Cell Differentiation	19 (30)	10
T-Cell Signaling	37 (48)	3

Equations

$$x_i(t+1)=f_i(\mathbf{x}(t)),$$

$$f_{i}(\mathbf{x})=\left\{\begin{array}[]{ll}1,&\sum_{j\in N_{in}(i)}{a_{ij}x_{j}}\geq\tau_{i}\\ 0,&\mbox{else}\end{array}\right.$$

$$f_{i}(\mathbf{x})=\left\{\begin{array}[]{ll}a_{l},&\mbox{if }x_{j_{1}}\neq b_{1},\ldots,x_{j_{l-1}}\neq b_{l-1},x_{jl}=b_{l}\\ a_{d+1},&\mbox{else}\end{array}\right.$$

$$Y[k]=\left\{v\notin X[k-1]:|N_{in}(v)\cap X[k-1]|\ge \tau (v)\right\}$$

$$\begin{array}[]{ll}\mbox{minimize}&|S|\\ \mbox{s.t.}&S\in\mathcal{C}\end{array}$$

$$(p_{11}\vee \cdots \vee p_{1r_1})\wedge \cdots \wedge (p_{l1}\vee \cdots \vee p_{l{r}_l}),$$

$$x_{v_1^1}(t+1)x_{v_2^1}(t+1)\cdots x_{v_m^1}(t+1)=x_{v_1^1}(t)\cdots x_{v_m^1}(t)+1(\mathrm{mod}{2}^m).$$

$$x_{v_i^2}(t+1)=\underset{s=1}{\overset{r_i}{\bigvee }}{y}_{is}(t),$$

$$f_i(\mathbf{x})=(y_{11}\vee \cdots \vee y_{1r_{1i}})\wedge \cdots \wedge (y_{l_{i}1}\vee \cdots \vee y_{l_i{r}_{li}}),$$

$$\hat{V}=V\cup \{a_{i_s}:s=1,\dots ,l\}.$$

$$\hat{E}$$

$$\tau(a_{i_{m}})=\left\{\begin{array}[]{ll}r_{mi},&x_{i}^{\ast}=0\\ 1,&x_{i}^{\ast}=1\end{array}\right.\quad\tau(i)=\left\{\begin{array}[]{ll}1,&x_{i}^{\ast}=0\\ l_{i},&x_{i}^{\ast}=1\end{array}\right.$$

$$\begin{array}[]{ll}\mbox{minimize}&\sum_{i=1}^{n}{s_{i}}\\ \mathbf{s},\mathbf{e}&\\ \mbox{s.t.}&\sum_{(i,j)\in\hat{E}}{e_{ij}}\geq\overline{\tau}_{i}(1-s_{i})\ \forall i\in\hat{V}\\ &e_{ij}+e_{ji}=1\ \forall i\neq j\\ &e_{ij}+e_{jl}+e_{li}\leq 2\ \forall i,j,l\mbox{ distinct}\\ &e_{ij}\in\{0,1\},s_{i}\in\{0,1\}\end{array}$$

$$\hat{\tau}(i)=\left\{\begin{array}[]{ll}\tau(i)+|N_{in}(i)\cap\mathcal{I}|,&x_{i}^{\ast}=1\\ \tau(i)+|N_{in}(i)\cap\mathcal{E}|,&x_{i}^{\ast}=0\end{array}\right.$$

$$|{\hat{N}}_{in}(i)\cap X[k-1]|\ge \tau (i)+|N_{in}(i)\cap \mathcal{I}|,$$

$$\begin{gathered} |{\hat{N}}_{in}(i)\cap {\mathcal{E}}^1\cap X[k-1]|+|{\hat{N}}_{in}(i)\cap {\mathcal{E}}^0\cap X[k-1]| \\ +|{\hat{N}}_{in}(i)\cap {\mathcal{I}}^1\cap X[k-1]|+|{\hat{N}}_{in}(i)\cap {\mathcal{I}}^0\cap X[k-1]|\ge \tau (i)+|N_{in}(i)\cap \mathcal{I}|. \end{gathered}$$

$$|N_{in}(i)\cap {\mathcal{E}}^1\cap X[k-1]|+|N_{in}(i)\cap {\mathcal{I}}^0\cap X[k-1]|\ge \tau (i)+|N_{in}(i)\cap \mathcal{I}|.$$

$$\sum _{j\in N_{in}(i)}{a}_{ij}{x}_j(t)$$

$$\begin{gathered} \hat{E}=\{(i,j):i\in {\mathcal{E}}^1\cup {\mathcal{I}}^0,j\in {\mathcal{E}}^1\cup {\mathcal{I}}^1\} \\ \cup \{(i,j):i\in {\mathcal{E}}^0\cup {\mathcal{I}}^1,j\in {\mathcal{E}}^0\cup {\mathcal{I}}^1\} \end{gathered}$$

$$|{\mathcal{E}}^1\cap \overline{X}[k-1]|+|{\mathcal{I}}^0\cap \overline{X}[k-1]|\ge \overline{\tau }(i).$$

$${\overline{\alpha }}_k$$

$${\alpha }_k$$

$${\overline{\beta }}_k$$

$${\beta }_k$$

$${\mathcal{E}}^1\cap \overline{X}[k-1]=\{s\in {\mathcal{E}}^1:\overline{\tau }(s)<{\overline{\alpha }}_{k-1}\}\cup (S\cap {\mathcal{E}}^1),$$

$$\{s\in {\mathcal{E}}^1:\overline{\tau }(s)<{\alpha }_{k-1}\}\subseteq \{s\in {\mathcal{E}}^1:\overline{\tau }(s)<{\overline{\alpha }}_{k-1}\}.$$

$$|(({\mathcal{E}}^1\cap X[k-1])\setminus ({\mathcal{E}}^1\cap \overline{X}[k-1]))|\le |(({\mathcal{E}}^1\cap \overline{X}[k-1])\setminus ({\mathcal{E}}^1\cap X[k-1]))|$$

$$\frac{R(k,m)}{n}=\frac{(2k{)}^m}{(k+1{)}^m}\le 2^m=n,$$

$$|S|\le \sum _{i=1}^l|S\cap V_i|+1\le \sum _{i=1}^mϵ(m)|S_i^{\ast }|+1.$$

$$\frac{|S|}{|S^{\ast }|}$$

$$\underset{j_s\in {\Omega }_i}{\bigvee }\left(\left(\underset{j_r\in {\Omega }_i^c\cap \{1,\dots ,s\}}{\bigwedge }(\neg x_{a_r})\right)\wedge x_{j_s}\right).$$

$$\frac{|S|}{|S^{\ast }|}\le \sum _{i=1}^m\frac{|S_i|}{|S_i^{\ast }|}\le (\log n{)}^2,$$

Full text

I Introduction

Biological processes, including gene expression and metabolism, are driven by complex interactions between basic building blocks. These interactions are often modeled as networked dynamical processes, in which nodes represent genes or proteins, links represent regulation of one component by another, and the node states describe the level of expression of each gene or protein. One widely-studied modeling approach is to assign a Boolean (on or off) state to each node, while describing the state of a node at each time step as a Boolean function of the neighbor states at the previous time step [1]. This approach provides biologically relevant insights as well as computational tractability.

The finite set of possible states implies that, in a deterministic network, the Boolean dynamics will eventually converge to a sequence of states that repeat infinitely, which is denoted as an attractor of the network [2]. An attractor could be a single state (a fixed point) or a cycle consisting of multiple states. It has been shown that attractors have biological interpretations in different contexts, for example as different types of differentiated stem cells [3], states of disease progression [4], or stages of the cell cycle [5, 6].

Reprogramming a stem cell, or driving a cell from a diseased to a healthy state, can be interpreted as applying a control input to steer the network to a desired attractor [7]. Control can be applied to a biological regulatory network by targeting a subset of genes to activate or repress, e.g., through drug therapies. This has been interpreted in the Boolean network framework as pinning a set of genes to a fixed state, corresponding to the desired attractor [8]. The targeted genes then influence the dynamics of their neighbors, eventually steering the entire network state towards the desired attractor.

This approach to control requires selecting a subset of genes that are sufficient to reach the desired attractor from an arbitrary, potentially pathological, initial state. In order to ensure minimal invasiveness and reduce cost, this set of genes should be as small as possible. There are, however, computational challenges associated with selecting a set of targeted genes. First, the number of such sets is exponential in the network size, making exhaustive search impractical. Second, verifying that any given set guarantees convergence to an attractor requires, in the worst case, evaluating convergence from an exponential number of possible initial states. Third, regulatory networks are noisy environments, creating uncertainty in the system model. As a result, existing algorithms for gene selection are either based on approximations from linear system theory (e.g., input selection for controllability) [9], which do not capture the dynamical properties of the regulatory network, or are based on heuristics that inherently cannot provide guarantees on the minimality of the chosen set or the convergence to the desired attractor [10].

In this paper, we propose combinatorial algorithms for selecting a subset of genes to control in order to guarantee convergence to a desired attractor. Our approach is to formulate sufficient conditions for convergence to an attractor that we prove are equivalent to a target set selection (TSS) problem [11]. While TSS is also computationally hard, we identify additional network structures that are common in biological networks and can be exploited to develop efficient approximation algorithms. We make the following specific contributions:

• •

We formulate the problem of selecting a minimum-size subset of nodes in order to guarantee convergence to a desired attractor. We prove a negative result, namely, that there is no approximation guarantee possible for this problem under arbitrary node dynamics unless $P=NP$.

• •

We construct a sufficient condition for convergence to a desired attractor and prove that selecting a minimum-size set that satisfies this condition can be mapped to a TSS problem.

• •

We study the resulting TSS problem under two widely-occurring classes of regulatory networks, namely networks with threshold dynamics and hierarchical structure, and networks with nested canalyzing dynamics. For each type of network, we formulate polynomial-time algorithms that exploit the network structure to provide provable optimality bounds.

• •

We generalize our approach to Boolean networks with probabilistic and asynchronous dynamics. We also show that our approach can be used to select input nodes in order to guarantee convergence to a cyclic attractor.

• •

We evaluate our approach on several real-world biological network datasets as well as randomly generated topologies. We find that our proposed approach requires fewer input nodes to achieve a desired attractor compared to existing heursitics.

The paper is organized as follows. Section II presents related work. Section III presents the system model and definitions, as well as background on the TSS problem. Section IV contains the problem formulation and our proposed gene selection algorithms. Section V generalizes our approach to probabilistic and asynchronous networks. Section VI contains our numerical study. Section VII concludes the paper.

II Related Work

Boolean networks were developed as a computationally tractable approximation to ODE models of biological processes [12, 13, 14, 15]. The concept of attractors was introduced by Waddington [16] and further investigated by Kauffman [2]. The biological relevance of attractors has been confirmed by studies including [3, 4]. Methodologies for inferring regulatory networks from gene expression data have been proposed, e.g., [12]. Generalizations to probabilistic networks were introduced in [shmulevich2002probabilistic]. These works, however, do not consider the problem of selecting a subset of genes to control a regulatory network.

Existing works have modeled therapeutic interventions, such as drug therapies, as inputs to the regulatory network, with the goal of informing possible new treatments [10, 17, 18]. Approaches based on breaking cycles in the network topology were proposed in [17]. Since cycles are very common (indeed, over-represented compared to random networks with similar degree distributions [1]), these approaches may be overly conservative. Heuristics such as genetic algorithms [10] have also been proposed, but do not provide any guarantees on the optimality of the chosen set or on whether convergence to a desired attractor is guaranteed from any initial state.

The problem of controlling a Boolean network is related to selecting input nodes to control a networked dynamical system. Existing works, however, typically consider a linear system with continuous state variables through methods such as controllability analysis [9, 19], which are not applicable to nonlinear biological networks with a limited range of possible control signals.

The target set selection (TSS) problem was first identified in the social networking community [11, 20, 21], and is known to be NP-hard and difficult to approximate [20]. Recent efforts to develop approximation algorithms have focused on special cases of the network topology, such as trees, sparse graphs, and cliques [20, 21]. To the best of our knowledge, application of TSS to biological networks, as well as algorithms that exploit the structural properties of biological structures to reduce the complexity of TSS, have not been studied.

III Model and Background

This section presents the gene regulatory network model, followed by background on the target set selection problem.

III-A Regulatory Network Model

A regulatory network is modeled as a graph $G=(V,E)$ with node set $V$ equal to the set of genes and $E$ denoting the set of edges. The number of vertices $|V|=n$. The network topology is assumed to be directed, with an edge $(i,j)$ implying that node $i$ regulates node $j$. The in-degree of node $i$, denoted $N_{in}(i)=|\{j:(j,i)\in E\}|$ is equal to the number of nodes that regulate $i$, while the out-degree, defined as $N_{out}(i)=|\{j:(i,j)\in E\}|$ is equal to the number of nodes that are regulated by $i$. The graph may contain edges between a node $i$ and itself. For any subset $A\subseteq V$, we let $G(A)=(A,E(A))$, where $E(A)=E\cap A\times A$, denote the subgraph induced $A$.

As a preliminary, for a graph $G=(V,E)$, if $V=V_{1}\cup\cdots\cup V_{m}$ is a disjoint partition of the node set, then we define the graph contraction $\overline{G}$ to be a graph with vertex set $\{1,\ldots,m\}$ and an edge from $i$ to $j$ if the exists an edge $(u,v)$ with $u\in V_{i}$ and $v\in V_{j}$. Note that this construction could result in multiple edges between the same nodes in $\overline{G}$.

Each node $i$ has a Boolean, discrete-time state variable $x_{i}(t)\in\{0,1\}$. Let $\mathbf{x}(t)\in\{0,1\}^{n}$ denote the vector of node states. The state dynamics of node $i$ are given by

[TABLE]

where $f_{i}:\{0,1\}^{n}\rightarrow\{0,1\}$ is a function that determines the state of node $i$ at time $(t+1)$ as a function of its neighbors’ states. The function $f_{i}$ satisfies $f_{i}(\mathbf{x})=f_{i}(\mathbf{x}^{\prime})$ whenever $\mathbf{x}_{j}=\mathbf{x}_{j}^{\prime}$ for all $j\in N_{in}(i)$. The dynamics are written in network form as $\mathbf{x}(t+1)=f(\mathbf{x}(t))$, where $f(\mathbf{x})=(f_{1}(\mathbf{x}),\ldots,f_{n}(\mathbf{x}))$. We define the notations $\vee$, $\wedge$, and $\neg$ denote Boolean OR, AND, and NOT, respectively.

We now discuss two important special cases of the Boolean dynamics $f_{i}$. A node $i$ has threshold dynamics if $f_{i}$ satisfies

[TABLE]

where $\tau_{i}$ is a real-valued threshold and $a_{ij}$ are real-valued coefficients. A positive value of $a_{ij}$ represents an excitatory link from $j$ to $i$, while a negative value represents an inhibitory link.

A node $i$ has nested canalyzing dynamics if $f_{i}$ is defined as follows. Let $d$ denote the in-degree of node $i$, and let $j_{1},\ldots,j_{d}$ be an ordering of $N_{in}(i)$. Let $a_{1},\ldots,a_{d},a_{d+1}\in\{0,1\}$ and $b_{1},\ldots,b_{d}\in\{0,1\}$. Then the nested canalyzing dynamics are defined by

[TABLE]

In words, under nested canalyzing dynamics, there is a ranking of inputs to node $i$. If the top-ranked neighbor $j_{1}$ is in state $b_{1}$, then node $i$ moves to state $a_{1}$. Otherwise, if the neighbor $j_{2}$ is in state $b_{2}$, then $i$ moves to state $a_{2}$, and so on. If none of the conditions are met, then node $i$ reverts to a default state $a_{d+1}$. An attractor of a Boolean network is defined as follows.

Definition 1

An attractor of length $r$ is a sequence of states $\mathbf{x}^{1},\ldots,\mathbf{x}^{r}$ such that $\mathbf{x}^{l}=f(\mathbf{x}^{l-1})$ for $l=2,\ldots,r$ and $\mathbf{x}^{1}=f(\mathbf{x}^{r})$.

An attractor is a set of states that repeat in the Boolean network, so that any network that reaches one of the states in the attractor will remain in the attractor. It can be shown that, for any initial state $\mathbf{x}(0)$, there exists a finite time $T$ such that $\mathbf{x}(T)$ belongs to an attractor.

Attractors can be further classified as fixed points (where $r=1$), limit cycles (where $r$ is small relative to the network size), and chaotic (where $r$ is large).

Lastly, the effect of supplying inputs is described as follows. The set of input nodes is denoted $S$. For any node $i\in S$, there is a variable $s_{i}\in\{0,1\}$ such that $x_{i}(t)\equiv s_{i}$ for all $t$, i.e., each input node is pinned to a fixed state for all time $t$.

III-B The Target Set Selection (TSS) Problem

The target set selection problem is defined as follows. Let $G=(V,E)$ be a graph, and suppose that each node $v\in V$ is assigned a threshold $\tau(v)$. Let $S$ be a subset of nodes.

Initialize set $X[0]=S$. At step $k$, for $k=1,2,\ldots,$ define a set $Y[k]$ by

[TABLE]

and set $X[k]=X[k-1]\cup Y[k]$.Clearly, $X[k]\subseteq X[l]$ when $k<l$.

This process converges when $X[k]=X[k+1]$ for some $k$. Let $X^{\ast}$ be the set $X[k]$ at the iteration when convergence occurs. If $X^{\ast}=V$, then the set $S$ is denoted as a target set of the graph. The target set selection problem is the problem of choosing a minimum-cardinality set for a given graph and set of thresholds.

Proposition 1** ([20])**

The TSS problem is NP-hard.

Although the target set selection problem is NP-hard, tractable algorithms have been found for specific graphs such as complete graphs and graphs with bounded tree-width [20].

IV Input Selection Problem Formulation

This section first formulates the minimum gene selection problem. We analyze the complexity of the problem and then identify a sufficient condition based on target set selection. We provide algorithms for the special cases of threshold and nested canalyzing dynamics.

IV-A Problem Formulation and Complexity

Consider a gene regulatory network defined by a graph $G=(V,E)$ and a set of Boolean functions $\{f_{i}:i\in V\}$. Let $\mathbf{x}^{\ast}\in\{0,1\}^{n}$ denote a fixed-point attractor of the network, i.e., a state satisfying $f(\mathbf{x}^{\ast})=\mathbf{x}^{\ast}$. The case where the attractor consists of multiple states will be considered in Section V.

Definition 2

We say that a set of inputs $S$ guarantees convergence to the desired attractor $\mathbf{x}^{\ast}$ if setting $x_{i}(t)\equiv x_{i}^{\ast}$ for all $t$ implies that $\mathbf{x}(T)=\mathbf{x}^{\ast}$ for $T$ sufficiently large, for any initial state $\mathbf{x}(0)$ with $x_{i}(0)=x_{i}^{\ast}$ when $i\in S$.

We let $\mathcal{C}$ denote the collection of input sets that guarantee convergence to $\mathbf{x}^{\ast}$. The minimum gene selection problem is then formulated as

[TABLE]

We first analyze the complexity of the problem, and find that non-trivial approximations are impossible unless $P=NP$.

Theorem 1

If there exists a function $\gamma:\mathbb{N}\rightarrow\mathbb{N}$ and a polynomial-time algorithm that takes as input an instance of (1) and is guaranteed to output a set $S$ satisfying $S\in\mathcal{C}$ and $|S|\leq\gamma(n)|S^{\ast}|$, where $|S^{\ast}|$ is the optimal solution to (1), then $P=NP$.

Proof:

The proof is by showing that, if there exists an algorithm that satisfies the conditions of the theorem, then that algorithm can also be used to solve 3-SAT, which is an NP-hard problem. An instance of the 3-SAT problem consists of determining whether, given a set of Boolean variables $\{q_{1},\ldots,q_{m}\}$, there exists a set of values for the $q_{i}$’s such that the relation

[TABLE]

where $p_{ij}\in\{q_{s},\neg q_{s}:s\in\{1,\ldots,m\}\}$ for all $i,j$, evaluates to true.

Suppose that an instance of 3-SAT is given. Construct a Boolean network as follows. Let $V=V_{1}\cup V_{2}\cup V_{3}$, where $V_{1}$ is indexed $v_{1}^{1},\ldots,v_{m}^{1}$, $V_{2}$ is indexed $v_{1}^{2},\ldots,v_{l}^{2}$, and $V_{3}$ is a singleton node $v^{3}$.

The edge set is defined as follows. We add an edge $(v_{i}^{1},v_{j}^{2})$ if $p_{js}\in\{q_{i},\neg q_{i}\}$ for some $s$. We add edges $(v_{i}^{1},v_{j}^{1})$ between all nodes in $V_{1}$. We include an edge $(v_{i}^{2},v^{3})$ for $i=1,\ldots,l$. Finally, we add an edge from $v^{3}$ to each other node.

For all nodes $v\in V$, we set $x_{v}(t+1)=1$ if $x_{v^{3}}(t)=1$. Otherwise, the dynamics are defined as follows. We choose the functions $f_{v_{i}^{1}}$ for $i=1,\ldots,m$ so that the binary string $x_{v_{1}^{1}}(t)x_{v_{2}^{1}}(t)\cdots x_{v_{m}^{1}}(t)$ satisfies

[TABLE]

We set

[TABLE]

where $y_{is}(t)=x_{j}(t)$ if $p_{is}=q_{j}$ and $y_{is}(t)=\neg x_{j}(t)$ if $p_{is}=\neg q_{j}$. Furthermore, we set $x_{v^{3}}(t+1)=1$ if $x_{v_{i}^{2}}(t)=1$ for all $i$.

Suppose first that there is a solution to (2), equal to $q_{1}=\overline{q}_{1},\ldots,q_{m}=\overline{q}_{m}$. By construction of the network, for any initial state, eventually there will exist time $T$ such that $x_{v_{i}^{1}}(T)=\overline{q}_{i}$ for all $i=1,\ldots,m$, and hence $x_{v^{3}}(T+2)=1$ and $x_{v}(t)=1$ for all $v\in V$ and $t\geq T+3$. Hence, if there is a solution to (2), then $S^{\ast}=\emptyset$.

On the other hand, suppose there is no solution to (2), and consider an initial state with $x_{v_{i}^{2}}(0)=0$ for all $i$, $x_{v^{3}}(0)=0$, and all other initial states arbitrary. Since there is no solution to (2), $x_{v_{i}^{2}}(t)=x_{v^{3}}(t)=0$ for all $t$, implying that the attractor is never reached when $S=\emptyset$. Convergence to the desired attractor can, however, be achieved by setting $S=\{v^{3}\}$. Hence $S^{\ast}=\{v^{3}\}$ is a minimum-size solution to (1) and $S^{\ast}=\emptyset$ iff the relation (2) is not satisfiable.

If there is a deterministic polynomial-time algorithm that returns a set $S$ satisfying $|S|\leq\gamma(n)|S^{\ast}|$ for some $\gamma$, then that algorithm must choose $S=\emptyset$ whenever $S^{\ast}=\emptyset$. Conversely, $S^{\ast}\neq\emptyset$, then $S\neq\emptyset$. Hence, we can construct a polynomial-time algorithm for 3-SAT by following the above procedure to construct a gene regulatory network, and outputting true if the algorithm returns $\emptyset$ and false otherwise. Thus, if there exists such an algorithm, then $P=NP$. ∎

The proof of Theorem 1 also implies that it is NP-hard to verify whether a given set of input nodes $S$ guarantees convergence to a desired attractor. If not, then it would be possible to check whether there exists a solution to 3-SAT (an NP-complete problem) by verifying whether $S=\emptyset$ guarantees convergence to $\mathbf{x}^{\ast}=\mathbf{1}$ in the graph constructed above.

IV-B Mapping to Target Set Selection

In order to develop efficient approximation algorithms for relevant special cases of (1), we first introduce a sufficient condition that is equivalent to a target set selection problem.

Given a gene regulatory network $G=(V,E)$ and dynamics $f_{i}$ for $i=1,\ldots,n$, we construct an extended network $\hat{G}=(\hat{V},\hat{E})$ as follows. For each Boolean function $f_{i}$, we can write $f_{i}$ in conjunctive normal form as

[TABLE]

where $y_{is}\in\{x_{j},\neg x_{j}\}$ for some $j\in N(i)$. The node set $\hat{V}$ is defined by

[TABLE]

The edge set $\hat{E}$ is defined by

[TABLE]

In words, we have an edge from each $a_{i_{s}}$ to $i$. We also have an edge from $j$ to $a_{i_{s}}$ if $j$ inhibits $i$ and $x_{i}^{\ast}\neq x_{j}^{\ast}$, and an edge from $j$ to $a_{i_{s}}$ if $j$ activates $i$ and $x_{i}^{\ast}=x_{j}^{\ast}$.

The thresholds for this augmented graph are given by

[TABLE]

The threshold that is chosen depends on whether each node is on or off in the desired attractor. We now prove that solving target set selection on this augmented graph is sufficient for ensuring convergence to a desired attractor.

Proposition 2

Suppose that $S\subseteq V$ is a solution to the target set selection problem on the augmented graph $\hat{G}$. Then $S\in\mathcal{C}$.

The proof is omitted due to space constraints. We observe that this condition also gives a polynomial-time algorithm for checking whether a given set of inputs guarantees convergence, namely, allowing the target set dynamics $X[k]$ to unfold for $2n$ iterations on the graph $\hat{G}$. While $V\subseteq X[n]$ implies convergence to the desired attractor from any initial state, the converse is not necessarily true.

We now compare this sufficient condition to a known sufficient condition from previous work [17].

Proposition 3

Suppose that a set $S$ satisfies $S\cap T\neq\emptyset$ for all cycles $T$ in the graph $G$ and each node is path-connected in $G$ to at least one node in $S$. Then the set $S$ is a target set for the augmented graph $\hat{G}$.

Proof:

Suppose that $S$ satisfies the conditions of the theorem, and yet $S$ is not a target set. Then there exists $i\in V$ such that $i\notin X[k]$ for any $k$. We must have $i\notin S$. If $i$ has no neighbors in $\hat{V}$, then $i$ is not connected to any input, a contradiction. Otherwise, by construction, there exists at least one $a_{i_{s}}$ such that $a_{i_{s}}\notin X[k]$, and therefore at least one neighbor $j\in N(i)$ such that $x_{j}\notin X^{\ast}$. If $j=i$, then there is a cycle $T$, consisting of the self-loop $(i,i)$, such that $S\cap T=\emptyset$, a contradiction.

Proceeding inductively, we maintain a set $U\subseteq V\setminus X^{\ast}$ where each node in $U$ is path-connected to $i$. Within $2n$ iterations, we must have a node that is added to $U$ twice, implying the existence of a cycle in the graph that is disjoint from $S$ and yielding a contradiction. ∎

Proposition 3 implies that the target set selection condition is weaker (easier to satisfy) than the current known approach of selecting a subset of nodes to break all cycles in the graph. On the other hand, this approach also requires adding new nodes and edges to the underlying the graph, and moreover, relies on solving the computationally difficult TSS problem. In the general case, this problem can be formulated as the integer linear program [11]

[TABLE]

The binary variables $\{s_{i}:i=1,\ldots,n\}$ satisfy $s_{i}=1$ iff $i\in S$. Hence, the solution to TSS can be obtained by solving (3) and selecting the set $S$ based on the $s_{i}$’s.

In the next subsections, we examine biologically relevant special cases and develop algorithms that exploit these additional structures.

IV-C Threshold Dynamics

We now analyze the target set selection formulation in the special case where the nodes have threshold dynamics. Specifically, we assume threshold dynamics where, for each node $i$, either $a_{ij}=1$ for all $j\in N_{out}(i)$ or $a_{ij}=-1$ for all $j\in N_{out}(i)$. Intuitively, all nodes are either purely excitatory ($a_{ij}=1)$ or purely inhibitory ($a_{ij}=-1$), and all nodes exert the same effect on their neighbors. Let $\mathcal{E}\triangleq\{i:a_{ij}=1\ \forall j\in N_{out}(i)\}$ and $\mathcal{I}\triangleq\{i:a_{ij}=-1\ \forall j\in N_{out}(i)\}$. For the given attractor $\mathbf{x}^{\ast}$, define $\mathcal{E}^{1}\triangleq\mathcal{E}\cap\{i:x_{i}^{\ast}=1\}$ and $\mathcal{E}^{0}\triangleq\mathcal{E}\cap\{i:x_{i}^{\ast}=0\}$, and define $\mathcal{I}^{1}$ and $\mathcal{I}^{0}$ in an analogous manner.

In the case of threshold dynamics, the target set selection instance is as follows. The graph $\hat{G}$ has node set $\hat{V}=V$. An edge $(i,j)$ exists from node $i$ to node $j$ if $x_{i}^{\ast}=x_{j}^{\ast}$ and node $i\in\mathcal{E}$, or if $x_{i}^{\ast}\neq x_{j}^{\ast}$ and $i\in\mathcal{I}$. The threshold $\hat{\tau}(i)$ is given by

[TABLE]

Lemma 1

Suppose that $S$ is a target set for the graph $\hat{G}$ with thresholds $\hat{\tau}$. Then there exists $T$ such that $\mathbf{x}(t)=\mathbf{x}^{\ast}$ for all $t\geq T$.

Proof:

We show that if $i\in X[k]$ for some $k>0$, then $x_{i}(t)$ converges to $x_{i}^{\ast}$. The proof is by induction on $k$, noting that $X[0]=S$. At time $k$, suppose that $i\in X[k]\setminus X[k-1]$, and hence the threshold condition is satisfied. We have that

[TABLE]

which is equivalent to

[TABLE]

Suppose that $x_{i}^{\ast}=1$; the case where $x_{i}^{\ast}=0$ is similar. Then (4) is equivalent to

[TABLE]

By inductive hypothesis, for $t$ sufficiently large, $x_{j}(t)=x_{j}^{\ast}$ for all $j\in X[k-1]$. Thus for $t$ sufficiently large, we have

[TABLE]

Hence $x_{i}(t)=1$ for $t$ sufficiently large, implying that the desired attractor is reached. The fact that $X^{\ast}=V$ completes the proof. ∎

Under threshold dynamics, formulating the input selection problem as TSS does not require adding any nodes to the graph $G$. We analyze algorithms for computing optimal target sets under network topologies that typically arise in biological networks. We first consider complete graphs, which arise as subgraphs of regulatory networks, and consider a generalization of known results to networks with both positive and negative edges.

If the graph $G$ is complete, then the edge set will be given by

[TABLE]

Lemma 2

Let $S$ be a minimum-size target set for a complete graph. Suppose that there exist two nodes $i$, $j$ such that $i$ and $j$ both lie in $\mathcal{E}^{1}$, $\mathcal{E}^{0}$, $\mathcal{I}^{1}$, or $\mathcal{I}^{0}$, $i\in S$, $j\notin S$, and $\hat{\tau}(i)<\hat{\tau}(j)$. Then $S\setminus\{i\}\cup\{j\}$ is also a solution to the target set selection problem.

Proof:

Let $X[k]$ denote the thresholding process on $\hat{G}$ when the initial set $X[0]=S$, and let $\overline{X}[k]$ denote the thresholding process when the initial set $\overline{X}[0]=\overline{S}=S\setminus\{i\}\cup\{j\}$. We will show that the result holds when $\{i,j\}\subseteq\mathcal{E}^{1}$; other cases are similar and omitted due to space constraints.

It suffices to show that $i\in\overline{X}[k]$ for some $k$. This will hold when

[TABLE]

Define

[TABLE]

Now, since it is assumed that $S$ is a target set and $j\notin S$, we must have $\alpha_{k}\geq\overline{\tau}(j)$ for some $k$. Hence, if we can show that $\overline{\alpha}_{k}\geq\alpha_{k}$ and $\overline{\beta}_{k}\geq\beta_{k}$ for all $k$, then (5) will be satisfied, since $\overline{\alpha}_{k}\geq\alpha_{k}\geq\overline{\tau}(j)\geq\overline{\tau}(i)$ for some $k$.

The proof that $\overline{\alpha}_{k}\geq\alpha_{k}$ and $\overline{\beta}_{k}\geq\beta_{k}$ is by induction on $k$. When $k=1$, $X[0]=S$ and $\overline{X}[0]=\overline{S}$, and hence $\alpha_{k}=\overline{\alpha}_{k}$ and $\overline{\beta}_{k}=\beta_{k}$ by definition of $S$ and $\overline{S}$. For larger values of $k$, we have that

[TABLE]

with similar identities for $\mathcal{E}^{0}$, $\mathcal{I}^{1}$, and $\mathcal{I}^{0}$. If $\overline{\alpha}_{k-1}\geq\alpha_{k-1}$ and $\overline{\beta}_{k-1}\geq\beta_{k-1}$, then $(\mathcal{I}^{0}\cap X[k])\subseteq(\mathcal{I}^{0}\cap\overline{X}[k]$, with similar identities for $\mathcal{E}^{0}$ and $\mathcal{I}^{1}$.

To show that $(\mathcal{E}^{1}\cap X[k])\subseteq(\mathcal{E}^{1}\cap\overline{X}[k])$, we first have that

[TABLE]

Hence $((\mathcal{E}^{1}\cap X[k-1])\setminus(\mathcal{E}^{1}\cap\overline{X}[k-1]))\subseteq\{i\}.$ If the above holds with equality, then $\overline{\alpha}_{k-1}<\overline{\tau}(i)$, and hence by inductive hypothesis $\alpha_{k-1}<\overline{\tau}(j)$, implying that $j\notin(\mathcal{E}^{1}\cap X[k-1]$. We therefore have that

[TABLE]

and hence $(\mathcal{E}^{1}\cap X[k])\subseteq(\mathcal{E}^{1}\cap\overline{X}[k])$. Thus $\overline{\alpha}_{k}\geq\alpha_{k}$ and $\overline{\beta}_{k}\geq\beta_{k}$ for all $k$, completing the proof. ∎

Lemma 2 implies that, when the graph is a clique, we can restrict the search space by ordering the vertices in $\mathcal{E}^{1}$, $\mathcal{E}^{0}$, $\mathcal{I}^{0}$ and $\mathcal{I}^{1}$ based on their thresholds, and choosing the $m(\mathcal{E}^{1})$ (resp. $m(\mathcal{E}^{0})$, $m(\mathcal{I}^{1})$, $m(\mathcal{I}^{0})$) vertices with largest threshold from $\mathcal{E}^{1}$ (resp. $\mathcal{E}^{0}$, $\mathcal{I}^{1}$, $\mathcal{I}^{0}$). For a network of $n$ nodes, there are no more than $n^{4}$ sets of this type, implying that the selection problem on a clique can be solved in $O(n^{4})$ time.

Next, we consider graphs that have a block cactus structure, in which the set of vertices $V$ can be partitioned as $V=V_{1}\cup\cdots\cup V_{m}$, where the subgraph induced by each $V_{i}$ is a clique, and the contraction of the graph around the $V_{i}$’s is a tree.

Proposition 4

There exists an $O(n^{4})$ algorithm for computing a minimum-size set of input genes in a block cactus graph.

Proof:

The proof is by induction on $m$. When $m=1$, the problem reduces to selecting a minimum-size input set for a clique. Suppose that the result holds up to $(m-1)$, and consider a tree of $m$ cliques. Without loss of generality, suppose that $V_{m}$ corresponds to a leaf in the tree, i.e., there is exactly one incoming edge that is not part of the clique. This assumption is without loss of generality because at least one node in the tree must be a leaf (a degree-one vertex), and hence we can always have $V_{m}$ as a leaf by reordering vertices.

We consider two cases on $V_{m}$. First, suppose that the only edge incident on $V_{m}$ is an incoming edge onto a node denoted $v\in V_{m}$. Define a new threshold vector for $V_{m}$ by $\hat{\tau}(v)=\overline{\tau}(v)-1$ and $\hat{\tau}(u)=\tau(u)$ for $u\neq v$. Let $S_{1}$ be a minimum-size input set for $V_{1}\cup\cdots\cup V_{m-1}$, and let $S_{2}$ be a minimum-size input set for $V_{m}$ when the threshold is equal to $\hat{\tau}(v)$. We then have that $S_{1}\cup S_{2}$ is a minimum-size input set for $G$, and the computation time is equal to $O((n-|V_{m}|)^{4})+O(|V_{m}|^{4})=O(n^{4})$.

Conversely, suppose that the only edge incident on $V_{m}$ is an outgoing edge. Define a new threshold vector for $V_{1}\cup\cdots\cup V_{m-1}$ as $\overline{\tau}(v)=\tau{\tau}(v)-1$, where $v$ is the node that has an incoming edge from $V_{m}$. Then the selection algorithm is equivalent to choosing a set of nodes to ensure that $V_{m}$ reaches the desired attractor and a set of nodes to ensure that $V_{1}\cup\cdots\cup V_{m-1}$ reaches the desired attractor with thresholds $\hat{\tau}$, requiring only $O(n^{4})$ time in total. ∎

The algorithm for selecting a minimum-size set of input nodes is described as Algorithm 1.

Networks that do not have block cactus structure can be addressed within this framework by grouping the network nodes into densely-connected clusters, denoted $V_{1},\ldots,V_{m}$ (e.g., via the methods in [22]). For any two nodes $i$ and $j$ in the same cluster that are not connected by an edge, add an edge and increment $\hat{\tau}(i)$ and $\hat{\tau}(j)$ by $1$. Then, find a set of edges $E^{\prime}$ to remove in order to remove cycles from the contracted graph; such a set of edges corresponds to a minimum arc feedback set. The thresholds are unchanged at this stage.

In addition to consisting of loosely connected components, biological networks also often have modular, hierarchical structure. These modular structures are believed to be derived from the functional organization of cells. In the following, we analyze gene selection algorithms on a model of hierarchical networks introduced in [23]. We first define the model as follows.

The hierarchical network is constructed iteratively. The network originates with a single hub node. At the first iteration, $k$ nodes are added and are connected to each other, creating a network $G_{1}$. At the $i$-th iteration, $k$ copies of $G_{i-1}$ are generated and connected to the hub node. An algorithm for constructing a minimum-size set of genes to control a hierarchical network. The graph is undirected, implying that for each edge $(i,j)$, there is a corresponding edge $(j,i)$.

The algorithm is recursive. For a hierarchical network with $d$ levels, each consisting of $m$ copies, the approach is to compute, for each copy, the number of input nodes required with and without the central hub node in $S$. For all the sub-graphs where the number of input nodes is the same under both cases, select a subset of input nodes to drive the subgraph to the desired state. If the nodes in these sub-graphs are insufficient to drive the hub node over its threshold, then add the central hub node to the input set, recompute all remaining thresholds, and compute sets of inputs to guarantee that the remaining subgraphs reach the desired state.

Proposition 5

Algorithm 2 selects a set $S$ satisfying $|S|\leq|S^{\ast}|\log{n}$, where $S^{\ast}$ is the minimum-size input set, in $O(n^{2})$ time.

Proof:

To show the complexity, let $R(k,m)$ denote the number of computations required to compute the optimal set in a graph with $k$ copies and $m$ iterations. We have that $R(k,m+1)=2kR(k,m)$, and hence $R(k,m)=(2k)^{m}$. At the same time, the number of nodes is equal to $n=(k+1)^{m}$. Hence we have

[TABLE]

and hence $R(k,m)=n^{2}$.

We then analyze the optimality gap. Let $\epsilon(m)$ denote the worst-case optimality gap in a network with $m$ levels of hierarchy. We then have that

[TABLE]

On the other hand, $|S^{\ast}|\geq\sum_{i=1}^{m}{|S_{i}^{\ast}|}$, and hence combining these expressions yields

[TABLE]

implying that $\epsilon(m+1)\leq\epsilon(m)+1$ and hence $\epsilon(m)\leq m$. Since $m=\frac{\log{n}}{\log{(k+1)}}\leq\log{n}$, we have the desired optimality bound. ∎

IV-D Nested Canalyzing Dynamics

We now present sufficient conditions for selecting genes to ensure convergence in networks with nested canalyzing dynamics, as defined in Section III-A. We first characterize the sufficient condition of Proposition 2 for this class of dynamics.

Lemma 3

For each node $i$, define $\Omega_{i}=\{r:a_{r}=x_{i}^{\ast}\}$, with $r_{i}=|\Omega_{i}|$. Then the following instance of TSS is sufficient to ensure convergence to the desired attractor. For each node $i$, define a collection of nodes $u_{i,1},\ldots,u_{i,r_{i}}$. Each node $u_{i,a_{s}}$ has an incoming edge from each node $j_{l}\in N(i)$ with $l<s$ and $a_{l}\neq x_{i}^{\ast}$, an incoming edge from $j_{s}$, and a threshold equal to the degree of $u_{i,j_{s}}$. Each node $i$ has threshold $1$ in the graph.

Proof:

The nested canalyzing dynamics are equivalent to the condition

[TABLE]

Applying the construction of Section IV-B yields the graph described in the statement of the lemma. ∎

The following corollary provides a condition that admits computationally tractable approximation algorithms.

Corollary 1

Let $s_{i}^{\ast}=\min{\{s:a_{i,s}=x_{i}^{\ast}\}}$. Consider an instance of the TSS problem defined by a graph $\hat{G}=(V,\hat{E})$, in which there is an edge $(j_{s},i)\in\hat{E}$ if $s\leq s_{i}^{\ast}$, where each node’s threshold is equal to the degree of the node. Then a solution to this instance of TSS is sufficient to ensure convergence to a desired attractor.

The instance of TSS defined by Corollary 1 has a desirable structure, namely each node has a threshold equal to its degree (a unanimous threshold), equivalent to a Boolean AND decision rule. In undirected graphs, it is known that this Boolean decision rule is equivalent to the vertex cover problem [24]. The following gives a necessary and sufficient condition for directed graphs

Proposition 6

The condition of Corollary 1 holds if and only if each cycle in $\hat{G}$ contains at least one node from $S$ and each node is connected to at least one node in $S$.

Proof:

First, suppose that a set $S$ does not satisfy the conditions of Corollary 1, and yet the two conditions of the proposition hold. Let $i$ be a node satisfying $X_{i}^{\ast}=0$. Then either $i$ is an isolated node, contradicting the assumption that all isolated nodes are in $S$, or there exists a neighbor, denoted $i_{1}$, satisfying $X_{i_{1}}^{\ast}=0$. Proceeding inductively, we obtain a set of nodes $i_{0},i_{1},\ldots,i_{r}$ that all satisfy $X_{i_{j}}^{\ast}=0$, and must either contain a cycle or have $X_{i_{r}}^{\ast}$ isolated, contradicting the conditions of the proposition.

Clearly if there is an isolated node $i\notin S$, then $X_{i}[k]\equiv 0$ for all $k$. Similarly, suppose that there is a set of edges $(i_{0},i_{1}),\ldots,(i_{m},i_{0})$ such that $S\cap\{i_{0},\ldots,i_{m}\}=\emptyset$. Then $X_{i_{l}}[0]=0$, and by induction $X_{i_{l}}[k]=0$ for all $k$ since there exists a neighbor $i_{l-1}$ satisfying $X_{i_{l-1}}[k-1]=0$. ∎

Note that this condition is the same as that of [17], but for the subgraph $\hat{G}$. Based on Proposition 6, we introduce an algorithm for selecting input genes under nested canalyzing dynamics as Algorithm 3.

Intuitively, Algorithm 3 is as follows. We first compute the maximal strongly connected subgraphs of $G$, which can be done in polynomial time, and contract with respect to these components to obtain a directed acyclic graph $\overline{G}$. It then suffices to ensure that each subgraph $G_{i}$ has no cycles that are disjoint from $S$, as well as ensuring that all nodes are connected to a node in $S$. This condition is ensured if each component is cycle-free and if all singleton isolated components (which are exactly the isolated nodes of $\hat{G}$) are in $S$. The optimality guarantees of this approach are given in Proposition 7.

Proposition 7

Let $S^{\ast}$ denote the optimal set for the sufficient condition of Corollary 1, and let $S$ denote the set returned by Algorithm 3. Then $|S|\leq(\log{n})^{2}|S^{\ast}|$.

Proof:

Define $S_{i}^{\ast}=S^{\ast}\cap[i]$, and $S_{i}=S\cap[i]$, so that $S^{\ast}=S_{1}^{\ast}\cup\cdots\cup S_{m}^{\ast}$ and $S=S_{1}\cup\cdots\cup S_{m}$. By [25], $|S_{i}|\leq(\log{n})^{2}|S_{i}^{\ast}|$, and hence

[TABLE]

as desired. ∎

V Generalizations to Probabilistic Graphs and Cyclic Attractors

In this section, we investigate two generalizations to the problem formulation. We first investigate probabilistic Boolean networks, followed by guaranteeing convergence to attractors with multiple states.

V-A Probabilistic Regulatory Networks

Probabilistic Boolean networks are an extension of Boolean regulatory networks to model the inherent uncertainty of biological systems. A Boolean regulatory network is defined by a graph $G=(V,E)$ and a set of $K$ update functions $f(\cdot,1),\ldots,f(\cdot,K)$ each of which maps $2^{|V|}$ into $2^{|V|}$. The Boolean network is also characterized by a random process $\xi(t)\in\{1,\ldots,K\}$, so that $\mathbf{x}(k+1)=f(\mathbf{x}(k),\xi(t))$.

A generalization of the approach of Section IV-B is as follows. Let

[TABLE]

be a CNF realization of the dynamics of node $i$ in topology $j$. Furthermore, define an extended function $f$ by

[TABLE]

where $f_{i}(\mathbf{x}_{-i},x_{i}^{\ast},j)$ refers to the value of $f_{i}(\cdot,j)$ when $x_{i}=x_{i}^{\ast}$ and all other indices are equal to $\mathbf{x}$. Define $\overline{f}_{i,0}$ by

[TABLE]

The definition of $f_{i,1}$ (resp. $f_{i,0}$) is chosen so that, if $\overline{f}_{i,1}(\mathbf{x}^{\ast})=1$, then $f_{i}(\mathbf{x},j)=1$ for at least one function $j$. Furthermore, if node $i$ achieves the desired attractor, then the function will remain in the desired state.

Proposition 8

Let $G=(V,E)$ be a Boolean network with update functions $f(\cdot,1),\ldots,f(\cdot,K)$. Construct an instance of TSS based on the approach of Section IV-B, using the functions $\overline{f}_{i,0}$ and $\overline{f}_{i,1}$. Then the resulting set $S$ is sufficient to guarantee convergence to an attractor $\mathbf{x}^{\ast}$, provided that $f(\mathbf{x}^{\ast},j)=\mathbf{x}^{\ast}$ for all $j$.

Proof:

The approach is to show by induction that, if $i\in X[k]$, then $x_{i}(t)$ eventually reaches $x_{i}^{\ast}$ regardless of the set of topologies. Suppose the result is true up to iteration $k$. By construction of $\overline{f}_{i,0}$, there exists at least one time step $T$ such that $x_{i}(T)=x_{i}^{\ast}$, and by construction of $\overline{f}_{i,1}$, $x_{i}(t)=x_{i}^{\ast}$ for all $t\geq T$. ∎

In the special case of threshold dynamics, we have the following corollary.

Corollary 2

For networks with threshold dynamics, if a node is excitatory or inhibitory under all functions $f^{(i)}$ and only the threshold varies, then the threshold $\overline{\tau}_{i}=\max{\{\tau_{i}^{(1)},\ldots,\tau_{i}^{(K)}\}}$ is sufficient to ensure convergence.

Another relevant class of models arises when different nodes may update their states asynchronously. The resulting network has $K=n$, where $n$ is the number of nodes. The functions $f(\cdot,i)$ are defined by

[TABLE]

so that only node $i$ updates its state and all other nodes maintain fixed state values.

Lemma 4

Suppose that the set $S$ satisfies the conditions of Proposition 2. Then the set $S$ is sufficient to guarantee convergence to a desired attractor under asynchronous dynamics.

Proof:

Let $i_{1},i_{2},\ldots,i_{n}$ denote the sequence in which nodes are activated by the process $X[k]$. Define $T_{1}$ as $T_{1}=\min{\{t:\xi(t)=i_{1}\}}$ and $T_{j}$ for $j\geq 2$ as

[TABLE]

We have that $x_{i_{1}}(T_{1})=x_{i_{1}}^{\ast}$. Proceeding inductively, at time $T_{j}$, a set of neighbors of $N_{i_{j}}$ has reached the desired attractor, which is sufficient to ensure that $X_{i_{j}}[k]=1$. Hence, $x_{i_{j}}(t)=x_{i_{j}}^{\ast}$ for $t\geq T_{j}$. ∎

V-B Convergence to Cyclic Attractors

We now remark on the gene selection problem to ensure convergence to a cyclic attractor, i.e., an attractor that consists of multiple states. We formulate a TSS-based condition that is analogous to the fixed point condition of Section IV-B. Let $(\mathbf{x}^{1},\ldots,\mathbf{x}^{p})$ denote the desired attractor.

The approach is to construct a network graph $\hat{G}=(\hat{V},\hat{E})$, where $\hat{V}=\hat{V}_{1}\cup\cdots\cup\hat{V}_{p}$. Using the CNF form, we have $\hat{V}_{a}=\{i_{a}:i\in V\}\cup\{j_{i_{s},a}:s=1,\ldots,l\}$ for $a=1,\ldots,p$. The edge set is defined by

[TABLE]

This definition is analogous to Section IV-B, except there is an edge from the $a$-th copy of $V$ to the $(a+1)$-th copy if the state of node $j$ in the $a$-th state of the attractor influences the state of node $i$ in the $(a+1)$-th state of the attractor. The thresholds are defined as in Section IV-B.

Proposition 9

If $S$ is a target set for the graph $\hat{G}$ with thresholds $\tau$, then controlling the set of genes $S$ ensures convergence to the desired cyclic attractor.

The proof is omitted due to space constraints.

VI Numerical Study

We conducted a numerical study of our approach using MatlabTM. The goals of our numerical study were two-fold. The first objective was to evaluate the behavior of our approach on a real-world biological network. The second objective was to observe trends in the number of input nodes required to converge to a desired attractor, as a function of parameters such as the class of network (e.g., scale-free or Erdos-Renyi graph), the average node degree, and the number of nodes in the network. For all threshold networks, a fixed-point attractor was computed by solving an integer linear program.

In order to complete the first objective, we obtained several biological regulatory networks from the Cell Collective website [26], which maintains an archive of biological networks including the topology and the Boolean dynamics of each node. We evaluated our approach on several networks, including the T-cell differentiation network, the yeast cell cycle network, an Apoptosis network, a model of the mammalian cell cycle, the T-cell signaling network, and the regulatory network of Bordatella Bronchiseptica. For the node dynamics, we used threshold dynamics with threshold [math], in which nodes were chosen as excitatory or inhibitory based on the published Boolean dynamics. The results are summarized in Table I.

From Table I, we observe that the number of inputs required is typically a small fraction of the number of network nodes, implying that only a few inputs are needed to guarantee convergence to a desired attractor. The exception to this rule is the T-Cell differentiation network, in which nearly half of the nodes must act as inputs.

We then evaluated the behavior of our algorithms on synthetic networks with different topologies. We first performed a comparison of classical random graph models with graph models that more closely approximate regulatory networks. We chose Erdos-Renyi and scale-free graphs for comparison. In an Erdos-Renyi graph, each node is connected to each other node with a fixed probability $p$. In a scale-free graph, each node is connected to $m$ randomly chosen nodes, where the probability of an edge is proportional to the degree of the node (preferential attachment model). As shown in Figure 1(a), we found that scale-free networks consistently require fewer input nodes than Erdos-Renyi random graphs to guarantee convergence. This may be due to the presence of high-degree hubs and increased clustering in scale-free networks. Indeed, the number of input nodes required by the scale-free network did not increase as a function of the network size.

We studied the effect of the degree distribution on the number of input nodes required. We considered scale-free networks in which the degree distribution is varied by changing $m$. For all cases we considered threshold dynamics with threshold [math] and edges randomly assigned as excitatory or inhibitory with probability $0.5$. We found that high-degree networks required fewer inputs, as a subset of well-connected hub nodes are sufficient to guarantee convergence (Figure 1(b)).

Finally, we compared our approach to a current state of the art approach (Figure 1(c)), which is based on selecting a minimum-size set of inputs such that all cycles contain at least one input [17]. The network considered was a scale-free graph with $m=0.2n$, where $n$ is the number of nodes. The cycle-based method consistently required more input nodes than our proposed TSS-based algorithm. This result agrees with the theoretical guarantees of Proposition 3.

VII Conclusions and Future Work

This paper investigated the problem of selecting input nodes to control biological regulatory networks. Under a Boolean network model, we formulated the problem of selecting a minimum-size set of inputs to guarantee convergence to a desired attractor, defined as a stable fixed point of the network dynamics, and showed that this problem cannot be approximated up to any provable bound unless P=NP. We showed that a sufficient condition for convergence can be mapped to an instance of the target set selection problem, which is defined as selecting a minimum-size set of nodes to ensure that all nodes are activated by a threshold dynamics.

We analyzed our sufficient condition under biologically relevant special cases of the network dynamics. For threshold dynamics with modular structure, we proposed polynomial-time exact algorithms for input selection. In networks with hierarchical structure, we introduced an $O(n^{2})$ algorithm that selects a minimum-size input set up to a provable bound of $\log{n}$. Finally, in networks with nested canalyzing dynamics, we showed that a sufficient condition for convergence to a desired attractor is ensuring that each cycle in a subgraph contains at least one input node, leading to polynomial-time algorithms with an optimality bound of 2. We proposed generalizations of our approach to asynchronous and probabilistic dynamics, as well as multi-state attractors.

We plan to investigate tighter sufficient conditions and exploit additional network structures to improve computation times. Furthermore, other control actions, such as time-varying interventions and changes in the network topology, will be considered in future. Finally, computing the number of distinct minimum-size input sets is an additional open research problem.