Robust Stabilizing Control of Perturbed Biological Networks via Coordinate Transformation and Algebraic Analysis

This article investigates robust stabilizing control of biological systems modeled by Boolean networks (BNs). A population of BNs is considered where a majority of BNs have the same BN dynamics, but some BNs are inflicted by mutations damaging particular nodes, leading to perturbed dynamics that prohibit global stabilization to the desired attractor. The proposed control strategy consists of two steps. First, the nominal BN is transformed and curtailed into a sub-BN via a simple coordinate transformation and network reduction associated with the desired attractor. The feedback vertex set (FVS) control is then applied to the reduced BN to determine the control inputs for the nominal BN. Next, the control inputs derived in the first step and mutated nodes are applied to the nominal BN so as to identify residual dynamics of perturbed BNs, and additional control inputs are selected according to the canalization effect of each node. The overall control inputs are applied to the BN population, so that the nominal BN converges to the desired attractor and perturbed BNs to their own attractors that are the closest possible to the desired attractor. The performance of the proposed robust control scheme is validated through numerical experiments on random BNs and a complex biological network.

biology [1]- [3].Though dynamics of BNs may seem too simplified compared with complicated aspects of biological systems, they can still represent the major characteristics of biological mechanisms, possibly after streamlining gene expressions via clustering [4] and/or reduction [5] of network dynamics.Indeed, Boolean variables having only 1 (on) and 0 (off) and their transitions described as logical rules on other variables epitomize underlying dynamics of biomolecular systems-activation and inhibition of gene expressions and their evolution governed by the regulatory nodes [6].
The state of a deterministic BN is a collection of the values of all the constituent nodes.Starting from an initial state, BNs always converge to equilibrium states termed attractors.An attractor is either a fixed point or cycle.Attractors represent key phenotypes of long-term behaviors possessed by the corresponding biological system [7], [8].Hence, stabilizing control is one of the most important research topics on BNs, wherein a control law is applied to a set of selected nodes so that the BN can be stabilized to a desired attractor.
In this article, we present a robust stabilizing control scheme for perturbed BNs where a set of nodes are selected as constant control inputs.We consider a cell population in which a majority of members have the same BN dynamics, whereas some cells are inflicted by mutations that perturb specific nodes, resulting in the change of BN dynamics.In particular, we focus on the kind of mutations by which a number of nodes are forced to have the opposite value of the desired attractor.The dynamics of the resultant BNs, termed perturbed BNs, becomes different from that of the nominal BN, and thereby, global stabilization of the perturbed BNs toward the desired attractor gets never achievable.Worse yet, it is often the case that perturbed BNs evolve into undesirable attractors, yielding fatal outcomes.For instance, instead of apoptosis, cancer cells can degenerate into proliferation by mutations on key tumor suppressor genes [9], [10].The best way to alleviate the adverse effect of mutations would be to apply additional control inputs so as to drive perturbed BNs toward some other attractors that have similar phenotypical characteristics to the desired attractor.Furthermore, as the number of control inputs is bounded by practical constraint, the stabilizing control policy must be developed in favor of suppressing not only the number of original control inputs for the nominal BN, but also that of additional control inputs targeting perturbed BNs.

A. Related Work
For the past decade, a number of efforts have been made to develop an elegant mathematical framework of analyzing and controlling BNs based on matrix operations called semi-tensor product (STP) [11]- [13].In view of the previous studies of STP, the present topic is related to simultaneous stabilization of a set of Boolean control networks (BCNs) and set stabilization, namely, a BCN is stabilized with respect to a subset of attractors.In [14], STP-based state and output feedback control laws were proposed for achieving simultaneous stabilization of multiple BCNs.In [15], ensemble controllability of BCNs was introduced whereby a family of BCNs is steered between two states via a free control sequence.In [16], on the other hand, set stabilization of BCNs was addressed in which a timeoptimal controller is presented to make a BCN converge to an invariant subset of states.In [17], the set stabilization problem for probabilistic BNs with randomly chosen Boolean functions was further investigated.
The proposed methodology is also related to pinning control of BNs in that only a fraction of nodes are selected as control inputs.In [18], pinned nodes are analytically selected in the framework of STP so as to achieve global stabilization of BNs.In [19], a pinning control method was presented for the disturbance decoupling problem of BNs.In [20], a single-input pinning control method was proposed to steer BNs toward a desired state via a free control sequence.In [21], finally, model reduction and pinning control were combined to make a BN converge toward a desired attractor.
Besides STP-based approaches, a number of previous studies exist on analysis and control pertaining to perturbations or mutations in biological systems with logical models.In rough terms, they can be classified according to the types of perturbation.In [22]- [24], the effect of function perturbations in BNs was studied, while gene or node related mutations or disturbances were studied in [25]- [27].In addition, edgerelated perturbations in BNs were considered in [28] and [29].

