Structure-based approach to identifying small sets of driver nodes in biological networks

Networks are one of the models used to encode the intricate interactions that underlie complex systems.^1,2 A network representation uses nodes to denote the components of the system and edges to represent their interactions. In some contexts, the only information given about a complex system is the network representation with its nodes and edges, which can be studied to gain valuable information.^3–5 In other cases, dynamics are also specified on the network to describe the processes that take place in the system. These dynamic processes are described by assigning a time-dependent variable (continuous or discrete) to each node. We can represent the states of the nodes of our system as the state vector $X(t)={x1(t),…,xn(t)}$⁠. For each node in the network, we represent its time evolution as $xi(t+1)=Fi(XIi(t))$⁠, where $Fi$ is the node update function for node $ni$⁠, specifying the state of the node in the next discrete time step. It is governed by time $t$ and $XIi(t)$⁠, which is a state vector of only the nodes that have edges directed at node $ni$⁠, i.e., they are the inputs $Ii$ of node $ni$⁠. The set $Ii$ is empty for source nodes, and the evolution of the source nodes is fully determined by an independent function. The set of node update equations are important to grasp how the system’s state changes in time.

In addition to studying complex systems’ autonomous dynamics, understanding how to influence these dynamics to reach important states is of great interest. Network control has many theoretical and practical facets to determine the best course of action for controlling a complex system.^6–17 Different control methods are applicable based on the information available (e.g., whether dynamic information is available or not) and on the goal (e.g., which state we want to reach). From here, control methods aim to find a set of elements and the corresponding actions on them required to attain the desired objective and then evaluate the practicality of this set.

We can apply these ideas to biology by representing biological systems as networks and modeling biological processes as dynamic information propagation on these networks. In biological networks at the cellular level, nodes represent biological macromolecules and edges represent the interactions between these macromolecules. These interactions underlie cell phenotypes and behaviors such as movement, cell division, and programed cell death. Dynamical models on a cellular network are used to describe the signal transduction, gene regulation, and metabolic processes of the biological system.^18–22 Thus, the models’ dynamics must capture the long-term steady states—attractors—of the cell (e.g., cell types).

When attempting to control complex biological systems, reaching an important attractor is often the main objective. As an example, the cell network in Fig. 1 represents a simplified version of the signaling network that underlies T cell large granular lymphocyte (T-LGL) leukemia.²³ There are two attractors, which represent two cell fates: a state of commitment to cell death (apoptosis) and a state of cell survival, which is pathological in this case. Understanding how to drive the system to the apoptosis attractor leads to predictions of therapeutic interventions.

One method for driving a network to any of its natural attractors (attractor control) is feedback vertex set (FVS) control.^13,14 In this control method, we assume that we can override the state of the variables of interest into a target state. FVS control depends on the connection between the feedback vertex set (FVS) and the attractors of the system. The FVS is a set of nodes in the network that contains nodes in every cycle of the network. Mochizuki et al. proved that, given a system governed by a general class of nonlinear dynamics on an underlying network structure, driving the state of every node in the FVS to its corresponding state in a target attractor is guaranteed to drive the system to this target attractor.¹³ Specifically, Mochizuki et al. demonstrated that the minimal FVS is a minimal set of “determining nodes” and mathematically proved this result. “Determining nodes” of the network are a subset of the nodes $J⊆{n1,…,nN}$ such that the convergence of the nodes in $J$ causes the convergence of every node’s trajectories. That is, if there are two attractors $X(t)$ and $X~(t)$ such that $XJ(t)−XJ~(t)→0$ in the limit $t→∞$⁠, i.e., for both attractors every node in $J$ converges to the same trajectory, then $X(t)−X~(t)→0$ in the limit $t→∞$⁠. Mochizuki et al. considered networks with no source nodes (nodes that do not receive any internal regulatory input). Zañudo et al. later showed that if the source nodes are not uniquely defined they must also be driven into their corresponding states in the target attractor.¹⁴ Source nodes typically signify an environmental signal or a receptor of that signal, so in a specific environment their states are uniquely defined and not necessary for driving the system to its attractors.

For our purposes, we consider the minimal FVS, which is the smallest set of nodes that accomplishes this task. This set is not necessarily unique. In Fig. 1(a), we label a possible minimal FVS in with dashed outlines; however, Ceramide is part of the same cycles as S1P, so it could be part of the minimal FVS instead of S1P. Similarly, DISC could replace FLIP in the minimal FVS. Figure 1 focuses on the minimal FVS of S1P and FLIP, but all minimal FVSs are probed in our broader analysis.

To better understand the connection between the FVS and the attractors of the system, it is useful to relate it to previous work connecting cycles in the network structure to the multistability in the system. Because the FVS consists of nodes in every feedback in the network, removing the FVS creates an acyclic network. As demonstrated by René Thomas, the cycles of a network are directly related to its multistability,^24,25 so this acyclic network only has one attractor. Driving a node in a cycle to a specific state will remove the multistability of that cycle. Thus, overriding the state of every node in the FVS into their states in a specific attractor reduces the state space of the network so there is a single attractor, namely, the desired one.¹³ In our biological systems, the act of overriding the FVS nodes can be accomplished through many means, most prominently through gene knockout or gene knock-in. These actions do not give full control over the state of the nodes, which is a difficult problem that is still being studied,⁷ but they do drive the gene’s concentration to its extremes, which will cause a response that reflects the process of overriding the node’s state in our dynamical models.

Controlling the entire FVS is a sufficient condition to drive the network into a given attractor, but it is not always necessary. Previous work has shown that control sets smaller than the FVS can drive the trajectory of a system into a desired attractor,^9,13,14,26 which means that partial control of the FVS (e.g., control of a subset of the FVS) can be sufficient for attractor control when restricted to a specific attractor or to a particular model.²⁷ This is important because the size of the FVS can often be larger than the size of combinatorial interventions that can be implemented in biological experiments, which is usually restricted to 1–3 nodes. For example, the gene regulatory network in ascidian embryos studied by Kobayashi et al. has a FVS size of five, and it is near the limit of what can be currently controlled experimentally.²⁸ Thus, it is of great interest to find FVS subsets that can drive a network to its natural attractors. At present, it is not known how to most effectively choose the correct FVS subsets, so in this work, we aim to rank various FVS subsets based on their ability to drive the system into a desired attractor.

Although the FVS can be found without knowledge of the dynamics of the system, in this work, we focus on Boolean dynamic models to evaluate the effectiveness of our methods. Boolean models characterize each node with two values, usually referred to as OFF (0) and ON (1) so that each node’s state becomes $xi(t)∈{0,1}$⁠. The update functions then become $Fi:{0,1}N→{0,1}$⁠, where $Fi$ is usually based on logical operations [Fig. 1(b)]. We evaluate the dynamics of these systems using a general asynchronous update scheme, where at each discrete time step, we stochastically choose a node $nj$ and update its state according to its update function, which makes the state vector $X(t+1)={x1(t),…,Fj(XIj(t)),…,xN(t)}$⁠. In biology, it is conventional that the node name (e.g., the molecule that the node represents) is used to denote the node’s state. While Boolean models represent an approximation, they have been shown to capture the attractors and qualitative dynamics of various systems.^18–20 Boolean models have a larger size (include more components) than continuous models. Larger networks will have larger FVSs, which will allow us to study cases in which a small fraction of the FVS is controlled. Furthermore, because the state space of Boolean models is finite, it is possible to comprehensively probe the effects of an intervention on the state space of the network. Thus, Boolean models provide a strong test bed for evaluating our methods.

