Control of Asymmetric Hopfield Networks and Application to Cancer Attractors

The asymmetric Hopfield model is used to simulate signaling dynamics in gene regulatory networks. The model allows for a direct mapping of a gene expression pattern into attractor states. We analyze different control strategies aimed at disrupting attractor patterns using selective local fields representing therapeutic interventions. The control strategies are based on the identification of signaling bottlenecks, which are single nodes or strongly connected clusters of nodes that have a large impact on the signaling. We provide a theorem with bounds on the minimum number of nodes that guarantee control of bottlenecks consisting of strongly connected components. The control strategies are applied to the identification of sets of proteins that, when inhibited, selectively disrupt the signaling of cancer cells while preserving the signaling of normal cells. We use an experimentally validated non-specific and an algorithmically-assembled specific B cell gene regulatory network reconstructed from gene expression data to model cancer signaling in lung and B cells, respectively. Among the potential targets identified here are TP53, FOXM1, BCL6 and SRC. This model could help in the rational design of novel robust therapeutic interventions based on our increasing knowledge of complex gene signaling networks.


Introduction
The vision behind systems biology is that complex interactions and emergent properties determine the behavior of biological systems. Many theoretical tools developed in the framework of spin glass models are well suited to describe emergent properties, and their application to large biological networks represents an approach that goes beyond pinpointing the behavior of a few genes or metabolites in a pathway. The Hopfield model [1] is a spin glass model that was introduced to describe neural networks, and that is solvable using mean field theory [2]. The asymmetric case, in which the interaction between the spins can be seen as directed, can also be exacty solved in some limits [3]. The model belongs to the class of attractor neural networks, in which the spins evolve towards stored attractor patterns, and it has been used to model biological processes of high current interest, such as the reprogramming of pluripotent stem cells [4]. Moreover, it has been suggested that a biological system in a chronic or therapyresistant disease state can be seen as a network that has become trapped in a pathological Hopfield attractor [5]. A similar class of models is represented by Random Boolean Networks [6], which were proposed by Kauffman to describe gene regulation and expression states in cells [7]. Differences and similarities between the Kauffman-type and Hopfield-type random networks have been studied for many years [8][9][10][11].
In this paper, we consider an asymmetric Hopfield model built from real (even if incomplete [12,13]) cellular networks, and we map the spin attractor states to gene expression data from normal and cancer cells. We will focus on the question of controling of a network's final state (after a transient period) using external local fields representing therapeutic interventions. To a major extent, the final determinant of cellular phenotype is the expression and activity pattern of all proteins within the cell, which is related to levels of mRNA transcripts. Microarrays measure genome-wide levels of mRNA expression that therefore can be considered a rough snapshot of the state of the cell. This state is relatively stable, reproducible, unique to cell types, and can differentiate cancer cells from normal cells, as well as differentiate between different types of cancer [14,15]. In fact, there is evidence that attractors exist in gene expression states, and that these attractors can be reached by different trajectories rather than only by a single transcriptional program [16]. While the dynamical attractors paradigm has been originally proposed in the context of cellular developement, the similarity between cellular ontogenesis, i.e. the developement of different cell types, and oncogenesis, i.e. the process under which normal cells are transformed into cancer cells, has been recently emphasized [17]. The main hypothesis of this paper is that cancer robustness is rooted in the dynamical robustness of signaling in an underlying cellular network. If the cancerous state of rapid, uncontrolled growth is an attractor state of the system [18], a goal of modeling therapeutic control could be to design complex therapeutic interventions based on drug combinations [19] that push the cell out of the cancer attractor basin [20].
Many authors have discussed the control of biological signaling networks using complex external perturbations. Calzolari and coworkers considered the effect of complex external signals on apoptosis signaling [21]. Agoston and coworkers [22] suggested that perturbing a complex biological network with partial inhibition of many targets could be more effective than the complete inhibition of a single target, and explicitly discussed the implications for multi-drug therapies [23]. In the traditional approach to control theory [24], the control of a dynamical system consists in finding the specific input temporal sequence required to drive the system to a desired output. This approach has been discussed in the context of Kauffmann Boolean networks [25] and their attractor states [26]. Several studies have focused on the intrinsic global properties of control and hierarchical organization in biological networks [27,28]. A recent study has focused on the minimum number of nodes that needs to be addressed to achieve the complete control of a network [29]. This study used a linear control framework, a matching algorithm [30] to find the minimum number of controllers, and a replica method to provide an analytic formulation consistent with the numerical study. Finally, Cornelius et al. [31] discussed how nonlinearity in network signaling allows reprogrammig a system to a desired attractor state even in the presence of contraints in the nodes that can be accessed by external control. This novel concept was explicitly applied to a T-cell survival signaling network to identify potential drug targets in T-LGL leukemia. The approach in the present paper is based on nonlinear signaling rules and takes advantage of some useful properties of the Hopfield formulation. In particular, by considering two attractor states we will show that the network separates into two types of domains which do not interact with each other. Moreover, the Hopfield framework allows for a direct mapping of a gene expression pattern into an attractor state of the signaling dynamics, facilitating the integration of genomic data in the modeling.
The paper is structured as follows. In Mathematical Model we summarize the model and review some of its key properties. Control Strategies describes general strategies aiming at selectively disrupting the signaling only in cells that are near a cancer attractor state. The strategies we have investigated use the concept of bottlenecks, which identify single nodes or strongly connected clusters of nodes that have a large impact on the signaling. In this section we also provide a theorem with bounds on the minimum number of nodes that guarantee control of a bottleneck consisting of a strongly connected component. This theorem is useful for practical applications since it helps to establish whether an exhaustive search for such minimal set of nodes is practical. In Cancer Signaling we apply the methods from Control Strategies to lung and B cell cancers. We use two different networks for this analysis. The first is an experimentally validated and non-specific network (that is, the observed interactions are compiled from many experiments conducted on heterogeneous cell types) obtained from a kinase interactome and phospho-protein database [32] combined with a database of interactions between transcription factors and their target genes [33]. The second network is cell-specific and was obtained using network reconstruction algorithms and transcriptional and post-translational data from mature human B cells [34]. The algorithmically reconstructed network is significantly more dense than the experimental one, and the same control strategies produce different results in the two cases. Finally, we close with Conclusions.