B. Contributions and Paper Organization
The proposed methodology of robust stabilization consists of two steps.In the first step, we employ the attractorspecific coordinate transformation and reduction technique introduced in [30] as the stabilizing control law for nominal BNs.Given the desired fixed point attractor, the controlled BN is transformed into an equivalent one having the minimum sum-of-product (SOP) expression and is reduced further while maintaining the desired attractor.The feedback vertex set (FVS) control [31], [32] is then applied to the reduced BN to identify the control inputs that can achieve global stabilization to the desired attractor.Since the identified control inputs can also deal with the stabilizing control problem for the given nominal BN, the number of control inputs is curtailed compared with the original FVS control and other eminent stabilizing control methods.
In the second step, we tackle the problem of robust stabilization by first providing the nominal BN with the control inputs derived in the first step and mutated node values to single out the dynamics of perturbed BNs.We then quantify the canalization effect [33] of remaining nodes, namely, investigating how many state variables are fixed in favor of the desired attractor by assigning constant control to each node.The additional control inputs are determined according to the total canalization number over perturbed BNs.To demonstrate the superiority and applicability of the proposed scheme, numerical experiments are conducted in which our stabilizing control is applied to random BNs and a real biological network.
In comparison with the prior work, this article has the following contributions.
1) In contrast with the STP-based approaches [14]- [21], this study aims to achieve robust stabilization against model uncertainties in consideration of real biological problems.Furthermore, though [14]- [21] presented analytical solutions to the stabilization problem via matrix computations, they have a serious drawback of exponential complexity.On the other hand, while the complexity of the proposed control scheme is also NP-complete in theory, its computational burden for tackling the stabilization problem of complex biological networks is manageable, mostly due to the network reduction undertaken in the first step of the method.2) Unlike the previous studies in the framework of STP [14]- [21], this study pursues constant control with neither state/output feedback nor free Boolean control sequences utilized.Constant control bears much practical strength, as it can be realized by gene knockout or pharmacological activation or inhibition of molecules [34].In contrast, applying time-varying or feedback control signals to real biological systems is almost infeasible.3) In spite of various approaches and perspectives on perturbations and mutations [22]- [29], no prior work exists on a control theoretic approach to minimizing adverse effects of gene mutations on the collection of BNs as presented in this article.Our study is worthwhile in pragmatic terms as well since it accomplishes robust stabilization of complex biological networks that are computationally intractable as shown in many previous studies.4) Though we adopt the stabilizing control law from our previous result [30], this work has a significant improvement over it.While [30] ensures global stabilization of a given BN toward a desired attractor, it cannot resolve the predicament that a number of nodes are permanently perturbed as opposed to the desired values.Moreover, the stabilizing control law of [30] is applicable only to a single BN.On the other hand, the proposed scheme can tackle the problem of simultaneously stabilizing a number of BNs with different dynamics as affected by mutations.In this way, each BN can be steered toward its own attractor having similar biological characteristics to the desired one.The rest of this article is organized as follows.In Section II, the stabilizing control of nominal BNs without mutation is addressed.In Section III, the main result on robust stabilization for perturbed BNs is proposed based on the coordinate transformation and algebraic manipulation.In Section IV, we validate the proposed scheme by conducting numerical experiments on a population of large-scale random BNs and a real biological network, the metastasis influence network.Finally, Section V concludes this article.

II. STABILIZING CONTROL OF NOMINAL BNS
The notations and stabilizing control law addressed in this section are adopted from our previous work [30].

A. Attractor-Specific Coordinate Transformation
Notation: N denotes the set of non-negative integers.For a finite set A, |A| ∈ N is the cardinality of A. For Boolean variables a, b ∈ {0, 1}, "a " denotes negation or logical NOT of a, "ab" denotes conjunction or logical AND of a and b, and "a + b" denotes disjunction or logical OR of a and b.
A BN having n state variables or nodes x 1 , . . ., x n ∈ {0, 1} is represented by where t ∈ N is the discrete time variable, x(t) is the state vector such that and f i : {0, 1} n → {0, 1}, i = 1, . . ., n, is a Boolean logic serving as the state transition equation of x i .For later usage, let N := {1, . . ., n}. x i (t) will be often denoted by i ∈ N whenever convenient.This BN is alternatively represented by a Boolean mapping Among fixed point attractors and cyclic ones possessed by F, assume that a fixed point attractor c = (c 1 , . . ., c n ) represents the state corresponding to the most desirable phenotype of the biological system characterized by F. T c , a simple coordinate transformation with respect to c, is introduced as follows: We utilize T c to transform F into an equivalent BN G such that Without loss of generality, each g i (y(t)) is supposed to be reduced to the minimum SOP expression after going through the transformation.As (3) signifies, T c is an attractor-specific transformation, whereby the desired attractor c of F is now represented by one vector 1 n := (1, 1, . . ., 1) in G.The latter feature provides a shortcut to solving global stabilization of the original BN F. Since 1 n is a fixed point of G and g i (y(t)) has the SOP expression, each g i (y(t)) contains at least one product (or AND) term consisting only of state variables with no negation operator attached; otherwise, logic 1 could not be a fixed value of the state transition equation y i (t +1) = g i (y(t)).We call such terms plus products.Those terms having state variables with negation operator are called non-plus products.Denote by P i and P i the set of plus products and non-plus products in g i (y(t)), respectively.P i is regarded as a logical embodiment of the convergence tendency of the original BN F toward the i th node value c i of c.The correlation between P i and c i can be understood intuitively.Clearly, the convergence of G to 1 n (or of F to c) is attainable only if one or more terms of P i turn on to 1. Non-plus products in P i do not contribute to the convergence, because as G approaches 1 n , all terms of P i become 0 owing to the existence of state variables with negation operator.