The reduced T-LGL network in Fig. 1(a) was described with a Boolean model by Saadatpour et al. [Fig. 1(b)]²³ The network’s two fixed-point attractors— apoptosis and survival—are represented in Fig. 1(c), where white represents a node that is OFF and dark gray represents a node that is ON. We know controlling both of the minimal FVS nodes, S1P and FLIP, can drive the network into either of its attractors, but how does controlling only one of these two nodes (indicated by an dashed-dotted outline) affect the dynamics? Controlling S1P alone can still drive the system to either of the two attractors [Fig. 1(d)]. However, controlling FLIP alone does not fix the state of any other node and both original attractors are still possible [Fig. 1(e)]. This difference between FVS members motivated us to search for an answer to the following questions: Are there network (topological) measures that differentiate S1P from FLIP and thus can predict that S1P outperforms FLIP in driving the system into an attractor? In a general network, how can we identify which FVS subsets are most likely to control the network?

We evaluated the ability of seven topological measures to identify important FVS subsets. These measures capture a node’s influence on the network, determined by either the local or global interactions of the node within the network. This influence tells us how much we expect the node to control the other nodes in the network. The seven metrics are out-degree, average inverse distance, PRINCE,^29,30 a modified version of PRINCE, CheiRank,³¹ involvement in the strongly connected component (SCC), and involvement in the network’s positive cycles. These metrics fall into three categories: centrality measures, propagation measures, and cycle-based measures.

The centrality measures, i.e., out-degree (a local measure) and average inverse distance (a global measure) provide a baseline of how important each node is [see the  Appendix and Fig. 2(a)]. PRINCE propagates a fixed perturbation of a target node through the network in a similar way as mass flow; our modification relaxes this assumption. CheiRank is a measure of the “source-ness” of each node [see the  Appendix and Fig. 2(b)]. The cycle-based metrics aim to measure how much a subset contributes to the multistability of the network [see the  Appendix and Fig. 2(c)].

When we introduce an intervention on a subset of the FVS, we measure the intervention’s effect on the attractors and their basins to understand how well it controls the system. First, based on the intervention, we sort the attractors into targets and non-targets [Fig. 3(b)]. Applying the intervention to the system restricts the state space of the system [Fig. 3(c)], which changes the basins of attraction of the system’s original attractors and may also introduce new attractor(s) [Fig. 3(d)]. When the goal is to drive the system into a target attractor, a successful intervention increases the basin of attraction of the target attractor and decreases the basins of the other attractors. We propose two metrics that measure these two aspects of control. To Control is the measure of how well control of a FVS subset drives the network to the target attractor whereas Away Control measures how well the same subset drives the network away from non-targeted attractors. The basis of both metrics is the normalized percent change in the relevant basin of attraction compared to the unperturbed system (see the  Appendix). These two values may be equal but are not necessarily the same. For example, if an intervention decreases the basin of a non-target attractor while driving the system to a new attractor, it will not increase the basin of the target attractor, so the intervention will have a large Away Control but a small To Control.

As our aim is to quantify the effectiveness of the FVS subsets identified from each of our topological metrics, we perform a systematic evaluation of all interventions, i.e., for a FVS subset of size $L$⁠, we consider all $2L$ combinations of fixing each node in the ON or OFF state. For each of these $2L$ interventions, we sort the attractors into target and non-target attractor(s). Each attractor wherein the states of every node are consistent with the intervention states becomes a member of the target attractors, while the other attractors are classified as non-target. Interventions consistent with every attractor are not informative, that is, an intervention that classifies every attractor as a target does not give any valuable information. Similarly, interventions that are not consistent with any attractors are only partially informative, that is, an intervention that classifies every attractor as a non-target can only give us information on how well the intervention can drive the system away from the attractors of the system (see the  Appendix). We refer to interventions that classify some attractors as targets and others as non-target attractors (so both sets are non-empty) as fully informative interventions. For every fully informative intervention, we determine the values of To Control and Away Control; if the intervention is partially informative, we only determine the value of Away Control.

For each FVS subset of size $L$⁠, there are $2L$ interventions, each with a value for To Control and Away Control. To summarize the values of all interventions associated with a FVS subset with a single quantity, we define an aggregate value of To Control and Away Control for a FVS subset as the maximal value over all its (partially) informative interventions. Once calculated, the value of these control metrics (control values) are compared with the topological metrics to determine how effective each metric is at identifying subsets that can control the network. We define a FVS subset as successful if the aggregate value of the control metric of interest (To Control or Away Control) is larger than a success threshold, which we choose to be 0.9.

We relate each of the seven topological metrics to the two control metrics using a logistic regression [Fig. 4(a)]. After creating the logistic regression, the sorting strength of each topological metric when classifying successful FVS subsets is determined based on the area under the precision–recall curve (AUPRC) metric [see the  Appendix and Fig. 4(b)]. If the AUPRC value is equal to the fraction of positive data points, the topological metric cannot sort the FVS subsets better than randomly sorting subsets. We define the AUPRC predictive threshold as the AUPRC value that is greater than halfway between the fraction of positive data points and the maximum AUPRC value of 1. Thus, a regressor variable is predictive if it sorts the binarized values of a control metric with an AUPRC value above the AUPRC predictive threshold. For example, if the fraction of positive data points is 0.75, the AUPRC predictive threshold is 0.875, so an AUPRC greater than 0.875 would be considered predictive. Calculating these AUPRC values for each of our topological metrics across different networks and subset sizes allows us to determine how well they can identify successful FVS subsets, or in other words, how well they predict a FVS subset’s ability to drive the system into the target attractor(s) or away from the non-target attractor(s).

We investigated an array of four different networks, corresponding to well-established Boolean models of biological systems, namely, the T-LGL network,³² NSCLC network,³³ an FA/BRCA variant network,³⁴ and a Helper T Cell Differentiation network.³⁵ For each of our networks, we determine the states of the source nodes. Because the source nodes are uniquely defined, we are only concerned with overriding the nodes within the FVS of each network and do not need to consider the full FC control presented in Zañudo et al.¹⁴ This is done to remove any effect the source nodes have on the attractor trajectory, so we can study only the FVS node’s effect on driving the system. For these networks, we need to identify their minimal FVS, which are not necessarily unique, but identifying a minimal FVS is an NP-hard problem. Instead, and following Zañudo et al.,¹⁴ we identify near-minimal FVSs using one of the efficient methods available (the simulated annealing algorithm introduced by Galinier et al.³⁶) We analyzed all FVS subsets of one to three nodes from multiple near-minimal FVSs we identified. The size of these networks range from 28 to 60 nodes, and the size of the near-minimal FVSs range from 8 to 17 nodes (see Table S4 in the supplementary material).

To characterize the predictive power of each of the seven topological metrics, we determined their associated AUPRC values for each network and FVS subset size and used these AUPRC values to determine the predictive power of each topological metric in the task of classifying successful control subsets. We summarize these results in Table I.