Mathematical Model
We define the adjacency matrix of a network G composed of N nodes as where j?i denotes a directed edge from node j to node i. The set of nodes in the network G is indicated by V (G) and the set of directed edges is indicated by Table 1 for a list of mathematical symbols used in the text.) The spin of node i at time t is s i (t)~+1, and indicates an expresssed (z1) or not expressed ({1) gene. We encode an arbitrary attractor statẽ j j~(j 1 ,j 2 ,:::,j N ) with j i~+ 1 by defining the coupling matrix [1] The total field at node i is then h i~h ext i z P j J ij s j , where h ext i is the external field applied to node i, which will be discussed below. The discrete-time update scheme is defined as where T §0 is an effective temperature. For the remainder of the paper, we consider the case of T~0 so that s i~s ign(h i ), and the spin is chosen randomly from +1 if h i~0 . For convenience, we take t[ and Dt~1. Nodes can be updated synchronously, and synchronous updating can lead to limit cycles [9]. Nodes can also be updated separately and in random order (anynchronous updating), which does not result in limit cycles. All results presented in this paper use the synchronous update scheme.
Source nodes (nodes with zero indegree) are fixed to their initial states by a small external field so that s q (t)~s q (0) for all q[Q, where Q is the set of source nodes. However, the source nodes flip if directly targeted by an external field. Biologically, genes at the ''top'' of a network are assumed to be controlled by elements outside of the network.
In application, two attractors are needed. Define these states as j j n andj j c , the normal state and cancer state, respectively. The magnetization along attractor state a is Note that if m a (t)~+1,s s(t)~+j j a . We also define the steady state magnetization along state a as m a ?~l im There are two ways to model normal and cancer cells. One way is to simply define a different coupling matrix for each attractor state a, Alternatively, both attractor states can be encoded in the same coupling matrix, Systems using Eqs. 6 and 7 will be referred to as the one attractor state (p~1) and two attractor state (p~2) systems, respectively. Eqs. 6 and 7 are particular cases of the general Hopfield form [1] where p is the number of attractor states, often taken to be large. An interesting property emerges when p~2, however. Consider a simple network composed of two nodes, with only one edge 1?2 with attractor statesj j n andj j c , and T~0. The only nonzero entry of the matrix J ij is Note that ifj j n~+j j c , J 21~2 j n 2 j n 1 . In either case, by Eq. 3 we have that is, the spin of node 2 at a given time step will be driven to match the attractor state of node 1 at the previous time step. However, if j n 1~+ j c 1 and j n 2~+ j c 2 , J 21~0 . This gives s 2 (t)~z 1 with probability 1 = 2 In this case, node 2 receives no input from node 1. Nodes 1 and 2 have become effectively disconnected. This motivates new designations for node types. We define similarity nodes as nodes with j n i~j c i , and differential nodes as nodes with j n i~{ j c i . We also define the set of similarity nodes S~i : j n i~j c i È É and the set of differential nodes D~i : j n i È {j c i g. Connections between two similarity nodes or two differential nodes remain in the network, whereas connections that link nodes of different types transmit no signals. The effective deletion of edges between nodes means that the original network fully separates into two subnetworks: one composed entirely of similarity nodes (the similarity network) and another composed entirely of differential nodes (the differential network), each of which can be composed of one or more separate weakly connected components (see Fig. 1). With this separation, new source nodes (effective sources) can be exposed in both the similarity and differential networks. For the remainder of this article, Q is the set of both source and effective source nodes in a given network.

Control Strategies
The strategies presented below focus on selecting the best single nodes or small clusters of nodes to control, ranked by how much they individually change m a ? . In application, however, controlling many nodes is necessary to achieve a sufficiently changed m a ? .   The effects of controlling a set of nodes can be more than the sum of the effects of controlling individual nodes, and predicting the truly optimal set of nodes to target is computationally difficult.
Here, we discuss heuristic strategies for controlling large networks where the combinatorial approach is impractical. For both p~1 and p~2, simulating a cancer cell means that s s(0)~zj j c , and likewise for normal cells. Although the normal and cancer states are mathematically interchangeable, biologically we seek to decrease m c ? as much as possible while leaving m n ? &z1. By ''network control'' we thus mean driving the system away from its initial state ofs s(0)~j j c withh h ext . Controlling individual nodes is achieved by applying a strong field (stronger than the magnitude of the field due to the node's upstream neighbors) to a set of targeted nodes T so that This ensures that the drug field can always overcome the field from neighboring nodes.
In application, similarity nodes are never deliberately directly targeted, since changing their state would adversely affect both normal and cancer cells. Roughly 70% of the nodes in the networks surveyed are similarity nodes, so the search space is reduced. For p~2, the effective edge deletion means that only the differential network in cancer cells needs to be simulated to determine the effectiveness ofh h ext . For p~1, however, there may be some similarity nodes that receive signals from upstream differential nodes. In this case, the full effect ofh h ext can be determined only by simulating all differential nodes as well as any similarity nodes downstream of differential nodes. All following discussion assumes that all nodes examined are differential, and therefore targetable, for both p~1 and p~2. The existence of similarity nodes for p~1 only limits the set of targetable nodes.
Directed acyclic networks. Full control of a directed acyclic network is achieved by forcing s q~{ j c q for all q[Q. This guarantees m c ?~{ 1. Suppose that nodes q[Q in an acyclic network have always been fixed away from the cancer state, that is, Because there are no cycles present, all upstream paths of sufficent length terminate at a source. Because the spin of all nodes q[Q point away from the cancer attractor state, all nodes downstream must also point away from the cancer attractor state. Thus, for acyclic networks, forcing s q~{ j c q guarantees m c ?~{ 1. The complications that arise from cycles are discussed in the next subsubsection. However, controlling nodes in Q may not be the most efficient way to push the system away from the cancer basin of attraction and, depending on the control limitations, it may not be possible. If minimizing the number of controllers is required, searching for the most important bottlenecks is a better strategy.
Consider a directed network G and an initially identical copy, G 0~G : If removing node i (and all connections to and from i) from G 0 decreases the indegree of at least one node j[V (G 0 ), j=i, to less than half of its indegree in network G, (4) Remove all nodes j from the network G 0 . If additional nodes in G 0 have their indegree reduced to below half of their indegree in G, go to step 3. Otherwise, stop. The impact of the bottleneck i, I(i), is defined as where DX D is the cardinality of the set X : The impact of a bottleneck is the minimum number of nodes that are guaranteed to switch away from the cancer state when the bottleneck is forced away from the cancer state. Figure 1. Network segregation for two attractor states (p~2). Every edge that connects a similarity node to a differential node or a differential node to a similarity node transmits no signal. This means that the signaling in the right network shown above is identical to that of the left network. Because the goal is to leave normal cells unaltered while damaging cancer cells as much as possible, all similarity nodes can be safely ignored, and searches and simulations only need to be done on the differential subnetwork. doi:10.1371/journal.pone.0105842.g001 The impact is used to rank the size 1 bottlenecks by importance, with the most important as those with the largest impact. In application, when searching for nodes to control, any size 1 bottleneck fig that appears in the bottleneck control set of a different size 1 bottleneck fjg can be ignored, since fixing j to the normal state fixes i to the normal state as well. Note that the definition given above in terms of G and G 0 avoids miscounting in the impact of a bottleneck.
The network in Fig. 2, for example, has three sources (nodes 1, 2 and 3), but one important bottleneck (node 6). If full damage, i.e. m c ?~{ 1, is required, then control of all source nodes is necessary. If minimizing the number of directly targeted nodes is important and m c ? w{1 can be tolerated, then control of the bottleneck node 6 is a better choice.
Directed cycle-rich networks. Not all networks can be fully controlled at T~0 by controlling the source nodes, however. If there is a cycle present, paths of infinite length exist and the final state of the system may depend on the initial state, causing parts of the network to be hysteretic. Controlling only sources in a general directed network thus does not guarantee m c ?~{ 1 unless the system begins with s i~{ j c i . Define a cycle cluster, C, as a strongly connected subnetwork of a network G: The network in Fig. 3, for example, has one cycle cluster with nodes V (C)~4,5,6,7 f g . If the network begins with s s(0)~j j c , forcing both source nodes away from the cancer state does nothing to the nodes downsteam of node 3 (see Fig. 4). This is because the indegree deg { (4)~4, and a majority of the nodes connecting to node 4 are in the cancer attractor state. At T~0, cycle clusters with high connectivity tend to block incoming signals from outside of the cluster, resulting in an insurmountable activation barrier.
The most effective single node to control in this network is any one of nodes 4, 6 or 7. Forcing any of these away from the cancer attractor state will eventually cause the entire cycle cluster to flip away from the cancer state, and all nodes downstream will flip as well, as shown in Fig. 4. The cycle cluster here acts as a sort of large, hysteretic bottleneck. We now generalize the concept of bottlenecks.
Define a size k bottleneck in a network G to be a cycle cluster B with DV(B)D~k which, when removed from G, reduces the indegree of at least one node j[V (G), j[ =V B ð Þ to less than half of its original indegree. Other than now using the set of nodes V (B) rather than a single node set, the above algorithm for finding the bottleneck control set remains unchanged. In In general, however, more than one node in a cycle cluster may need to be targeted to control the entire cycle cluster. Fig. 5 shows a cycle cluster (composed of nodes 2-10) that cannot be controlled by targeting any single node. The precise value of n crit for a given cycle cluster C depends on its topology as well as the edges connecting nodes from outside of C to the nodes inside of C, and finding Z(C,G) can be difficult. We present a theorem that puts bounds on n crit to help determine whether a search for Z(C,G) is practical.
Theorem: Suppose a network G contains a cycle cluster C: Define the set of externally influenced nodes the set of intruder connections and the reduced set of critical nodes If N~DV (C)D and m: min where where Proof: First, prove the lower limit of Eq. 18. Let C be a cycle cluster in a network G with R(C,G)~f1g. (A cycle cluster in a network with DR(C,G)Dw0 will have the same or higher activation barrier for any node in the cluster than the same cycle cluster in a network with R~f1g. Since we are examining the lower limit of Eq. 18, we consider the case with the lowest activation barrier. Any externally influenced nodes cause n crit to either increase or remain the same.) For any node i to be able to flip away from the cancer state (although not necessarily remain there), we must have that h i~{ aj c i for a §0, meaning that at least half of the nodes upstream of i must point away from the cancer state. The node i requiring the smallest number of upstream nodes to be in the normal state is the node that satisfies deg { (i)~m: Controlling less than m=2 nodes will leave all uncontrolled nodes with a field in the cancer direction, and no more flips will occur. Thus, n crit §q m 2 r: For the upper limit of Eq. 18, consider a complete clique on N nodes, C~K N (that is, A ij~1 for all i, j[V (K N ), including self loops) in a network G: First, let there be no connections to any nodes in C from outside of C so that R(C,G)~f1g. For odd N, forcing (Nz1)=2 nodes away from the cancer state will result in the field for all nodes i: After one time step, all nodes will flip away from the cancer state. For even N, forcing N=2 nodes away from the cancer state will result in the field for all nodes i: At the next time step, the unfixed nodes will pick randomly between the normal and cancer state. If at least one of these nodes makes the transition away from the cancer state, the field at all other nodes will point away from the cancer direction. The system will then require one more time step to completely settle to s i~{ j c i :. Thus, we have that for C~K N in a network G with R(C,G)~f1g, K N with s i (0)~j c i gives the largest activation barrier for any cycle cluster on N nodes with R(C,G)~f1g to switch away from the cancer attractor state. A general cycle cluster C with any topology on N nodes with R(C,G)~f1g in a network G will have deg { (i)ƒN for all nodes i, and so we have the upper bound Finally, combining the upper limit in Eq. 26 with the lower limit from Eq. 20 gives Eq. 18. & There can be more than one Z red for a given cycle cluster. Note that the tightest constraints on n crit in Eq. 18 come from using the Z red with the largest overlap with R. If finding Z red is too difficult, an overestimate for the upper limit of n crit can be made by assuming that R\Z red~f 1g so that The cycle cluster in Fig. 5 has N~9, R~f2,9g, m~1, and one of the reduced sets of critical nodes is Z red~f 9,10g, so 1ƒn crit ƒ6: It can be shown through an exhaustive search that for this network n crit~2 , and the set of critical nodes is Z~f9,10g (see Fig. 6). Here, Z~Z red , although this is not always the case. Because the cycle cluster has 9 nodes and 1ƒn crit ƒ6, at most P 6 n~1 9 n

~4
65 simulations are needed to find at least one solution for Z(C,G). However, the maximum number of simulations required to find Z(C,G) increases exponentially and for larger networks the problem quickly becomes intractable. One heuristic strategy for controlling cycle clusters is to look for size k 0 vDV (C)D bottlenecks inside of C: Bottlenecks of size k&1 and average indegree Sdeg { (B)T%k can contain high impact size k 0 bottlenecks, where k 0 vk. Size k §1 bottlenecks need to be compared to find the best set of nodes to target to reduce m c ? . Simply comparing the impact is insufficent because a cycle cluster with a large impact could also have a large n crit , requiring much more effort than its impact merits. Define the critical efficiency of a bottleneck B as If the critical efficiency of a cycle cluster is much smaller than the impacts of size 1 bottlenecks from outside of the cycle cluster, the the cycle cluster can be safely ignored.
For some cycle clusters, however, not all of the nodes need to be controlled in order for a large portion of the nodes in the cycle cluster's control set to flip. Define the optimal efficiency of a bottleneck B as All strategies presented above are designed to select the best individual or small group of nodes to target. Significant changes in the biological networks' magnetization require targeting many nodes, however. Brute force searches on the effect of larger combinations of nodes are typically impossible because the required number of simulations scales exponentially with the number of nodes. A crude Monte Carlo search is also numerically expensive, since it is difficult to sample an appreciable portion of the available space. One alternative is to take advantage of the bottlenecks that can be easily found, and rank all size k §1 bottlenecks B i in an ordered list U such that where for all B i ,B j [U and fix the bottlenecks in the list in order. This is called the efficiency-ranked strategy. If all size kw1 bottlenecks are ignored, it is called the pure efficiency-ranked strategy, and if size kw1 bottlenecks are included it is called the mixed efficiencyranked strategy. An effective polynomial-time algorithm for finding the top z nodes to fix, which we call the best+1 strategy (equivalent to a greedy algorithm), works as follows: (1) Begin with a seed set of nodes to fix, F ; (2) Test the effect of fixing F |i for all allowed nodes i[ =F ; (3) F /F |i best , where i best is the best node from all i sampled; (4) If DF Dvz, go to step (2). Otherwise, stop. The seed set of nodes could be the single highest impact size 1 bottleneck in the network, or it could be the best set of n nodes (where nvz) found from a brute force search.

Cancer Signaling
In application to biological systems, we assume that the magnetization of cell type a is related to the viability of cell type a, that is, the fraction of cells of type a that survives a drug treatment. It is reasonable to assume that the viability of cell type a, v a (m a ? ), is a monotonically increasing function of m a ? . Because the exact relationship is not known, we analyze the effect of external perturbations in terms of the final magnetizations.
We need to use as few controllers as possible to sufficiently reduce m c ? while leaving m n ? &z1. In practical applications, however, one is limited in the set of druggable targets. All classes of drugs are constrained to act only on a specific set of biological components. For example, one class of drugs that is currently under intense research is protein kinase inhibitors [35]. In this case one has two constraints: the only nodes that can be targeted are those that correspond to kinases, and they can only be inhibited, i.e. turned off. We will use the example of kinase inhibitors to show how control is affected by such types of constraints. In the real systems studied, many differential nodes have only similarity nodes upstream and downstream of them, while the remaining differential nodes form one large cluster. This is not important for p~1, but the effective edge deletion for p~2 results in many Figure 5. A network with a cycle cluster C, composed of nodes 2-10, that cannot be controlled at T~0 by controlling any single node. Here, the set of externally influenced nodes is R(C,G)~f2,9g, the set of intruder connections is W (C,G)~f(1,2), (1,9)g, the reduced set of critical nodes is Z red (C,G)~f9,10g, the minimum indegree is m~1 and the number of nodes in the cycle cluster is N~9: By Eq. 18, this gives the bounds of the critical number of nodes to be 1ƒn crit ƒ6. doi:10.1371/journal.pone.0105842.g005 Figure 6. Magnetization for network from Fig. 5, averaged over 10,000 runs. There is no single node to target that will control the cycle cluster, but fixing nodes 9 and 10 results in full control of the cycle cluster, leaving only node 1 in the cancer state. This means Z(C,G)~f9,10g and n crit~2 . doi:10.1371/journal.pone.0105842.g006 islets, which are nodes i with A ij~Aji~0 for all i=j (self-loops allowed). Controlling islets requires targeting each islet individually. For p~2, we concentrate on controlling only the largest weakly connected differential subnetwork. All final magnetizations are normalized by the total number of nodes in the full network, even if the simulations are only conducted on small portion of the network.
The data files for all networks and attractors analyzed below can be found in Supporting Information.

Lung Cell Network
The network used to simulate lung cells was built by combining the kinase interactome from PhosphoPOINT [32] with the transcription factor interactome from TRANSFAC [33]. Both of these are general networks that were constructed by compiling many observed pairwise interactions between components, meaning that if j?i, at least one of the proteins encoded by gene j has been directly observed interacting with gene i in experiments. This bottom-up approach means that some edges may be missing, but those present are reliable. Because of this, the network is sparse (*0:057% complete, see Table 2), resulting in the formation of many islets for p~2. Note also that this network presents a clear hierarchical structure, characteristic of biological networks [36,37], with many ''sink'' nodes [38] that are targets of  transcription factors and a relatively large cycle cluster originating from the kinase interactome. It is important to note that this is a non-specific network, whereas real gene regulatory networks can experience a sort of ''rewiring'' for a single cell type under various internal conditions [39]. In this analysis, we assume that the difference in topology between a normal and a cancer cell's regulatory network is negligible. The methods described here can be applied to more specialized networks for specific cell types and cancer types as these networks become more widely avaliable. In our signaling model, the IMR-90 cell line [40,41] was used for the normal attractor state, and the two cancer attractor states examined were from the A549 (adenocarcinoma) [42][43][44][45][46] and NCI-H358 (bronchioalveolar carcinoma) [42,43] cell lines. Gene expression measurements from all referenced studies for a given cell line were averaged together to create a single attractor. The resulting magnetization curves for A549 and NCI-H358 are very similar, so the following analysis addresses only A549. The full network contains 9073 nodes, but only 1175 of them are differential nodes in the IMR-90/A549 model. In the unconstrained p~1 case, all 1175 differential nodes are candidates for targeting. Exhaustively searching for the best pair of nodes to control requires investigating 689725 combinations simulated on the full network of 9073 nodes. However, 1094 of the 1175 nodes are sinks (i.e. nodes i with outdegree deg z (i)~0, ignoring self loops) and therefore have I(i)~e opt (i)~1, which can be safely ignored. The search space is thus reduced to 81 nodes, and finding even the best triplet of nodes exhaustively is possible. Including constraints, only 31 nodes are differential kinases with j c i~z 1. This reduces the search space at the cost of increasing the minimum achievable m c ? . There is one important cycle cluster in the full network, and it is composed of 401 nodes. This cycle cluster has an impact of 7948 for p~1, giving a critical efficiency of at least *19:8, and 1ƒn crit ƒ401 by Eq. 27. The optimal efficiency for this cycle cluster is e opt~2 9, but this is achieved for fixing the first bottleneck in the cluster. Additionally, this node is the highest impact size 1 bottleneck in the full network, and so the mixed efficiency-ranked results are identical to the pure efficiency-ranked results for the unconstrained p~1 lung network. The mixed efficiency-ranked strategy was thus ignored in this case. Fig. 7 shows the results for the unconstrained p~1 model of the IMR-90/A549 lung cell network. (All simulations were performed using MATLAB on a desktop computer. Running the simulations to make all curves shown below required approximately 12 hours.) The unconstrained p~1 system has the largest search space, so the Monte Carlo strategy performs poorly. The best+1 strategy is the most effective strategy for controlling this network. The seed set of nodes used here was simply the size 1 bottleneck with the largest impact. Note that best+1 works better than effeciency-ranked. This is because best+1 includes the synergistic effects of fixing multiple nodes, while efficiency-ranked assumes that there is no overlap between the set of nodes downstream from multiple bottlenecks. Importantly, however, the efficiency-ranked method works nearly as well as best+1 and much better than Monte Carlo, both of which are more computationally expensive than the efficiency-ranked strategy. Fig. 8 shows the results for the unconstrained p~2 model of the IMR-90/A549 lung cell network. The search space for p~2 is much smaller than that for p~1. The largest weakly connected differential subnetwork contains only 506 nodes (see Table 3), and the remaining differential nodes are islets or are in subnetworks composed of two nodes and are therefore unnecessary to consider. Of these 506 nodes, 450 are sinks. Fig. 9 shows the largest weakly connected component of the differential subnetwork, and the top five bottlenecks in the unconstrained case are shown in red. If limiting the search to differential kinases with j c i~z 1 and ignoring all sinks, p~2 has 19 possible targets. There is only one cycle cluster in the largest differential subnetwork, containing 6 nodes. Like the p~1 case, the optimal efficiency occurs when targeting the first node, which is the highest impact size 1 bottleneck. Because the mixed efficiency-ranked strategy gives the same results as the pure efficiency-ranked strategy, only the pure strategy was examined. The Monte Carlo strategy fares better in the unconstrained p~2 case because the search space is smaller. Additionally, the efficiency-ranked strategy does worse against the best+1 strategy for p~2 than it did for p~1. This is because the effective edge deletion decreases the average indegree of the network and makes nodes easier to control indirectly. When many upstream bottlenecks are controlled, some of the downstream bottlenecks in the efficiency-ranked list can be indirectly controlled. Thus, controlling these nodes directly results in no change in the magnetization. This gives the plateaus shown for fixing nodes 9-10 and 12-15, for example.
The only case in which an exhaustive search is possible is for p~2 with constraints, which is shown in Fig. 10. Note that the polynomial-time best+1 strategy identifies the same set of nodes as the exponential-time exhaustive search. This is not surprising, however, since the constraints limit the available search space. This means that the Monte Carlo also does well. The efficiencyranked method performs worst. The efficiency-ranked strategy is designed to be a heuristic strategy that scales gently, however, and is not expected to work well in such a small space when compared with more computationally expensive methods.

B Cell Network
The B cell network was derived from the B cell interactome of Ref. [34]. The reconstruction method used in Ref. [34] removes edges from an initially complete network depending on pairwise gene expression correlation. Additionally, the original B cell network contains many protein-protein interactions (PPIs) as well as transcription factor-gene interactions (TFGIs). TFGIs have definite directionality: a transcription factor encoded by one gene affects the expression level of its target gene(s). PPIs, however, do not have obvious directionality. We first filtered these PPIs by checking if the genes encoding these proteins interacted according to the PhosphoPOINT/TRANSFAC network of the previous section, and if so, kept the edge as directed. If the remaining PPIs are ignored, the results for the B cell are similar to those of the lung cell network. We found more interesting results when keeping the remaining PPIs as undirected, as is discussed below.
Because of the network construction algorithm and the inclusion of many undirected edges, the B cell network is more dense (*0.290% complete, see Table 2) than the lung cell network. This Table 3. Properties of the largest weakly connected differential subnetworks for all cell types. higher density leads to many more cycles than the lung cell network, and many of these cycles overlap to form one very large cycle cluster containing *66% of nodes in the full network. All gene expression data used for B cell attractors was taken from Ref. [47]. We analyzed two types of normal B cells (naïve and memory) and three types of B cell cancers (diffuse large B-cell lymphoma (DLBCL), follicular lymphoma, and EBV-immortalized lymphoblastoma), giving six combinations in total. We present results for only the naïve/DLBCL combination below, but Tables 3 and 4 list the properties of all normal/cancer combinations. Again, all gene expression measurements for a given cell type were averaged together to produce a single attractor. The full B cell network is composed of 4364 nodes. For p~1, there is one cycle cluster C composed of 2886 nodes. This cycle cluster has 1ƒn crit (C) ƒ1460, I(C)~4353, and 3:0ƒe crit (C)ƒ4353: Finding Z(C) was deemed too difficult. Fig.11 shows the results for the unconstrained p~1 case. Again, the pure efficiency-ranked strategy gave the same results as the mixed efficiency-ranked strategy, so only the pure strategy was analyzed. As shown in Fig. 11, the Monte Carlo strategy is outperformed by both the efficiency-ranked and best+1 strategies. The synergistic effects of fixing multiple bottlenecks slowly becomes apparent as the best+1 and efficiency-ranked curves separate. Fig. 12 shows the results for the unconstrained p~2 case. The largest weakly connected subnetwork contains one cycle cluster with 351 nodes, with 1ƒn crit ƒ208. Although finding a set of critical nodes is difficult, the optimal efficiency for this cycle cluster is 62.2 for fixing 10 bottlenecks in the cycle cluster. This makes targeting the cycle cluster worthwhile. The efficiency of this set of 10 nodes is larger than the efficiencies of the first 10 nodes from the pure efficiency-ranked strategy, so the m c ? from the mixed strategy drops earlier than the pure strategy. Both strategies quickly identify a small set of nodes capable of controlling a significant portion of the differential network, however, and the same result is obtained for fixing more than 10 nodes. The best+1 strategy finds a smaller set of nodes that controls a similar fraction of the cycle cluster, and fixing more than 7 nodes results in only incremental decreases in m c ? . The Monte Carlo strategy performs poorly, never finding a set of nodes adequate to control a significant fraction of the nodes in the cycle cluster.

Conclusions
Signaling models for large and complex biological networks are becoming important tools for designing new therapeutic methods for complex diseases such as cancer. Even if our knowledge of biological networks is incomplete, rapid progress is currently being made using reconstruction methods that use large amounts of publicly available omic data [12,13]. The Hopfield model we use in our approach allows mapping of gene expression patterns of normal and cancer cells into stored attractor states of the signaling dynamics in directed networks. The role of each node in disrupting the network signaling can therefore be explicitly analyzed to identify isolated genes or sets of strongly connected genes that are selective in their action. We have introduced the concept of size k bottlnecks to identify such genes. This concept led to the formulation of several heuristic strategies, such as the efficiencyranked and best+1 strategy to find nodes that reduce the overlap of the cell network with a cancer attractor. Using this approach, we have located small sets of nodes in lung and B cancer cells which, when forced away from their initial states with local magnetic fields (representing targeted drugs), disrupt the signaling of the cancer cells while leaving normal cells in their original state. For networks with few targetable nodes, exhaustive searches or Monte Carlo searches can locate effective sets of nodes. For larger networks, however, these strategies become too cumbersome and our heuristic strategies represent a feasible alternative. For tree-like networks, the pure efficiency-ranked strategy works well, whereas the mixed efficiency-ranked strategy could be a better choice for networks with high-impact cycle clusters.
We make two important assumptions in applying this analysis to real biological systems. First, we assume that genes are either fully off or fully on, with no intermediate state. Modelling the state of a neuron as ''all-or-none'' has long been accepted as a reasonable assumption [48], which provided the spin glass framework for the   Hopfield model. While similar switch-like behavior in gene regulatory networks has been proposed as an explanation of, for example, segmentation in Drosophila embryos [49], assigning a Boolean value to gene expression may be overly simplistic in many cases. A model which uses spins with more than two projections could prove to be more realistic and predictive. Second, we assume that all nodes update their status with a single timescale and with a single interaction strength. If the signaling timescale t ij for each edge in the biological network is known, simulations could be conducted in which a signal traveling along an edge (j,i) reaches its target after t ij time steps. This would amount to a synchronous update schedule with a ''queue'' of signals moving between nodes. Likewise, our model gives equal weight to all edges (aside from edges that are effectively deleted in the p~2 case), whereas real gene regulatory networks exhibit a spectrum of interaction strengths. This can easily be integrated with our model by using a weighted, directed adjacency matrix. However, doing this would surely require a change in control strategy. Despite these issues, our model shows promise. Some of the genes identified in Table 4 are consistent with current clinical and cancer biology knowledge. For instance, in the lung cancer list we found a well known tumor suppressor gene (TP53) [50] that is frequently mutated in many cancer types including lung cancer [51]. Mutations in PBX1 have recently been detected in nonsmall-cell lung cancer and this gene is now being considered as a target for therapy and prognosis [52]. MAP3K3 and MAP3K14 are in the MAPK/ERK pathway which is a target of many novel therapeutic agents [53], and SRC is a well known oncogene and a candidate target in lung cancer [54]. BCL6 (B-cell lymphoma 6) is the most common oncogene in DLBCL, and it is known that its expression can predict prognosis and response to drug therapy [55][56][57]. BCL6 is also frequently mutated in follicular lymphoma [58,59]. Our analysis identified BCL6 as an important drug target for both DLBCL and follicular lymphoma using either naive or memory B-cells as a control for both p~1 and p~2. RBL2 disregulation has been recently associated with many types of lymphoma [60][61][62]. FOXM1 is a potential therapeutic target in mature B cell tumors [63] and ATF2 has been recently found to be highly disregulated in lymphoma [64,65]. Besides BCL6 discussed above, the N/D list for DLBCL contains genes (MEF2A [66], NCOA1 [67,68], TGIF1 [69][70][71], NFATC3 [72]) that are all known to have a functional role in cancer, even if they have not been associated to the specific B-cell cancer types we have considered. Our predictions are for the immortalized cell lines we have selected, some of which are commonly used for in-vitro testing in many laboratories. RNAi and targeted drugs could then be used in these cell lines against the top scoring genes in Table 4 to test the disruption of survival or proliferative capacity. If experimentally validated, our analysis based on attractor states and bottlenecks could be applied to patient-derived cancer cells by integrating in the model patient gene expression data to identify patient-specific targets.
The above unconstrained searches assume that there exists some set of ''miracle drugs'' which can turn any gene ''on'' and ''off'' at will. This limitation can be patially taken into account by using constrained searches that limit the nodes that can be addressed. However, even the constrained search results are unrealistic, since most drugs directly target more than one gene. Inhibitors, for example, could target differential nodes with j c i~{ 1 and j n i~z 1, which would damage only normal cells. Additionally, drugs would not be restricted to target only differential nodes, and certain combinations could be toxic to both normal and cancer cells. Few cancer treatments involve the use of a single drug, and the synergistic effects of combining multiple drugs adds yet another level of complication to finding an effective treatment [27]. On the other hand, the intrinsic nonlinearity of a cellular signaling network, with its inherent structure of attractor states, enhances control [31] so that a properly selected set of druggable targets might be sufficient for robust control.

Supporting Information
Table S1 Lung cell network. The column labeled ''Source EzID'' contains the Entrez IDs of transcription factors and kinases, and ''Target EzID'' contains the Entrez IDs of the genes targeted by the transcription factor or kinase to its left. (TXT) Table S2 IMR-90/A549 attractors for lung cell network. The column labeled ''EzID'' contains the Entrez ID of the genes. The second and third columns are the normal and cancer attractor, respectively. (TXT) Table S3 IMR-90/NCI-H358 attractors for lung cell network. The column labeled ''EzID'' contains the Entrez ID of the genes. The second and third columns are the normal and cancer attractor, respectively.     contains the Entrez ID of the genes. The second and third columns are the normal and cancer attractor, respectively. (TXT) Table S10 Naïve/follicular lymphoma attractors for B cell network. The column labeled ''EzID'' contains the Entrez ID of the genes. The second and third columns are the normal and cancer attractor, respectively. (TXT)