B. Network Reduction and FVS Control
An interesting property of G is that even though all terms except for one plus product are removed from each state transition equation, the reduced BN still maintains 1 n as its fixed point.Furthermore, the control inputs solving the global stabilization problem of the reduced BN with respect to 1 n serve as a solution to the same problem for G [30,Th. 1].Equivalently, the latter implies that the original BN F can be globally stabilized to c by applying the control inputs derived with respect to the reduced BN.
We now present a network reduction scheme that yields the reduced BN H from G such that where z(t) := (z 1 (t), . . ., z n (t)) and h i (z(t)) are derived as follows.
1) Among P i of each y i , i = 1, . . ., n, take one plus product, termed h i (y(t)), according to the selection algorithm addressed in [30, Algorithm 1]. 2) Set h i (z(t)) to be the state transition function of z i in which z i replaces y i , i = 1, . . ., n.In other words, H is derived from G simply by discarding all terms except one plus product from g i (y(t)) and by substituting y(t) with z(t).The selection of the plus product is made in such a way that it minimizes the number of feedback loops associated with the plus product and the rest of state variables [30, Algorithm 1].
The FVS control [31], [32] is applied to H to induce the set of minimum FVSs where 2 N denotes the power set of N. The set notation is needed to depict the general case where more than one minimum FVS are derived with respect to H .We take an arbitrary element with Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
still maintained as an attractor of H , the latter control globally stabilizes H to 1 n .As addressed earlier, H also achieves global stabilization of F to c.The values of control inputs for F are determined by conducting the inverse transformation is defined as follows: Since H is derived with respect to the reduced BN H , it is highly likely that | H | is much less than | F |, where F denotes the minimum FVS derived with respect to F. The superiority of this stabilization scheme over the original FVS control and other eminent stabilizing control schemes is fully validated in [30].
Example 1: Consider a synthetic BN F with ten nodes (n = 10) whose state transition equations are given by Suppose that among six attractors of F c := (0, 1, 0, 0, 1, 1, 0, 0, 0, 1) with 6.2% basin of attraction represents the most desirable phenotype.Then, the coordinate transformation is applied according to (3) so as to derive the transformed BN G such that One can identify plus and non-plus products of each state variable.For instance Activating the plus product selection algorithm [30], we reduce G to H by removing all terms except one plus product for each y i , which is underlined in (9).The network graphs of G and H are illustrated in the left part of Fig. 1.Finally, FVS control to H yields the solution Since c 2 = 1 and c 9 = 0, the control inputs for the original BN F are {x 2 = 1, x 9 = 0}.♦