The top of Table I shows the percentage of cases (out of 12 total cases, where each case corresponds to a single FVS subset size of one of the networks) in which the indicated metric was predictive according to the AUPRC. Propagation metrics had the highest % of cases of all seven metrics for both To Control and Away Control and cycle-based metrics had the lowest % of cases for both control metrics. The bottom of Table I shows the rank of the AUPRC values for each of the seven metrics averaged over all 12 cases, where a lower rank corresponds to a higher AUPRC value. Consistent with what we found with the % of cases, propagation metrics had the lowest average rank of all seven metrics for both control metrics. Although the order of centrality and cycle-based metrics was not consistent for both control metrics, cycle-based metrics had the highest average rank of all metrics for the To Control metric. In summary, we found that propagation metrics had the highest predictive power out of all metrics (highest % of cases in which they were predictive and lowest average AUPRC rank) and that cycle-based metrics had the lowest predictive power out of all metrics (lowest % of cases in which they were predictive and highest average AUPRC rank for the To Control metric).

The number of FVS subsets that have a high value in a topological metric but a low control value reflects the metric’s inability to fully capture the dynamics of the system. We hypothesized that intersections of the FVS subsets predicted to successfully control the network by each topological metric would reduce the number of FVS subsets that have a high value in a topological metric but a low control value. To do this intersection, we set a fixed percentile cutoff, and for each topological metric, we identified the set of FVS subsets that was above the percentile cutoff (see the  Appendix and Fig. 5). The sets of above-cutoff FVS subsets [box in Fig. 5(a)] for each topological metric of interest were then intersected with each other resulting in a set of FVS subsets that are hypothesized to better control the network and contain a lower percentage of false positives than each metric individually. We use the term intersection metrics to refer to the metrics defined by intersecting topological metrics, and we denote them according to the metrics that are included in the intersection. In an intersection metric, we assign to each FVS subset the highest percentile cutoff value at which the FVS subset appears in the intersected topological metrics (see the  Appendix). To evaluate how well an intersection metric performs as more metrics are included in the intersection and to determine which combinations of the most predictive metrics best approximate FVS subsets’ level of control over the network dynamics, we tested three different intersection metrics: (1) an intersection of all seven metrics, (2) an intersection of just the propagation metrics, and (3) an intersection of the Modified PRINCE and CheiRank metrics. As an example, Fig. 6(a) shows the percentile cutoff value and Away Control value for the FVS subsets of the intersection metric that includes the three propagation metrics.

To analyze each intersection metric, we took a similar approach to what we did in Sec. III A. We obtained the AUPRC values for each intersection metric, the percentage of AUPRC values above the predictive threshold (Table II top), and the average rank of the AUPRC values among the three selected intersection metrics and all seven topological metrics (Table II bottom). The intersection metrics are predictive (have an AUPRC above the predictive threshold) in more cases than the non-propagation metrics, as reflected by their higher percentage of cases (see Table I for the non-propagation metrics and Table II for the intersection metrics). The intersection metrics always have an average rank that is lower than that of the non-propagation metrics (Table S1), and they perform similarly to the propagation metrics (Table II), which performed the best among all individual topological metrics (Table I). Among the intersection metrics, the Modified PRINCE and CheiRank intersection metric is predictive in the highest percentage of cases (75% in both To Control and Away Control), it has the lowest average rank for Away Control, and it has the second lowest average rank for To Control. Notably, and contrary to our expectations, the intersection metrics do not outperform the Modified PRINCE metric, which has the first or second highest percentage of cases in which it is predictive and the lowest average rank of AUPRC values among all metrics (Table II).

Although the Modified PRINCE and CheiRank intersection metric performed the best among the intersection metrics and the propagation intersection metric performed the worst, we noticed that relying on the AUPRC to rank metrics can mask a property that is desirable in our setting: the ability of a topological metric to be predictive in the high precision and low recall regime. This regime is important because it is where we expect topological metrics to excel. For each FVS subset we expect that having a top rank in a predictive topological metric is sufficient but not necessary to have a high control metric value; this corresponds to the high precision and low recall regime. Note that the reason we do not necessarily expect topological metrics to fully capture the control metric values is that the control values depend on the dynamic information in the model and not solely on the network topology.

For example, in the two-node, the Away Control case for the NSCLC network, the AUPRC for the propagation intersection metric is 0.669, which is below the AUPRC predictive threshold of 0.698, but on the precision–recall curve [Fig. 6(c)], the propagation intersection metric has the highest recall before the precision drops from 1.0, i.e., the propagation intersection metric identifies the most successful Away Control FVS subsets before encountering an unsuccessful Away Control FVS subset. Despite performing the worst according to the AUPRC values, the propagation intersection metric performs the best among all studied metrics (has the lowest average rank, 2.83) when comparing metrics according to the number of the successful FVS subsets they identify before the precision drops below 0.95 (supplementary material, Table S3 in the supplementary material). The propagation intersection metric outperforms Modified PRINCE, which performed the second best (2.96), and the Modified PRINCE and CheiRank intersection metric, which performed the third best (3.33). Thus, while the AUPRC is a good indicator of how well every subset is sorted based on its binarized control value, it can fail to capture how well a metric performs in the high precision and low recall section of the precision–recall curve, as with the propagation intersection metric in this case. To focus on this regime of the precision–recall curve, we find the percentile cutoff value for which the precision drops below 0.95 and use this 0.95 precision percentile cutoff value, which we refer to as the 0.95 precision point, to identify the high precision and low recall FVS subsets.

To illustrate how the intersection metrics perform in identifying successful FVS subsets when compared to other intersection metrics, we plot the percentile cutoffs of each intersection metric against the percentile cutoffs of the CheiRank in Figs. 6(d)–6(f). Successful FVS subsets are represented by filled circles and unsuccessful FVS subsets by open shapes. We use horizontal and vertical dotted lines to mark the 0.95 precision point for CheiRank and the intersection metric, respectively. These two lines split each figure panel into four quadrants, which we can use to evaluate how the FVS subsets identified by the intersection metric differ from those of CheiRank. The top right quadrant corresponds to the FVS subsets identified by both CheiRank and the intersection metric, the bottom right quadrant corresponds to the FVS subsets that CheiRank identifies but are absent from the intersection metric, and the top left quadrant corresponds to the new FVS subsets that the intersection metric identifies that were not identified with CheiRank alone. As shown in Figs. 6(d)–6(f), the number of FVS subsets with high percentile cutoff value and a binarized control value of one (control value $<0.9$⁠) increases as more metrics are added to the intersection (25 for the Modified PRINCE and CheiRank intersection metric to 42 for the propagation intersection metric). The reason for this is that the unsuccessful FVS subsets (open data points) have a larger decrease in percentile cutoff (shift down farther in the plot) than the successful FVS subsets (filled data points), so more FVS subsets are in the top section before the precision drops below 0.95. Thus, as we increase the number of metrics in our intersection metric, high percentile cutoff and high control FVS subsets remain above the 0.95 precision point value (horizontal line) while high-percentile cutoff and low-control FVS subsets shift to a lower percentile cutoff value, which lowers the 0.95 precision point cutoff value and increases the number of FVS subsets above cutoff, creating a higher recall before the precision decreases.