A. Problem Statement
In this study, we envisage that the population of biological networks modeled by F is vulnerable to mutations, which perturb the value of specific nodes.Since stability toward the desired state c is our main interest, we examine the kind of mutations that fix some nodes to the complementary values of their positions in c, e.g., fixing the i th node to c i .Since G is equivalent with F, we describe perturbed BNs in the formulation of G. Let d, 1 ≤ d < n, be the number of mutations that occur in the given population of biological networks and denote each perturbed BN by G j , j = 1, . . ., d. G j is supposed to have a mutation at the σ j th node, namely, y σ j is fixed to 0 as a result of mutation.Thus, we can represent each G j as Remark 1: Though the abovementioned formulation that a perturbed BN has only one mutation is maintained throughout this section, it is only for the brevity of notations.The proposed scheme imposes no restriction on the number of mutated nodes in a perturbed BN.As will be shown in Section IV, we conduct the numerical experiment on random and biological BNs by setting that each perturbed BN contains up to two mutated nodes.Note that the dynamics of each perturbed BN may significantly diverge from the nominal BN under the influence of multiple mutations.Thus, from a control perspective, G and G 1 , . . ., G d should be regarded as The main objective is to determine a set of constant control inputs so that while G is stabilized to c, each G j , j = 1, . . ., d, is stabilized to its own attractor as near as possible to c.In this article, we use the Hamming distance between two attractors as the proximity measure.Taking into account the requirement that the convergence of the nominal BN F toward c be always ensured, we propose a stabilizing control policy wherein a proper set of m control inputs H ∈ derived with respect to the reduced BN H is applied by default, and a number of additional control inputs are selected to realize robust stabilization for perturbed BNs.Since the number of control inputs is usually subject to practical bounds, we stipulate that the number of total control inputs does not exceed a prescribed limit.Let us formalize, as follows, the robust stabilization problem under the aforementioned conditions.
so that applying H ∪ G stabilizes every perturbed BN G j , j = 1, . . ., d, to its own attractor as near as possible to c where the proximity is measured in the average Hamming distance.
As the foregoing problem statement attests, robust stabilization proposed in this study does not necessarily imply global stabilization of every perturbed BN.G j may globally converge to its own fixed point attractor, or it may reach one of the multiple attractors including cyclic ones.In such a case, the average Hamming distance between the desired attractor and all the possible attractors is used as the proximity measure.If | | ≥ 2, the use of an arbitrary control input set in ensures global stabilization of the nominal BN G.In the case of robust stabilization, by contrast, the efficiency of the solution in terms of the number of control inputs can vary depending on the selected control input set H ∈ .Hence, the robust stabilization problem considers the selection of not only additional control inputs but also the control input set for the nominal BN (provided that | | ≥ 2).
Assumption 1: We will develop our methodology by assuming that Under the abovementioned assumption, the mutation does not occur to the control inputs H .This is not a burdensome requirement, since for a complex BN with many nodes, usually multiple minimum FVSs are derived.Hence, given a profile of possible mutations, one can select H that does not contain any mutated node.The latter scheme will be demonstrated in our numerical experiment.Even though the control inputs H are applied, the existence of the mutated node y σ j prohibits the global stabilization of G j to 1 n .Furthermore, owing to the canalization effect [35] of the mutated node y σ j , some other nodes may be forced to reach 0 or may not be stabilized to either 0 or 1.Hence, selecting additional control inputs from such nodes would be the most reasonable strategy.To this end, we first introduce a notation S 1 for S := {i 1 , . . ., i k } ⊂ N to indicate that all the nodes in the index set S are assigned logic 1, or In association with S 1 , let G(S 1 ) denote the BN G wherein the nodes of S are evaluated to 1.Note that all the subsequent canalization propagated from S 1 is reflected by G(S 1 ).According to this notation and the previous result [30], G( where mutated node values are marked in bold.Clearly, the control inputs 1 H derived in the prior work [30] cannot solve the present robust stabilization problem, since only a subset of nodes converges to the desired value (1's) in two perturbed BNs-80% in G 1 and 50% in G 2 , resulting in average 65%.The network graphs of G 1 ( 1 H ) and G 2 ( 1 H ) are drawn in the right part of Fig. 1.Note that we eliminate all the incoming and outgoing edges of H , σ j ( j = 1, 2), and the nodes fixed to 1.Those nodes not fixed to 1 and their remaining edges clearly illustrate the adverse effect of the mutation.♦ Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

B. Determining Additional Control Inputs via Canalization
As shown in Example 2, G j ( 1 H ) allows us to identify which nodes in G j deviate from the desired value 1 as the result of mutation.To depict the latter information in detail, we divide the node set where σ j is the mutated node, V j z is the set of nodes whose state transition equations take logic 0 or non-constant Boolean logic, and V j n logic 1.By definition, H ⊂ V j n .We also define as the union of V j z values.Example 3: In Example 2, given σ 1 = 3 and σ 2 = 8, V j z , V j n , and V z are derived as ♦ State variables in V z are upset to 0 or become unstable by the mutations.As mentioned before, a reasonable strategy for robust stabilization is to quantify the favorable effect of setting each state variable in V z to be 1 and to determine the additional control input based on the quantification.In doing so, the following two aspects must be taken into consideration.
1) First, determining additional control inputs is conducted in a recursive sense.Once an additional control input is chosen, we derive the quantification of the favorable effect in the next step after reflecting the latest chosen control input (along with the previously selected control inputs).2) Next, the favorable effect of an additional control input is evaluated over all the perturbed BNs, as our objective is to achieve robust control.If choosing a candidate node makes a perturbed BN converge to a fixed point that is close to 1 n while making another perturbed BN to another fixed point very far away from 1 n , it is unlikely that such a node will be selected as the next control input.To formalize the favorable effect of an additional control input, let us define the canalization number of each state variable in V z , adapting the definition addressed in [11].
Definition 1: Given the current control input set S ⊂ N with H ⊆ S and y i ∈ V z \ S, where V z is derived from G j ( 1 H ), j = 1, . . ., d (see (13) and ( 14)), C i, j |S ∈ N is the canalization number of y i on G j (S 1 ) with respect to 1, i.e., fixing is the summation of the canalization numbers of y i over j = 1, . . ., d. G j (S 1 ) represents the perturbed BN G j to which the control inputs S 1 are provided.Thus, C i, j |S specifies how many nodes are fixed to logic 1 by applying the additional control input y i = 1 to G j (S 1 ).Here, one must count all successive canalization effect caused by y i = 1 in computing C i, j |S , since more than one state transition equation can be fixed to 1 in a sequential way evolving from y i = 1.Also, only state variables belonging to V j z \ S are considered in deriving C i, j |S , as those in V j n already converge to 1 by the control inputs S. Clearly, the greater the value of C i, j |S is, the closer to 1 n the resultant attractor of G j is positioned when the additional control input y i = 1 is furnished.Since the favorable effect of the additional control input should be delivered over all perturbed BNs, it would be the best policy to recursively select the additional control inputs G := {φ m+1 , . . ., φ m * } in the order of C i |S .Once G is determined, H for G can be substituted by another solution set among under the assumption of | | ≥ 2. If a member of , any of which succeeds in global stabilization of G, has some entries belonging to G , employing such a solution set could reduce the number of overall control inputs.Thus, if γ ∈ exists for which γ ∩ G = ∅, the original H is replaced by γ .If there exist more than one such a solution set, γ having the greatest |γ ∩ G | is chosen as the control inputs for G.
Combining our discussions so far, we address in formal terms the algorithm for robust stabilization as follows.

For y i ∈ V z \ S, derive C i
|S by applying y i := 1 to each G j (S 1 ) and computing the canalization number C i, j |S .4. Among V z \ S, select an additional control input φ m+k such that The solution to the robust stabilization problem is given by H ∪ G .
As the following theorem attests, the abovementioned algorithm provides the solution to the robust stabilization problems for perturbed BNs G 1 , . . ., G d .
Theorem 1: Given the nominal BN G and perturbed BNs G j , j = 1, . . ., d, assume that the set of minimum FVSs ⊂ 2 N is derived with respect to the reduced BN H , where |γ | = m, ∀γ ∈ .Then, the control input set H ∪ G derived in Algorithm 1 solves the robust stabilization problem for G j 's with the maximum allowable number of control inputs m * > m.
Proof: First, assume that m * = m + 1.By Steps 4-6 of Algorithm 1, G = {φ m+1 } where applying y φ m+1 = 1 along with 1  H makes G 1 , . . ., G d converge to their attractors that have the smallest sum of Hamming distances with 1 n .Hence, H ∪ G solves the robust stabilization problem for G j 's.Next, assume that m * > m + 1.For proving by induction, suppose that at the end of the kth iteration, 1 ≤ k < m * − m, the additional control input set G := {φ m+1 , . . ., φ m+k } is derived by Algorithm 1 such that applying y φ m+i := 1, i = 1, . . ., k, along with 1 H makes G 1 , . . ., G d converge to their attractors yielding the smallest sum of Hamming distances with 1 n .Since k < m * − m, the iteration must continue by returning to Step 3 of Algorithm 1, where S = H ∪ G .Clearly, in view of Step 4, the next chosen control input y φ m+k+1 = 1 makes the resultant BNs converge to the attractors that are as close to 1 n as possible in terms of the summation of Hamming distances.Hence, G := {φ m+1 , . . ., φ m+k , φ m+k+1 } serves as the optimal additional control input set with the cardinality k + 1.This completes the proof.
Since derivation of FVSs in Step 1 is NP-complete [36], Algorithm 1 falls into a class of NP-completeness.However, as demonstrated in [30], the computational load for Step 1 is manageable even for complex BNs so long as the networks have realistic average in-degree-about two for biological networks [37].Also, since the size of V j z in the canalized network G j ( 1 H ) is much reduced compared with n, the computational load for obtaining C i |S in Step 3 is not serious either as will be validated in our case study.II.As boldfaced 0's in the attractors illustrate, the sum of the Hamming distance between the attractor of a perturbed BN and the desired one 1 10 is identical for two sets of control inputs.As this leads to average 85% stabilization of the state variables (17 out of 20), a significant improvement is obtained compared with the case of applying only H (65% stabilization as addressed in Example 2).The actual control inputs to the original BN F are obtained by the inverse transformation T −1 c (refer to (7) and ( 8)).
♦ The network graphs of G 1 ( 1 H ) and G 2 ( 1 H ) in Fig. 1 exemplify the efficiency of determining additional control inputs φ m+k .By Step 3 of Algorithm 1, we first extract V z , of which convergence to the desired value is distracted by the mutations ("nodes not fixed to 1" in Fig. 1).By Step 4, furthermore, we recursively search for the one whose designation as the additional control input contributes most significantly to the convergence of the remaining nodes to the desired values.This two-step stabilization strategy can greatly reduce the number of control inputs compared with other methods that do not consider the influence of mutations in deriving the control inputs.
Though only the summation of the Hamming distance is utilized in Algorithm 1, one can adjust the proximity measure in accordance with the detailed purpose of stabilizing control.For example, if the dispersion degree of the attractors around the desired one is critical, the associated measure such as the standard deviation can be included in Algorithm 1. Also, if the convergence of specific nodes is of great importance, e.g., apoptosis node in cancer cells, the proximity measure may be augmented by assigning much weights to those nodes.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