As an example of the improvement brought by intersection metrics, consider the cluster of high CheiRank percentile cutoff and low control value FVS subsets (open diamonds). As more metrics are added to the intersection in each subsequent panel, these FVS subsets’ percentile cutoffs decrease (i.e., they are shifted down on the y-axis) farther than the FVS subsets with a similar percentile cutoff value and successful control value. Most of the FVS subsets marked with open diamonds (8 out of 10) contain the node JAK, which is part of a 2-node negative feedback loop. Negative feedback loops do not contribute to the multistability of the network, so their control has a weak contribution to driving the system into or away from an attractor. Since CheiRank does not consider edge signs, it assigns JAK the 4th highest value among all nodes. JAK ranks 15th for PRINCE, which does consider edge signs, so in Fig. 6(c), we see these open diamonds have much lower percentile cutoff values. As more topological metrics are added to the intersection metric, every FVS subset’s percentile cutoff will decrease or stay the same. In particular, when a FVS subset is bottom ranking for a newly introduced topological metric, the decrease in its percentile cutoff will be larger than that of subsets that are top ranking in the new metric. This also occurs for the other high percentile cutoff FVS subsets that do not include JAK (2 out of 10). These subsets always contain either CTLA4 or TCR, which form an isolated negative cycle, so when the positive cycles metric is included in the intersection, these FVS subsets also greatly decrease their percentile cutoff.

Based on the above results, intersection metrics appear to be more predictive than the individual metrics in the high precision and low recall regime, and we attribute this to the ability of the metrics in the intersection to complement each other’s weakness. We further support these results with Tables S2 and S3 in the supplementary material where we indicate the 0.95 precision point for all four networks (Table S2 in the supplementary material) and the number of FVS subsets above the 0.95 precision point (Table S3 in the supplementary material) for subsets of size 1–3 using all three intersection metrics and CheiRank. Overall, these data indicate that the percentile cutoff decreases and more FVS subsets are identified as more metrics are added to the intersection metric. However, it is possible for the intersection to become too strict if one of the metrics in the intersection sorts the data poorly. This negatively affects the number of subsets identified, so when using the all metric intersection fewer FVS subsets are identified than in the propagation intersection. These results show that the intersection metrics have higher predictive power than any individual topological metric in the high precision and low recall regime.

Knowing that at certain percentile cutoffs the intersection metrics can better identify successful subsets than individual metrics, we used our four networks to determine which percentile cutoff was the most appropriate for picking successful FVS subsets. We aimed to find a percentile cutoff with high precision while maximizing the number of identified FVS subsets. This involves a trade-off because as the precision increases the number of FVS subsets decreases, as shown in Fig. 7. For the rest of the manuscript, we chose the point where the precision first crosses 0.9 from below (dashed line on Fig. 7), which we refer to as the 0.9 precision crossing point. These crossing points were identified for all four networks on FVS subsets of size 1, 2, and 3, for the propagation intersection metric, separately for To Control and Away Control.

Based on the full set of these propagation intersection metric 0.9 precision crossing points for each control measure, two percentile cutoff values were chosen, which we expect will be able to identify successful FVS subsets when applied to other Boolean networks. The maximum of the 0.9 precision crossing points indicates the value that is most likely to achieve a high precision, but it is possible that it does not detect any FVS subsets. To increase the likelihood of identifying a FVS subset, we also evaluated the third quartile percentile cutoff because it is less stringent, but is still expected to have high precision. For To Control, we found that the 93 percentile was the maximum cutoff value and the 83 percentile was the third quartile. For Away Control, these two values were the 90 and 72.75 percentile cutoffs, respectively. The Away Control values are lower because it is easier for a subset to achieve high Away Control than To Control, so a less stringent percentile cutoff performs better.

To verify that our method does identify a sufficient number of successful FVS subsets with high precision, we reapplied the percentile cutoffs identified by our method to the four networks used to derive these percentile cutoffs. For each control measure, both cutoffs were applied on the propagation intersection metric (Table III) and the all metric intersection metric (Table S5 in the supplementary material). In Table III, the number of identified FVS subsets of size 2 with percentile cutoff higher than these values was recorded as well as the precision in terms of being successful FVS subsets; data for FVS subsets of size 1 or 3 are presented in Table S5. In Table III, the FVS subsets identified using the maximum percentile cutoff value have a precision of 100% in all cases except for the To Control case of the FA/BRCA variant network, and have a precision higher than 85%/95% in the To Control/Away Control cases, respectively (Table S5 in the supplementary material). As expected, the third quartile identifies more subsets and has a lower precision than the maximum, but overall using the third quartile percentiles still identifies FVS subsets that achieve a precision above 60% for To Control and a precision above 75% for Away Control. As a point of reference, we also show the size and precision of the identified FVS subsets for both a balanced and unbalanced logistic regression (see the  Appendix). When compared to the logistic regressions, finding the FVS subsets above these percentile cutoffs is more stringent but also typically more precise because the logistic regression does not focus on the high precision and low recall regime.

While using these percentile cutoffs has a higher precision than the logistic regressions in Table III, this is not true in the To Control case of this variant of the FA/BRCA network. This outlier results from the node ICL being crucial for driving the system to the correct attractor, but of all the single nodes whose interventions are fully informative, ICL has the lowest topological values. Thus, if a two- or three-node subset has a high percentile cutoff value, it likely would not contain ICL and would not drive the system to the correct attractor, so the subsets with percentile cutoffs above our percentile cutoff values would not be successful. Logistic regression is still able to identify some successful FVS subsets because the logistic fit leverages the observation that lower percentile cutoffs tend to have higher control values for this model. This observation cannot be incorporated in our percentile cutoff method, which assumes that FVS subsets with high percentile cutoffs are correlated with high control values. This outlier demonstrates that there are complexities of the dynamics of a network that may not be easily captured by topological metrics. Nevertheless, because applying the percentile cutoff values to identify FVS subsets on our networks was able to identify a sufficient number of FVS subsets that have a precision higher than 85% when using the maximum percentile cutoff, we hypothesize that the percentile cutoff values identified by our method would have a similar performance when applied to other networks.

To test this hypothesis, we applied the percentile cutoffs identified by our method to identify FVS subsets that are expected to have high control values in a second array of Boolean models of biological networks. These networks were: a geroconversion network,³⁷ two more variants of the FA/BRCA network, and two variants of a MAPK network³⁸ (see Table S4 in the supplementary material). We identified and found the precision of the subsets using our maximum and third quartile percentile cutoff values (see Table S6 in the supplementary material). Table IV shows the results of the propagation intersection metric for two-node FVS subsets. Many of the top-ranking FVS subsets had a binarized control value of 1 on their respective networks. However, this second array of network models had fewer nodes than our original array, so the most restrictive choice often did not result in any FVS subsets being identified (see Table S6 in the supplementary material). So when applying this method to another network, it may be useful to start with the third quartile percentile values (83 for To Control and 72.75 for Away Control) despite their potentially lower precision. The results of Table S6 in the supplementary material show that the sets of FVS subsets with percentile cutoffs above the third quartile percentile cutoff values are of sufficient size and high precision on other networks too. This shows that our percentile cutoff method is able to identify FVS subsets in the high precision and low recall regime using only the structure of the network.