A. Random Complex BNs
For the purpose of evaluating the performance of the proposed method, we implement the algorithm and conduct numerical experiments using 4000 randomly generated BNs in which each ensemble of 1000 BNs contains 20, 30, 40, and 50 nodes, respectively.The authors' website [38] provides the Python codes of Algorithm 1 and random BNs generation.For all random BNs, we apply 1) the proposed method wherein one additional control input is derived besides H for the nominal BN, and 2) the original FVS control deriving the minimum FVS of the nominal BN for the sake of comparison.
The experimental results on robust stabilization of random BNs are shown in Fig. 2. Each ensemble of BNs has two kinds of perturbed BNs (d = 2) in each of which two nodes degenerate into adverse values by mutation."Score" Fig. 3. Average stabilization score for each ensemble of 1000 BNs having, respectively, 20, 30, 40, and 50 nodes and two kinds of perturbed BNs (d = 2).Instead of the original FVS, we apply all the possible subsets of the original FVS that have the same size as the solution of the proposed algorithm, and take their average scores.in Fig. 2(a) denotes the ratio of the number of successfully stabilized nodes to the total node number except mutated ones.If, for instance, 15 nodes are stabilized to desired values in a random BN with 20 nodes, the score is 15/(20 − 2) = 0.833.Also, if the original FVS method yields multiple solution sets, the maximum value among their scores is selected.
As Fig. 2 illustrates, the proposed control scheme exhibits greater performance than the original FVS method in terms of both score and computational time, even though the FVS method expends more control inputs than the proposed scheme.The performance enhancement becomes evident as the node number increases.In the case of BNs with 40 nodes, the proposed algorithm has a slightly better score (0.956) than the FVS method (0.921), but the number of its control inputs is less than that of the FVS method by more than one-3.252versus 4.558 as marked in Fig. 2(b).The superiority of the proposed algorithm over the FVS method is obvious when comparing their computational load.Referring to Fig. 2(c), it takes mere average 1.07 s to run the proposed algorithm, while the FVS method needs 12.316 s to complete its execution.
This tendency of performance improvement is even more striking when the node number increases to 50.While the stabilization score is comparable with each other (0.963 versus 0.931), the proposed algorithm shows much better performance in terms of the control input size and computational time: it needs only 3.392/4.893= 69.32% the number of control inputs employed by the FVS method and 1.547/100.269= 1.54% computational time.
For a comparative study, we investigate the stabilization outcome of the original FVS method under the constraint that its control input size is trimmed to that of the proposed algorithm.To this end, we derive average stabilization scores of all the subsets of the original FVS, which have the same cardinality as the solution of the proposed algorithm.Referring to Fig. 3, the proposed scheme shows greater stabilization scores than the trimmed original FVS method by about 19% over all types of random BNs.This result elicits that the proposed algorithm can serve as a proper robust stabilization strategy for accommodating practical constraint on the number of control inputs.Figs. 2 and 3 affirm that the proposed algorithm derives efficient constant controls for solving the robust stabilization problem of large-scale BNs.We expect from these results that in comparison with the FVS method, the proposed algorithm will retain its performance and scalability as the complexity of BNs escalates further.The latter attests to the capability of the proposed scheme that depending on the practical constraint on the considered biological network, it can present various control input sets with incremental cardinalities and performance.