To confirm that this method identifies successful FVS subsets with high precision, we compared our results to random samples of node subsets. Two different random samples were generated: a random sample of every possible subset and a random sample of FVS subsets. To approximate the precision value of all subsets, we use bootstrapping in each random sample (see the  Appendix).Twenty-five randomly selected single node subsets were chosen, or every individual FVS node was chosen if there were less than 25 total options. For 2 and 3 node subsets, 100 random subsets were picked when looking at all possible subsets and 50 random FVS subsets were picked when only looking at the FVS subsets. After a random sample was chosen, each subset’s $2L$ possible interventions were simulated to determine their To Control and Away Control values. We use the random subsets’ binarized control values to determine their precision and compare it to the FVS subsets identified using the intersection metrics.

The distribution of precision values from bootstrapping is plotted as a box-and-whisker plot together with the precisions obtained when using the maximum percentile cutoff and third quartile percentile cutoff determined in the previous analysis of the FVS subsets, shown in Fig. 8 for the original array of networks and Fig. 9 for the new array of networks. In our networks, the precisions of the FVS subsets identified using our method significantly outperformed the random samples’ precisions. This is true for each intersection metric, and even works partially well for CheiRank. The only exceptions were the To Control cases of both the first variant of the FA/BRCA network and geroconversion network cases and when the intersection did not identify any subsets (indicated by an open shape).

In the following, we analyze the FVS subsets identified using our percentile cutoff method and discuss their biological significance for two of the analyzed networks. We also provide a similar analysis of the third variant of the FA/BRCA network in the supplementary material Text I.

T-cell large granular lymphocytic (T-LGL) leukemia is a disease that is caused by activated cytotoxic T-cells surviving and possibly proliferating instead of undergoing activation induced cell death. Zhang et al. constructed a signal transduction network of activation induced cell death and its deregulation in T-LGL leukemia and modeled it with a Boolean model of this process.³² Here, we use a modified version of the model.²⁶ This T-LGL network contains 60 nodes and 141 edges; we identified several near-minimal FVSs, all of size 12. The model has three attractors, two survival attractors with overlapping node states that include the OFF state of the node Apoptosis, and one apoptosis attractor which includes the ON state of the node Apoptosis. Using our method, we identified multiple FVS subsets; in this follow-up analysis we aim to evaluate which of the $2L$ interventions for each FVS subset drives the system out of the survival attractors and into the apoptosis attractor. As a state with Apoptosis=1 reflects the biological commitment to apoptosis even if it is not identical to the system’s apoptosis attractor, we can analyze a FVS subset’s effectiveness either by the state of Apoptosis after each simulation or by the To Control and Away Control values.

The top five single node FVS subsets predicted by the propagation intersection metric were IL2RB, NFKB, RAS, S1P, and TBET. The intervention expected to be successful in driving the system into the apoptosis attractor is to set the target node into its state corresponding to the apoptosis attractor. This means the OFF state for S1P. Indeed, knockout of S1P drives the system away from the survival state and into the apoptosis attractor a majority of the time (Table V). However, the situation is less clear for the other four nodes in this list because they are ON in all attractors; these nodes are partially informative when used alone and become fully informative when used in combination with S1P. After simulating both the knockout and constitutive expression of each of these four nodes in the Boolean model, we found that knockouts of NFKB, IL2RB, or RAS increases the basin of the survival state to $<90$%. The resulting attractor is not close to the apoptosis attractor but it does have the Apoptosis node in the ON state (Table V). The effectiveness of S1P, NFKB, and RAS knockout in a practical setting can be corroborated with experimental evidence.^39–41 TBET has a smaller Apoptosis basin of attraction than the other four identified interventions, however, it did drive the system to new attractors that are reasonably different from the original survival attractor and have Apoptosis oscillating. We also simulated the fully informative single node intervention with the next highest percentile cutoff, FLIP, and found that

knocking it out does not meaningfully change the basins compared to the unperturbed (wild-type) system. This indicates that fully informative FVS subsets might not always be better options than partially informative ones in terms of their ability to control the network.

We also identified the top 5 two and three node FVS subsets of the propagation intersection metric. We simulated all four (for two nodes) or eight (for three nodes) combinations of knockouts and constitutive expressions for each subset. We identified the intervention that best drove the system to the apoptosis state and recorded this intervention in Table V. The two node interventions drive the system to an apoptosis state more consistently than the single node interventions, but many of the interventions were not driving the system to the wild-type apoptosis attractor.

Furthermore, the best performing combined intervention is not always a combination of the best performing single interventions. For example, the intervention set {NFKB = 1, S1P = 0} is more effective than the set {NFKB = 0, S1P = 0}. Once NFKB = 1 and S1P = 0 are locked in, the system will almost always be driven to the wild-type apoptosis attractor. Upon further analysis of the Boolean model, we discovered the mechanisms that determine which state of NFKB will be more effective. Depending on the state of NFKB, two different subgraphs determine the attractor [Fig. 10(a)]. Certain state configurations of these subgraphs are stable motifs, meaning that these state configurations can be maintained regardless of the rest of the network. Stable motifs’ role in driving a Boolean system to its attractors has been studied in detail.^26,42 In the wild-type system, NFKB naturally turns ON, so the S1P subgraph determines the attractor [Fig. 10(b)]. Turning S1P OFF drives the system into the Apoptosis state. Conversely, when NFKB is fixed OFF, the TBET subgraph determines the state of the system [Fig. 10(c)]. Fixing TBET, IL2RB, or JAK ON drives the system to the Apoptosis state. However, fixing S1P OFF cannot return the system to the original wild-type apoptosis attractor and instead drives to new survival and apoptosis attractors. This example illustrates why the state of the nodes in a FVS subset with the highest ability to drive the system towards or away from an attractor of interest does not always correspond to the one matching the nodes states of the target attractor. In cases with nodes whose state is the same in every attractor, such as NFKB, either fixing the node in the state of the attractors or in the state opposite the attractors could be better at achieving high control values depending on the model.

In the three node case, setting every node in the FVS subset to its state in the wild-type apoptosis attractor is the most successful intervention for reaching an attractor with Apoptosis ON. The top ranking FVS subsets almost all achieve a To Control value of 1 and all achieve an Away Control value of 1, and every intervention contains the combination {NFKB = 1, S1P = 0}. Thus, when any node in the subset is distinct between attractors, such as S1P, the ambiguity in choosing the states for other nodes in the subset that are not distinct between attractors, such as NFKB, is no longer an issue.

We identified seven topological metrics that had the potential to sort FVS subsets by their ability to drive the dynamic trajectory of a biological network. Verifying the effectiveness of each intervention in multiple Boolean models of biological networks, we determined the predictive power of each topological metric. We found that the centrality measures performed adequately. This moderate performance may be because centrality measures only glean local or path-based information from the network and ignore the cycles, which are crucial for understanding the multistability in the network. To our surprise, we found that the cycle-based metrics performed poorly. This indicates that only using the cycle structure of the network does not successfully differentiate between FVS subsets. This poor predictive power may be because the FVS already implements knowledge of the network’s cycle structure and reusing the cycle structure in the cycle-based metrics results in diminishing returns. We demonstrated that the three propagation metrics, which utilize both path information and cycle information, are the best at ranking FVS subsets based on their ability to control the network.

We further demonstrated that taking an intersection of the top ranking subsets of these three propagation metrics is even more successful. Using these intersections, we focused on the high precision and low recall regime because we expect FVS subsets with high intersection metric values to be sufficient but not necessary for having a high control value. While an intersection of all seven metrics captures more structural properties of the network than the propagation metrics, the all metric intersection metric is often too strict because it includes topological metrics with little predictive ability. Based on the analysis of the studied Boolean models, we find that the propagation intersection metric is less stringent, more precise, and identifies more intervention FVS subsets than the all metric intersection metric in the high precision and low recall regime.