B. Biological Networks
To validate the applicability of the proposed scheme, we consider the robust stabilization problem of two real biological systems-an influence network representing the metastatic process of cells [39] and the mitogen-activated protein kinase (MAPK) signaling network on cancer cell fate decision [40].
1) Metastasis Influence Network: Fig. 4 is the connectivity graph of the metastasis influence network adapted from [39].Rectangles in Fig. 4 symbolize ingredient nodes corresponding to biochemical species (proteins, microRNAs (miRNAs), processes, and so on), and edges indicate activating (red) or inhibitory (blue) influences of one node onto others; Table VII in Appendix A shows Boolean logical rules of ingredient nodes in the metastasis influence network.This BN has total 32 state variables (n = 32) among which there exist two input nodes ECMicroenv and DNAdamage, and three output nodes Metastasis, CellCycleArrest, and Apoptosis.Input and output nodes are positioned in the upper-and downstream of the BN in Fig. 4, respectively.Here, ECMicroenv represents whether the effect of the extracellular microenvironment turns on or off.Furthermore, DNAdamage indicates whether or not a DNA damage occurs to the cell.
Out of nine attractors of this BN [39], we designate Apoptotic Attractor 3 as the desired attractor c, as it corresponds to the state that the metastasis influence network converges to programed cell death (Apoptosis = 1), while metastasis is prevented (Metastasis = 0).The first row entitled "Nominal BN" in Fig. 5 shows the binary values of Apoptotic Attractor 3, wherein the external inputs are set to be DNAdamage = 1 and ECMicroenv = 1 ("ECM" in Fig. 5).
Among various mutants of the metastasis influence network identified in [39], we designate MCF7 and T47D as perturbed BNs (d = 2), respectively, of which information is summarized in Table III.TWIST1 is excited by mutation in MCF7, and AKT2 in T47D.Both mutants lead to breast cancer with the Carcinoma histological type.Detailed biological implications of TWIST1 and AKT2 genes are addressed in the Supplementary Material.
When the aforementioned mutations are applied to the nominal metastasis influence network, the resultant BNs, G 1 for Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Binary values of the desired attractor, Apoptotic Attractor 3, of the metastasis influence network, the degraded attractors of perturbed BN 1 (TWIST1 = 1) and perturbed BN 2 (AKT2 = 1) [39, Table S4], and stabilized attractors of perturbed BNs.MCF7 and G 2 for T47D, yield long-term behaviors as shown in "Perturbed BN 1" and "Perturbed BN 2" rows of Fig. 5. Fractions in the entries are average activation ratios of the genes.For instance, Dkk1 = 0.53 in perturbed BN 1 implies that Dkk1 gene converges to 1 (including 1's in cycles) from 53% all the possible initial states of perturbed BN 1. Referring to Fig. 5, perturbed BNs 1 and 2 have stabilization score 7/31 = 0.226 and 8/31 = 0.258, respectively, in regard to desired attractor c.Clearly, both mutations inflict serious degradation to long-term behaviors of the networks.In particular, both perturbed BNs are not able to programed cell as Apoptosis node is upset to 1, and perturbed BN 2 even fails to suppress metastasis, as Metastasis node is forced to turn on.
Table IV shows the result of robust stabilization for the population of the metastasis influence network with two perturbed BNs.Like the previous numerical experiments on random BNs, one additional control input is allowed to be used besides those ones for the nominal BN.For a comparative study, we apply to the same problem the original FVS control, the stable motif method [41], and the target control based on the logical domain of influence (LDOI) [34].The proposed method results in one minimum FVS H = {miR200, p53} for the nominal BN and the additional control input G = {CDH1}.The original FVS control produces a single control input set with eight elements, and the target control based on LDOI results in eight control input sets with one to four elements.Though the original FVS control indeed yields two control input sets, one of them turns out infeasible as it contains the mutated node TWIST1.Similarly, we discard all the solution sets of the LDOI method containing at least one mutated node.On the other hand, the stable motif method fails to derive any solution set even after exhaustive runs.
The superiority of the proposed method over the original FVS control and the other methods is obvious.Compared with the FVS control, the proposed method yields only 3/8 = 37.5% control inputs and takes almost the same computational time, while achieving 3.6% higher performance in terms of the maximum stabilization score.As Jaccard similarity between two solution sets is merely 0.22, it can be said that the proposed scheme succeeds in deriving different combinations of control inputs that are suitable for tackling the robust stabilization problem at hand.The performance improvement is more striking when evaluated in terms of the normalized stabilization score, which is obtained by applying the same number of control inputs as the proposed method (three) among the derived solution set of the FVS control (eight) and by averaging the results.While the normalized score of the FVS control is 0.661, that of the proposed method is 0.919 (identical to the maximum score), meaning that the proposed method produces the stabilization score 39.0% higher than the FVS control under the same condition on the number of control inputs.
In comparison with the LDOI method, the proposed scheme also shows greater performance in terms of both robustness and computation time-16.3% and 26.1% higher Fig. 6.MAPK signaling network modeled by a BN [40], where red edges with arrows represent activation and blue ones with bars inhibition.performances in terms of the maximum and normalized score, respectively, while taking only 63.3% computation time.Furthermore, it is analyzed that the failure of the stable motif method stems from the high density of the metastasis influence network which impedes the derivation of feasible sets of stable motifs in the application of the method.
Binary values of the attractor in each stabilized perturbed BN are shown in the last two rows of Fig. 5.One can see that perturbed BN 1 successfully inhibits metastasis, while perturbed BN 2 not only inhibits metastasis but also regains programed cell death.The performance of robust stabilization may be enhanced further at the cost of increasing the number of additional control inputs.It turns out that when three additional control inputs are allowed, the proposed scheme accomplishes 100% robust stabilization, namely, all the state variables except mutated ones converge to the desired values in both perturbed BNs.The corresponding additional control inputs are derived as {CDH1, CDH2, miR34}.
Remark 2: Unlike the case of random BNs, it takes slightly more time than the FVS control to run the proposed scheme for the metastasis influence network (12.6 versus 12.1 [s]).This protraction is caused by exceptionally high density of the metastasis influence network as mentioned before, since its average in-degree is 5.2, and among 32 state variables, seven nodes have the maximum in-degree 8, whereas random BNs employed in the previous experiment have average indegree 1.64.In particular, high in-degree of a BN gives rise to much computational burden on transforming the state transition equation of G into the minimum SOP form.We find that among total computation time of 12.6 s, about 12 s is used for deriving the minimum SOP form.While we utilize SymPy Python package [42] in implementing the transformation (refer to the authors' website [38] for source codes), it is expected that the computation time may be curtailed by applying more efficient algorithms.
2) MAPK Signaling Network: Fig. 6 is the connectivity graph of the MAPK signaling network taken from [39].The corresponding Boolean logical rules of ingredient nodes are provided in Table VIII of Appendix A. This BN has total 53 state variables (n = 53) among which there exist four inputs nodes and three output nodes that are positioned in the upper-and downstream of the BN in Fig. 6, respectively.
In this experiment, we take as the desired attractor a fixed point of the MAPK signaling network whose binary values are shown in the first row of Fig. 7.Note that the phenotype of this attractor is set in favor of inducing proper cancer cell fate decision as Proliferation = 0, Apoptosis = 1, and Fig. 7. Binary values of the desired attractor of the MAPK signaling network, the degraded attractors of perturbed BN 1 (AKT = 1, ERK = 1) and perturbed BN 2 (FRS2 = 1, p53 = 0) [40], and stabilized attractors of perturbed BNs. and 639-V as perturbed BNs (d = 2), each of which contains two mutations and causes bladder cancer.For detailed biological implications of each mutant, refer to the Supplementary Material.The stable state behavior of the perturbed BNs is summarized in the second and third rows of Fig. 7.
Table VI shows the result of the numerical experiment for the MAPK signaling networks.Similar to the first experiment, the proposed scheme produces a better maximum stabilization score than the other methods with less or compatible computational time.In terms of the normalized score (i.e., having the same number of control inputs), the superiority of the proposed scheme is more remarkable-14.7%,16.3%, and 5.9% higher than, respectively, the FVS control, stable motif method, and target control based on LDOI.Especially, the computation time of the proposed scheme is only 24.3/812.1 = 3.0% that of the FVS control.It is observed that the MAPK signaling network has more nodes than the previous metastasis influence network (53 versus 32), while its average in-degree is much lower (2.1 versus 5.2).Hence, it can be said that the execution speed of the proposed scheme is dominated mainly by the average in-degree of the considered BN.

V. CONCLUSION
This article addresses a novel robust stabilizing control strategy aiming at complex biological systems based on logical models with uncertainties.Throughout the attractor-specific network reduction and FVS control, we search for additional control inputs which, along with those ones for the nominal BN, take all the perturbed BNs to their own attractors as near as possible to the desired fixed point measured in terms of the Hamming distance.The practicality of the proposed scheme is remarkable since it requires neither algebraic nor structural constraint on the considered BN.The algorithm for deriving complete control inputs needs exponential complexity in theory with respect to the number of nodes of BNs.But, we ascertain that it can still be applied to complex BNs, since the actual computational load is manageable as verified in the presented numerical experiments on complex random BNs and a real biological system with high in-degrees.
Though the Hamming distance as the proximity measure embodies primary biological similarity of the stabilized perturbed BN to the desired long-term behavior, it can be replaced by other criteria according to characteristics of the underlying biological system as mentioned before.For instance, one can assign unbalanced high scores on output nodes, which determine key phenotypes of the cell, thus pursuing output stabilization or target control [27], [34] subject to mutations.The latter problem will serve as an interesting topic for a future study.

APPENDIX BOOLEAN LOGICAL RULES OF BIOLOGICAL NETWORKS
We tabulate, in Tables VII and VIII, Boolean logical rules of the activity level of each node in the metastasis influence network and MAPK signaling network, respectively.

Fig. 1 .
Fig. 1.Illustration of applying the proposed scheme to G in Example 1. Robust Stabilization Problem for BNs: Given a population of the nominal BN G and perturbed BNs G 1 , . . ., G d , let ⊂ 2 N be the set of minimum FVSs derived with respect to the reduced BN H , where |γ | = m, ∀γ ∈ , and let m * , m < m * < n, be the maximum allowable number of control inputs.Then, determine 1) H := {φ 1 , . . ., φ m } ∈ , the set of m control inputs for G, and 2) m * − m additional control inputs