To explore how well the propagation intersection metric identifies FVS subsets with high control value, we simulate the top ranking FVS subsets identified by the propagation intersection metric. We found that the top ranking FVS subsets have high control values for two arrays of Boolean networks: the array of networks we used to study our topological and intersection metrics and a second array of networks. We also verified that these identified FVS subsets perform better than randomly chosen node subsets on the same two arrays of networks. Our in-depth analysis of the interventions predicted for two specific models led to general insights. Looking at the T-LGL network, we highlight an ambiguity when trying to identify the intervention state with the highest control value for nodes that have the same state in every attractor. We find that this ambiguity is eliminated when a node whose state is distinct between attractors is added to the intervention. Looking at the 3rd variant of the FA/BRCA network, we confirm that if a FVS subset performs well, then supersets of this subset tend to achieve control values that are equal or greater than the original FVS subset. Using this analysis, we support our hypothesis that the top ranking FVS subsets identified by the propagation intersection metric are successful at driving networks to and away from their attractors.

In summary, our approach successfully identifies few-node combinatorial interventions that can drive a network into a desired attractor and away from undesired attractors. For many of the networks in our study we were able to identify two or three-node interventions that controlled the network. The largest benefit of this approach is that it is based solely on the topology of the system, so it can be used for systems whose dynamics are poorly characterized or to identify interventions that are robust to the dynamic details of the system. Our method has a wide breadth of applications among biological as well as non-biological networks. Being able to identify crucial, attractor-driving nodes in dynamical systems has practical implications such as identifying targets of restorative or preventative interventions; however, it is important to note that this approach only identifies which variables would need to be controlled and assumes we can override the state of these variables. Overall, our study contributes to the body of work documenting how certain dynamical behaviors of a network model, and the methods to elicit them, can be fully determined given only topological knowledge of the network.

See the supplementary material for expanded datasets and descriptions of the analyzed networks. Table S1 provides the expanded AUPRC ranks from Table II including the average rank of the AUPRC values for every intersection metric and topological metric. Tables S2 and S3 provide the data used for the analysis in Sec. III B. Table S4 includes the key properties of the two arrays of networks we analyzed. Tables S5 and S6 are excel files that include the full data for the all metric intersection metric and propagation intersection metric on the first array of networks (Table S5) and the second array of networks (Table S6). Text I includes an analysis of the third variant of the FA/BRCA network, which also includes Table S7.

This work was supported by NSF Grant Nos. IIS 1814405 and MCB 1715826.

The authors have no conflicts to disclose.

J.G.T.Z. conceptualized the idea and J.G.T.Z. and E.N. developed the methodology with input from R.A. R.A. acquired funding and provided project administration. E.N. programed the software with help from J.G.T.Z. E.N. investigated the problem through conducting simulations and provided formal analysis of the resulting data under the supervision of both R.A. and J.G.T.Z. E.N. visualized the data and wrote the original draft, and J.G.T.Z. and R.A. reviewed and edited the final work. J.G.T.Z. and R.A. contributed equally to this work as senior authors.

The data that support the findings of this study are available from the corresponding author upon reasonable request. The code used to generate these findings is available in GitHub at https://github.com/EliNewby/FVSSubsets, Ref. 43.

First, we introduce key notations. The number of nodes in the network is denoted by $N$⁠, and the number of attractors is denoted by $M$⁠. Individual attractors are denoted $Ai$⁠, where $i$ goes from $1$ to $M$⁠. Individual nodes are denoted by $n$ in the network. The state of node $n$ in attractor $i$ is represented as $sni$⁠. A FVS subset of size $L$ is denoted by $S⊆FVS={n1,n2,…nL}$⁠. An intervention of a FVS subset is denoted by $P={p1(n1),…,pL(nL)}$⁠. Here, $pj(nj)$ denotes an individual intervention on node $nj$ (e.g., $P={p1(S1P),p2(NFKB)}={0,1}$ means S1P is being driven to its OFF state and NFKB is driven to its ON state).

A Boolean model is a discrete dynamic model that characterizes each node of a network with a state variable that can take one of two values. These values are 0 and 1, interpreted as OFF (inactive) and ON (active), respectively. Each node is characterized with an update function usually expressed by logical operations. This update function takes as input the state values of other nodes in the network, and its output determines the node’s state at the next time step. We use a stochastic asynchronous update scheme in which at every time step a randomly selected node’s state is updated by evaluating its update function. An asynchronous update scheme updates nodes individually and in a stochastic manner, as opposed to a synchronous update.

Attractors are sections of the state space that the system can enter but cannot exit. They can be either simple attractors (e.g., fixed points) or complex attractors (e.g., stable oscillations). Along with each attractor, it is valuable to identify the attractor’s basin of attraction—the region of state space from which every initial state will lead to that specific attractor.

Given that minimal and near-minimal FVS are non-unique,^13,14,36 we identify multiple near-minimal FVS, obtaining a superset of every node that could be in any of the near-minimal FVS. When finding multi-node subsets, we only analyze subsets that are part of the same near-minimal FVS. In other words, we ignore subsets that contain nodes that appear in separate near-minimal FVSs.

In the following, we describe each topological metric. Of note, the computational complexity of the propagation metrics is higher than that of the centrality metrics, but their slight computational disadvantage in our models is much less than the informational advantage we observe when using the propagation metrics.

The sum of the out-degrees is a local measure that measures the number of nodes the subset is directly connected to

The distance from node $ni$ to node $nj$ (⁠$lij$⁠) is a measure of how close these nodes are to each other, but to account for unreachable nodes, the inverse distance $1/(1+lij)$ is used. The average of these inverse distances from a node to every other node in the network defines the distance metric. When working with a subset of more than one node (⁠$L>1$⁠), the distance $lij$ used in the equation is the distance to the closest node in the subset, so it is defined as

Edge sign is not incorporated into the distance metric, thus it cannot distinguish between positive and negative paths.

The PRINCE Propagation metric reflects how much information from a node propagates through the system. It is formulated by applying a constant perturbation on the specified node subset and then allowing the perturbation’s information to be distributed to the node’s neighbors according to a discrete time dynamic process. This initial perturbation vector $Y→$ is a discrete value vector of size $N$⁠, where the value corresponding to each node is either 1 (fixing the node ON), 0 (no effect on the node), or $−$1 (fixing the node OFF).

The value propagated through an edge depends on the edge’s sign, i.e., a negative edge propagates a negative value. This dynamic process attributes every node a value between $−$1 and 1. The value is the node’s information score, which represents how much information from the perturbation is received. An absolute value of 1 means we have full information about the node, while a value of 0 means that we do not know anything about the given node.

The information scores of every node in the system are based on a normalization scheme and a propagation drop-off value $α$⁠. The normalization scheme mimics mass flow. The out-going edges split the information evenly between every successor. The amount of intake is also distributed evenly to every predecessor. The propagation dropoff value causes the information from an input node to decrease the farther it travels. The propagation is calculated based on the algorithm

where $π(t)$ is the propagation vector at time t, $W′=D1−1/2WD2−1/2$ is the normalized adjacency matrix, in which $D1$ is a diagonal matrix where element $D1(i,i)$ is the out-degree of node $i$⁠, $D2$ is a diagonal matrix where element $D2(i,i)$ is the in-degree of node $i$⁠, and $W$ is the signed adjacency matrix where element $Wij$ is 1 if there is a positive edge from node $i$ to node $j$⁠, $−$1 if the edge is negative, and 0 if there is no edge between nodes $i$ and $j$⁠. The propagation dropoff value is denoted by $α$⁠, which we set to $α=0.9$ for our calculations following the value used by Santolini et al.³⁰ 

The long term behavior of the system can be characterized by its steady state $π∗$⁠, which is such that $π∗=αW′Tπ∗+(1−α)Y→$⁠. To calculate $π∗$⁠, we follow Ref. 30, which shows $π∗=(1−α)Y→(I−αW′T)−1$⁠. To measure how much information we have on every node, the absolute value of $π∗$ is taken. The PRINCE Propagation value for a specific subset is the average of these absolute values of the information scores taken over every node in the network (⁠$π∗¯)$⁠.

Because the absolute values of the information scores are taken, an intervention of 1 and $−$1 return the same value in a one-node intervention, so only the inputs of 1 are tested. When considering a multi-node perturbation, some further modifications need to be used. Instead of a constant input of 1 to every node in the subset, all combinations of inputs of 1 and $−$1 are taken, and the maximum of these results is used as the PRINCE value for this subset.

We propose a variant of the PRINCE Propagation algorithm that only normalizes the propagation by the in-degree of each node. This better treats the initial perturbation as information instead of mass, so it can fully spread to all successors instead of being split among them. To keep the system convergent, the adjacency matrix is changed to $W′=WD2−1$⁠. Other than adjusting the adjacency matrix, the values for this metric are obtained the same way as in the original PRINCE Propagation. It still accounts for edge signs by allowing for the propagation value to be negative, and it also still treats multi-node interventions by considering every combination of input values.

CheiRank is an extension of PageRank, which consists of a modified random walk over the nodes of a network. CheiRank follows the same algorithm as PageRank but traverses the network’s edges in reverse. The random walk is augmented with a probability $α$⁠, where there is a probability of $1−α$ that the walker jumps to a random node instead of following an edge. The probability vector $π(t)$ of the random walk follows

where $A$ is the adjacency matrix of the network with each row normalized so it sums to 1, so there is an equal probability to travel to any neighbor nodes, and $E$ is an all 1 matrix.

In the long-time limit, the system will reach a unique steady state $π∗$⁠, that satisfies $π∗=[αAT+(1−α)E/Nn]π∗$⁠, where each value in $π∗$ is the CheiRank of that node, and it is the probability that the random walker is found in this specific node. Because CheiRank denotes a probability, when the subset has multiple nodes in it, the probabilities are combined multiplicatively in the following way: $CheiRank(n0,…nL)=1−(1−CheiRank(n0))∗⋯∗(1−CheiRank(nL)).$

The cycle-based metrics estimate the multistability reintroduced into the system after the FVS is removed and a subset of it is reintroduced. We used two topological measures as an estimate of multi-stability: the number of positive cycles (cycles with an even number of negative edges) and the size of the strongly connected component. Both measures return 0 for subsets containing nodes that are not in the FVS.

The number of positive cycles of a FVS subset $S$ is defined as the number of positive cycles left in the network after the FVS has been removed and the subset has been reintroduced,

The SCC of a FVS subset $S$ is defined as the size of the largest SCC after the FVS has been removed and the subset has been reintroduced,

Unlike the cycles, the signs of the edges are not incorporated in this metric.

To determine the attractors of the wild-type (unperturbed) system, we use the trapspace method in the BioLQM toolkit.⁴⁴ Technically, the objects identified by this method are the minimal trap spaces (subsets of the state space specified by fixing the value of a set of nodes which if entered cannot be exited), not the attractors, but there is often a one-to-one correspondence between them, and trap spaces can be identified efficiently.⁴⁵ The attractors include fixed points (steady states) in which the state of all the nodes is fixed and complex attractors, in which the system keeps revisiting a finite set of states. We calculate each attractor’s basin of attraction by performing 1000 simulations starting from random initial conditions and using a general asynchronous update. The simulation continues until an attractor is reached, and we count how many of the simulations reach each of the system’s attractors to determine the basins of attraction.

For each intervention, we perform 100 simulations from random initial conditions using a general asynchronous update. After each simulated intervention, we measure how close the final state of the system, $Sim(P→)$⁠, is to the attractors of the wild-type system, $Ai$⁠. We measure the closeness to each attractor through a normalized Hamming distance. The Hamming distance $h(S1,S2)$ measures the fraction of nodes with different values between two different states $S1$ and $S2$⁠. The normalized Hamming distance between the final state and every wild-type attractor is used to generate an attractor control value $C$ for each attractor defined as

The attractor control value is a linearly decreasing function with a cutoff $F$⁠, which indicates when the state is too far away from our attractor. The $F$ parameter is defined to be halfway between the two closest attractors in the network, $F=mini,jh(Ai,Aj)/2$⁠. This ensures that each final state in a simulation is attributed to one attractor and the other attractors get a value of 0 for that specific simulation.

For each simulation, we get an attractor control value for every attractor, but only one is nonzero. These values are averaged over the 100 simulations with different, random initial conditions. The average for each attractor is considered the basin of attraction for that attractor in the perturbed system, $B(Ai|P)=C(h(Ai,Sim(P));F)¯$⁠. The intervention basin vector is defined as $B(P)→=(B(A1|P),B(A2|P),…,B(AM|P))$⁠. As a point of comparison, we consolidated the wild-type basin values into the wild-type basin vector, $WTB→=(WTB(A1),WTB(A2),…,WTB(AM))$⁠, where $WTB(Ai)=B(Ai|0)=C(h(Ai,Sim(0));F)¯$ is the wild-type basin of each attractor, which is averaged over 1000 simulations.

To study the effects of an individual intervention, we separated the wild-type system’s attractors into two categories: target attractors and non-target attractors. These two categories are uniquely determined based on each specific intervention. A target attractor is an attractor wherein the intervention state of each node is the same as the node’s state in the attractor, and a non-target attractor is an attractor wherein the intervention state of any node is not in the attractor’s state.

We describe this using indicator functions. The indicator function $I(snji,pj(nj))$ compares the state of node $nj$ in a specific attractor, $snji$⁠, to the state the node is fixed into through the intervention, $pj(nj)$⁠. If they are the same, the function returns 1 and if they are different the function returns 0, so it is defined as

To extend this to multiple nodes, the indicator functions for every node in the intervention are multiplied together. These products are merged into a single identification vector of size $M$⁠, $I(P)→={∏j=1L[I(snj1,pj(nj))],…,∏j=1L[I(snjM,pj(nj))]}$⁠. The attractors for which the corresponding elements of the identification vector are 1 are the target attractors and the attractors for which the elements are 0 are the non-target attractors.

After labeling each attractor as a target or non-target, we filter out interventions that are considered uninformative. Interventions where every attractor is a target attractor are uninformative because we cannot attribute changes in an attractor’s basin to the intervention. Similarly, the To Control metric for an intervention is uninformative if every attractor is considered a non-target. Because we are not driving to any attractor, we cannot measure how often we reach a target attractor, but we can still measure if the intervention drives the system away from the attractors of the system, so we denote these interventions as partially informative. We denote an intervention as fully informative if the intervention has attractors that are targets and attractors that are non-targets.