Example 4 :
We apply Algorithm 1 to the robust stabilization problem for the population of BNs G in Example 1 in which two perturbed BNs G 1 and G 2 exist as addressed in Example 2. Assume that the maximum allowable number of control inputs is m * = 3.As H = {2, 9} with | H | = m = 2, one additional control input is available.Since Steps 1 and 2 of Algorithm 1 are already undertaken in Examples 1-3, let us turn to Step 3. From Example 3, and C 1,2 | H = 3 as fixing y 1 := 1 canalizes y 3 to 1, which in turn canalizes y 10 to 1; refer to (12) and the network graph of G 2 ( 1 H ) in the right part of Fig. 1.According to the result shown in Table I, we select G := {7} or {10} as the additional control input by Step 4 of

Algorithm 1 .
Since m * − m = 1 and | | = 1, Steps 5 and 6 are not activated.The final solution to the robust stabilization problem is H ∪ G = {2, 9, 7} or {2, 9, 10}.The resultant attractors of perturbed BNs with respect to each set of control inputs are shown in Table

Fig. 2 .
Fig. 2. Experimental results on robust stabilization of complex random BNs.(a) Average stabilization score for 4000 random BNs, each ensemble of 1000 BNs having, respectively, 20, 30, 40, and 50 nodes and two kinds of perturbed BNs (d = 2).(b) Average number of elements in the derived control input sets.(c) Average computational time required to derive the solutions for each ensemble of random BNs.The error bars in (a)-(c) denote the 95% confidence intervals of the average score, number of elements, and computation time, respectively; refer to [38] for their values.

Fig. 4 .
Fig.4.Metastasis influence network modeled by a BN[39], where red edges with arrows represent activation and blue ones with bars inhibition.

TABLE III INFORMATION
OF PERTURBED BNs OF THE METASTASIS INFLUENCE NETWORK

TABLE IV RESULTS
OF ROBUST STABILIZATION FOR THE METASTASIS INFLUENCE NETWORK (n = 32 AND AVERAGE IN-DEGREE = 5.2)

TABLE V INFORMATION
OF PERTURBED BNs OF THE MAPK SIGNALING NETWORK

TABLE VI RESULTS
OF ROBUST STABILIZATION FOR THE MAPK SIGNALING NETWORK (n = 53 AND AVERAGE IN-DEGREE = 2.1)

TABLE VII BOOLEAN
LOGICAL RULES OF THE METASTASIS INFLUENCE NETWORKGrowth_Arrest = 1. to V, designate

TABLE VIII BOOLEAN
LOGICAL RULES OF THE MAPK SIGNALING NETWORK