For every (fully or partially) informative intervention, we sum the wild-type basins of every target attractor into one wild-type target basin $WTTB(P)=I(P)→∙WTB→$ and sum the basins for every Non-target attractor into one wild-type non-target basin $WTNTB(P)=1−WTTB(P)=(1→−I(P)→)∙WTB→$⁠. Similarly, we combine the intervention attractor’s basins into an intervention target basin $TB(P)=I(P)→∙B(P)→$ and an intervention non-target basin $NTB(P)=1−TB(P)=(1→−I(P)→)∙B(P)→$⁠.

We then find the percent difference between these basins to understand how the intervention affects the system. The difference between the intervention target basin and wild-type target basin is the intervention-specific To Control, and the difference between the intervention non-target basin and wild-type non-target basin is the intervention-specific Away Control. Correspondingly, the To Control and Away Control values range between $−$1 and 1,

Once all $2L$ possible interventions are tested or filtered out for a subset, we aggregate the values into a singular aggregate To Control and Away Control values for the entire subset by taking the maximum of all available values. Taking a maximum of all the $2L$ interventions ensures that the aggregate To Control and Away Control values of the entire FVS is 1, indicating full control over the network

Logistic regressions are a type of binary classifiers. To convert our continuous FVS subset control values to a binarized value of a successful vs unsuccessful FVS subset, we impose a threshold. We used a threshold of 0.9 to binarize the control values and found there were no major changes in our results using thresholds in the range 0.5–0.95. The logistic regression fits a logistic function to the binarized data and results in a curve in which the x-values of the fit are the values of the topological metrics and the y-values indicate the probability of that topological value resulting in successful control (i.e., control above 0.9). Therefore, the y-value for the logistic function is between 0 and 1, and in an ideal situation, it starts at 0 for low x-values, increases to 1 as the x-value increases, reflecting how we expect topological values to be positively correlated with control.

This regression can be modified by applying weights to the negative and positive data points. By default, we create an unbalanced regression which weights every data point equally. We also test a balanced regression where data points are weighted by the number of points in their dataset, i.e., positive data points are weighted by the relative number of positive data points and negative data points are weighted by the relative number of negative data points. This means the positive dataset and negative dataset are weighted to correct for differences in their relative sizes.

We use the precision–recall curve and the area under the precision–recall curve (AUPRC) to evaluate the strength of the fit. We choose the AUPRC over the more common receiver operating characteristic (ROC) curve because there is often a large imbalance in the binarized control values, which limits the usefulness of the ROC curve for our analysis.⁴⁶ To generate a precision–recall curve, the entire range of a topological value is scanned. As each value is scanned over, the value of the logistic regression function is used to sort the data into four categories: true positives $TP$⁠, false positives $FP$⁠, false negatives $FN$⁠, and true negatives $TN$⁠. Then, based on these values, a precision $TP/(TP+FP)$ and recall $TP/(TP+FN)$ value can be calculated. As the precision and recall are calculated, they are plotted on a curve, so when the entire x-axis is scanned, the precision–recall curve is formed. The AUPRC is calculated from the precision–recall curve and consolidates the curve into a single value that can indicate how well the topological metric sorts the data.

The AUPRC is a positive value less than or equal to 1, where 1 indicates that the topological metric perfectly sorts between positive and negative binarized control values. When the AUPRC value is equal to the fraction of positive data points, it indicates that the topological metric does not sort the FVS subsets better than random. As performing better than random is not sufficient to confirm that the topological metric can predict which FVS subsets are successful, we introduced an AURPC predictive threshold, which is halfway between the random expectation and the maximum score of 1 [$1−(1−Positive Fraction)/2$].

For each topological metric, the FVS subsets are ordered in increasing order of the value of the metric. Then, this ordered list is converted into percentile values $f(S,m)$ for each FVS subset $S$⁠. For a set of topological metrics of interest, the sets of FVS subsets above a fixed threshold percentile $pc$ are intersected. We refer to metrics that intersect the top ranking FVS subsets of a set of topological metrics as an intersection metric, where each intersection metric is denoted by their set of topological metrics

We investigated three different intersection metrics: an intersection of all seven metrics, an intersection of only the propagation metrics, and an intersection of the CheiRank and Modified PRINCE metrics.

In an intersection metric, we assign to each FVS subset a percentile cutoff value $IP(S,Metrics)$ that is equal to the minimal percentile cutoff value of the FVS subset among the metrics being intersected. This is equivalent to choosing the maximal $pc$ value for which $S$ is found in $IM(pc|TopologicalMetrics)$ given the chosen metrics

The methods were implemented in python using various libraries and modules. The code is available at https://github.com/EliNewby/FVSSubsets. This code consists of two python files: FindBestSubsets.py and RunSimulations.py. FindBestSubsets.py finds the topological metric values and intersection metric values for the FVS subsets of a network topology and uses these values to rank the FVS subsets. RunSimulations takes a list of node sets and a Boolean model for a network and calculates the To Control and Away Control values for each node set. We also include a Jupyter notebook that includes an example showing the output of these two functions on the T-LGL network and that also reproduces the plots included in the figures of this paper.

To identify the near-minimal FVS and its subsets, a python code developed by Gang Yang was used,¹⁴ which utilizes the simulated annealing algorithm presented by Galinier et al.³⁶ This code is available at https://github.com/jgtz/FVS_python3. NetworkX was used to analyze the structure of the networks and to calculate the values of the topological metrics. NetworkX implemented some of the topological metrics we used. Custom Python code and NetworkX functions were used for the ones that were not implemented.

The simulations used to calculate control values were implemented using the bioLQM toolkit developed by the CoLoMoTo Consortium.⁴⁴ Using bioLQM, we implement our Boolean models of the networks, find their attractors, and simulate the system’s trajectories using the random asynchronous update mode.

The logistic regressions were created and analyzed using the scikit-learn module. From the linear_model submodule, the LogisticRegression function is used to generate a fit of the data. This function is implemented using the liblinear solver and a regularization strength of 100. The scikit-learn module was also used to obtain the precision–recall curve (metrics.precision_recall_curve) and to calculate the AUPRC (metrics.average_precision_score).

Bootstrapping was used to approximate the precision of the set of all node subsets from random samples of these subsets. Bootstrapping is a resampling method that generates a distribution from a sample set of subsets to make inferences about the value of an observable (precision in our case) in the set of all node subsets. The sample sets were randomly chosen samples of subsets. We tested two separate sets: the set of every node subset and the set of FVS subsets. For 1 node subsets, the random samples chosen are of size 25 or include all 1 node subsets if there are less than 25 total node subsets. For 2 and 3 node subsets, 100 subsets are chosen for the random sample of all node subsets, and 50 subsets are chosen for the random sample of FVS subsets. After picking a random sample, every subset in the sample was simulated to determine its To Control and Away Control value. For each random sample, we resample with replacement a resampled set of the same size. We resample 1000 times and calculate the precision of each resampled set. From this distribution of these 1000 precision values, we can approximate the precision value of the entire set of node subsets or FVS subsets depending on the original set that the random sample is generated from.