{'fraction_non_alphanumeric': 0.0355254115748902, 'fraction_numerical': 0.009531207983507126, 'mean_word_length': 4.436972059779078, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Dynamics-based Approach for the Target Control of Boolean Networks

We study the target control problem of asynchronous Boolean networks, to identify a set of nodes, the perturbation of which can drive the dynamics of the network from any initial state to the desired steady state (or attractor). We are particularly interested in temporary perturbations, which are applied for sufficient time and then released to retrieve the original dynamics. Temporary perturbations have the apparent advantage of averting unforeseen consequences, which might be induced by permanent perturbations. Despite the infamous state-space explosion problem, in this work, we develop an efficient method to compute the temporary target control for a given target attractor of a Boolean network. We apply our method to a number of real-life biological networks and compare its performance with the stable motif-based control method to demonstrate its efficacy and efficiency.


INTRODUCTION
Cell reprogramming has garnering attention for its therapeutic potential for treating the most devastating diseases characterised by diseased cells or a deficiency of certain cells. It is capable of reprogramming any kind of abundant cells in the body into the desired deficient cells to restore functions of the diseased organ [10,11,38]. It has shown promising benefits for clinical applications, such as cell and tissue engineering, regenerative medicine and drug discovery.
In their seminal work, Yamanaka et al. showed that human somatic cells can be converted to induced pluripotent stem cells (iP-SCs) by a cocktail of defined factors [40,44]. The generated iPSCs have the ability to further propagate and differentiate into many cell types. However, the application of iPSC reprogramming is often restricted, due to that: (1) the generated iPSCs have a risk of cancerous tumour formation [10,11]; (2) the iPSC reprogramming and differentiation process usually requires long time to produce sufficient cells for application, which leads to a significant experimental cost [11]; and (3) the iPSCs often encounter cell cycle arrest after differentiation, which makes it impossible to expand the number of cells for therapeutic transplantation [10]. The limitations of iPSC reprogramming reinforce the need of direct reprogramming, also called transdifferentiation. Direct reprogramming harnesses abundant somatic cells to regenerate defective cells by reprogramming the somatic cells directly into the desired cell type bypassing the pluripotent state. As a consequence, direct reprogramming can not only reduce the risk of tumourigenesis and teratoma formation, but also shorten the period of time for producing enough desired cells for therapeutic application.
A major challenge of cell reprogramming lies in the identification of effective target proteins or genes, the manipulation of which can trigger desired changes. Lengthy time commitment and high cost hinder the efficiency of experimental approaches, which perform brute-force tests of tunable parameters and record corresponding results [42]. This strongly motivates us to turn to mathematical modelling of biological systems, which allows us to identify key genes or pathways that can trigger desired changes using computational methods. Boolean network, first introduced by Kauffman [15], is a well-established modelling framework for gene regulatory networks and their associated signalling pathways, and it has apparent advantages compared to other modelling frameworks [1]. Boolean network provides a qualitative description of biological systems and thus evades the parametrisation problem, which often occurs in quantitative models, such as models of ordinary differential equations (ODEs). In Boolean networks, molecular species, such as genes and transcription factors, are described as Boolean variables. Each variable is assigned with a Boolean function, which determines the evolution of the node. Boolean functions characterise activation or inhibition regulations between molecular species. The dynamics of a Boolean network is assumed to evolve in discrete time steps, moving from one state to the next, under one of the updating schemes, such as synchronous or asynchronous. Under the synchronous scheme, all the nodes update their values simultaneously at each time step; while under the asynchronous scheme, only one node is randomly selected to update its value at each time step. We focus on the asynchronous updating scheme since it can capture the phenomenon that biological processes occur at different time scales. The steady-state behaviour of the dynamics is described as attractors, to one of which the system eventually settles down. Attractors are hypothesised to characterise cellular phenotypes [14]. Each attractor has a weak basin and a strong basin. The weak basin contains all the states that can reach this attractor, while the strong basin includes the states that can only reach this attractor and cannot reach any other attractors of the network. In the context of Boolean networks, cell reprogramming is interpreted as a control problem: modifying the parameters of a network to lead its dynamics towards a desired attractor.
Control theories have been employed to modulate the dynamics of complex networks in recent years. Due to the intrinsic nonlinearity of biological systems, control methods designed for linear systems, such as structure-based control methods [5,8,19], are not applicable -they can both overshoot and undershoot the number of control nodes for non-linear networks [9]. For nonlinear systems of ODEs, Fiedler et al. proved that the control of a feedback vertex set is sufficient to control the entire network [7,26,47]; and Cornelius et al. proposed a simulation-based method to predict instantaneous perturbations that can reprogram a cell from an undesired phenotype to a desired one. However, further study is required to figure out if these two methods can be lifted to control Boolean networks. Several methods based on semi-tensor product (STP) have been proposed to solve different control problems for Boolean control networks (BCNs) under the synchronous updating scheme [3,18,21,43,45,48,49,51]. For synchronous Boolean networks, Kim et al. developed a method to compute a small fraction of nodes, called 'control kernels', that can be modulated to govern the dynamics of the network [16]; and Moradi el al. developed an algorithm guided by forward dynamic programming to solve the control problem. However, all these methods are not directly applicable to asynchronous Boolean networks. To tackle this problem, we have developed several decomposition-based methods, which exploit both the structural and dynamical information, to cope with source-target control with instantaneous, temporary and permanent perturbations [22,23,31,32,39] and target control with instantaneous perturbations [2] for asynchronous Boolean networks. In view of the difficulties and expenses in conducting biological experiments, our methods compute the minimal control sets, which can be easily translated for wet-lab validation.
Cells in tissues and in culture normally exist as a population of cells, corresponding to different stable steady states [36]. There is a need of target control methods to compute a subset of nodes, the control of which can always drive the system from any initial state to a desired target attractor. The target control method developed in our previous work [2] adopts instantaneous perturbations, that are only applied instantaneously, but at a cost, rather larger number of control nodes are required than temporary and permanent perturbations [39]. Moreover, it is difficult to guarantee that all the perturbations take effect at the same time in biological experiments. Thus, target control with temporary perturbations is more appealing.
In this paper, we develop a target control method with temporary perturbations for asynchronous Boolean networks. Our idea is to find a control C = (0, 1), which is a tuple of two sets, such that the application of C -setting the value of a node, whose index is in 0 (or 1), to 0 (or 1) -can drive the network from any initial state s in the state space S to an intermediate state s ′ in the weak basin of the target attractor. We hold the control C for sufficient time and let the network evolve to a state in the strong basin of the target attractor. After that, the control can be released and the network will eventually and surely reach the target attractor. Since the network can take any state s ∈ S as an initial state, the possible intermediate states form a subset S ′ of S, called schema. According to our previous work [39], we know that all the intermediate states should fall into the weak basin of the target attractor. Therefore, we partition the weak basin into a set of mutually disjoint schemata. Each schema results in a candidate control, which is further minimised and verified. Clinical applications are highly time-sensitive, controlling more nodes may shorten the period of time for generating sufficient desired cells [11]. Hence, we integrate our method with a threshold ζ on the number of perturbations. By increasing ζ , we can obtain solutions with at most ζ perturbations. It is worth noting that more perturbations may cause a significant increase in the experimental cost, hence, the parameter ζ should be considered individually based on specific experimental settings.
We have implemented our method and compared its performance with the stable motif-based control (SMC) [46] on various real-life biological networks, as both methods focus on temporary target control of asynchronous Boolean networks. The results show that our method outperforms SMC in terms of the computational time for most of the networks. Both methods find a number of valid temporary controls, but our method is able to identify more controls with fewer perturbations for some networks. Another interesting observation is that the number of required perturbations is often quite small compared to the sizes of the networks. This agrees with the empirical findings that the control of few nodes can reprogram biological networks [27].

BACKGROUND AND NOTATIONS
In this section, we give preliminary notions of Boolean networks. Let [n] denote the set of positive integers {1, 2, . . . , n}.

Boolean networks
A Boolean network (BN) describes elements of a dynamical system with binary-valued nodes and interactions between elements with Boolean functions. It is formally defined as: represents the value that x i takes when the network is in state s.
For two states s, s ′ ∈ S, the Hamming distance between s and s ′ is denoted as hd(s, s ′ ).
and 0 and 1 are mutually disjoint (possibly empty) sets of indices of nodes of a Boolean network G. The size of the control C is defined as |C | = |0| + |1|. Given a state s ∈ S, the application of C to s, denoted as C(s), is defined as a state s ′ ∈ S, such that The control can be lifted to a subset of states S ′ ⊆ S. Given a control C = (0, 1), C(S ′ ) = S ′′ , where S ′′ = {s ′′ ∈ S |s ′′ = C(s ′ ), s ′ ∈ S ′ }. S ′′ includes all the intermediate states with respect to C. The application of C results in a new Boolean network, defined as follows.
The state space of G | C , denoted S | C , is derived by fixing the values of the variables in the set C to their respective values and is defined as

Dynamics of Boolean networks
In this section and the next section, we define several notions that can be interpreted on both G and G | C . We use the generic notion G = (X , F ) to represent either G = (X , F ) or G | C = (X ,F ). We assume that a Boolean network G = (X , F ) evolves in discrete time steps. It starts in an initial state s 0 and its state changes in every time step according to the update functions F . Different updating schemes lead to different dynamics of the network [25,50]. In this work, we are interested in the asynchronous updating scheme as it allows biological processes to happen at different classes of time scales and thus is more realistic.
The asynchronous evolution of G is a function ξ : It is worth noting that the asynchronous dynamics is non-deterministic. At each time step, only one node is randomly selected to update its value and a different choice may lead to a different next state s ′ ∈ ξ (j + 1). Henceforth, when we talk about the dynamics of G, we shall mean the asynchronous dynamics. The dynamics of a Boolean network can be described as a transition system (TS). Definition 2.4 (Transition system of Boolean networks). The transition system of a Boolean network G, denoted as TS, is a tuple (S, E), where the vertices are the set of states S and for any two states s and s ′ there is a directed edge from s to s ′ , denoted s → s ′ iff s ′ is a possible next state of s according to the asynchronous evolution function ξ of G.
Similarly, we denote the transition system of a Boolean network under control, G | C , as TS| C .

Attractors and basins
A path ρ from a state s to a state s ′ is a (possibly empty) sequence of transitions from s to s ′ in TS, denoted ρ = s → s 1 → . . . → s ′ . A path from a state s to a subset S ′ of S is a path from s to any state s ′ ∈ S ′ . An infinite path ρ from s, ρ = s → s 1 → . . ., is a sequence of infinite transitions from s. A state s ′ ∈ S appears infinitely often in ρ if for any i ≥ 0, there exists j ≥ i such that s j = s ′ . We assume every infinite path ρ is fair -for any state s ′ that appears infinitely often in ρ, every possible next state s ′′ of s ′ also appears infinitely often in ρ. For a state s ∈ S, reach(s) denotes the set of states s ′ such that there is a path from s to s ′ in TS. Definition 2.6 (Attractor). An attractor A of TS (or of G) is a minimal non-empty subset of states of S such that for every state s ∈ A, reach(s) = A.
Attractors are hypothesised to characterise the steady-state behaviour of the network. Any state which is not part of an attractor is a transient state. An attractor A of TS is said to be reachable from a state s if reach(s) ∩A ∅. The network starting at any initial state s 0 ∈ S will eventually end up in one of the attractors of TS and remain there forever unless perturbed. Under asynchronous updating scheme, there are singleton attractors and cyclic attractors. Cyclic attractors can be further classified into: (1) a simple loop, in which all the states form a loop and every state appears only once per traversal through the loop; and (2) a complex loop, which has intricate topology and includes several loops. Fig. 2 (a), (b) and (c) show a singleton attractor, a simple loop and a complex loop, respectively. Let A denote all the attractors of TS. For an attractor A, A ∈ A, we define its weak basin as bas W  Fig. 1(c). For attractor A 1 , its strong basin bas S TS (A 1 ) = {000, 001} is shown as the shaded grey region; its weak basin contains six states, i.e. bas W TS (A 1 ) = {000, 001, 101, 011, 100, 010}. We can see that only attractor A 1 is preserved in TS| C in Fig. 1(d).

The control problem
As described in the introduction, many devastating diseases, such as Parkinson's disease and Alzheimer's disease, are caused by a deficiency of particular cells. Cell reprogramming can transform abounding somatic cells into the desired cell type. In the context of Boolean networks, this process is, indeed, stirring the dynamics of the network from a source attractor to a desired target attractor. However, cells in culture and in situ are usually not isolated but exist in a population consisting of various cell phenotypes or even transient cell states. Hence, it is important to develop a target control method to identify key nodes that can guide the network towards a desired target attractor from any other distinct steady states or transient states. This can be defined as a target control problem: given a Boolean network G and a target attractor A t , finding a control C, the application of which can drive the network from any source state s ∈ S to A t . When the source state s is fixed, finding a control C to drive the network from s to A t , is a source-target control problem. Based on the application time of control, we have: (1) temporary control -perturbations are applied for a finite (possibly zero) number of steps and then released; (2) permanent control -perturbations are applied for all the following steps. When perturbations are applied instantaneously, we call it instantaneous control, which is a special case of temporary control. Temporary control has shown its apparent advantages in reducing the number of perturbations [39], thus in this work, we focus on temporary target control, formally defined as follows.
Definition 2.8 (Temporary target control). A temporary target control is a control C = (0, 1), such that there exists a t 0 > 0, for all t > t 0 , the network always reaches the target attractor A t on the application of C to any source state s ∈ S for t steps.

RESULTS
In this section, we shall develop a method to solve the temporary target control problem. First, we introduce the following lemma, which is crucial for the development of the method. Instead of presenting a formal proof for Lemma 3.1, we give an intuitive explanation below. Definition 2.3 shows that the application of a control C results in a new Boolean network G | C and the state space is restricted to S | C . To guarantee the inevitable reachability of A t , by the time we release the control, the network has to reach a state s in the strong basin of A t w.r.t. the original transition system TS, i.e. bas S TS (A t ), from which there only exist paths to A t . This requires the remaining strong basin in S | C , i.e. (bas S TS (A t ) ∩ S | C ), is a non-empty set; otherwise, it is not guaranteed to reach A t . Furthermore, the condition C(S) ⊆ bas S TS | C (bas S TS (A t ) ∩ S | C ) ensures any possible intermediate state s ′ ∈ C(S) is in the strong basin of the remaining strong basin (bas S TS (A t ) ∩ S | C ) in the transition system under control TS| C , so that the network will always evolve to the remaining strong basin. Once the network reaches the remaining strong basin, the control can be released and the network will evolve spontaneously towards the target attractor A t . Based on the definition of the weak basin, it is sufficient to search the weak basin bas W TS (A t ) for temporary target control. A noteworthy point is that temporary control needs to be released once the network reaches a state in (bas S TS (A t ) ∩ S | C ). On one hand, Lemma 3.1 guarantees that partial strong basin of A t in TS is preserved in TS| C , while it does not guarantee the presence of A t in TS| C . In that case, the control C has to be released at one point to recover the original TS, which at the same time retrieves A t . On the other hand, in clinic, it is preferable to eliminate human interventions to avoid unforeseen consequences. Concerning the timing to release the control, since it is hard to interpret theoretical time steps in diverse biological experiments, it would be more feasible for biologists to estimate the timing based on empirical knowledge and specific experimental settings.
Previously, we have developed efficient decomposition-based algorithms to compute the exact basins of an attractor, which exploit both the structural and dynamical properties of the network [31,32]. In the algorithm we develop here, we shall use these procedures to compute the weak basin and the strong basin of an attractor and refer them as Comp_WB and Comp_SB, respectively. Next, we define the projection of a state s ∈ S to a subset B of [n], which represents the indices of a subset of nodes X ′ ⊆ X as follows.  Intuitively, for any node x i , i ∈ 0, it has a value of 0 in any state s ∈ S ′ ; for any node x i , i ∈ 1, it has a value of 1 in any state s ∈ S ′ . The projection of S ′ to the don't-care-set D contains all combinations of binary strings with |D| bits. Thus, any schema S ′ is of size 2 |D| . Since the total number of nodes n = |0| + |1| + |D| is fixed, a larger schema implies more elements in D and fewer elements in 0 ∪ 1.

18:
if C j i Ω then // C i has not been checked.

21:
if isValid = true then 22: add C j i to L, ζ := min(ζ , |C j i |) 23: The notion of schema leads the way to find temporary target control. Each schema W i of the weak basin bas W TS (A t ) gives a candidate temporary target control C i = (0 i , 1 i ) for further optimisation and validation. A larger schema results in a smaller control set. To explore the entire weak basin bas W TS (A t ), we partition it into a set of mutually disjoint schemata W = {W 1 ,W 2 , . . . ,W m }, . For W i , the indices of its support variables in 0 i and 1 i form a candidate control C i = (0 i , 1 i ). Each candidate control C i is primarily optimised based on the properties of input nodes. Because input nodes do not have any predecessors, it is reasonable to assume that specified input nodes I s are redundant control nodes, while non-specified input nodes I ns are essential for control. For the remaining noninput nodes in C i , denoted C r i , we verify its subsets of size k based on Lemma 3.1 from k = 0 with an increment of 1, until we find a valid solution.
To further improve the efficiency of our method, we use binary decision diagram (BDD) as a symbolic representation of large state space. The size of a BDD is determined by both the set of states being represented and the chosen ordering of the variables. In BDD, a schema is represented as a cube and each state is the smallest cube, also called a minterm. To compute the largest schema S i of S is equivalent to compute the largest cube of S. The partitioning of the weak basin into schemata is then transformed into a cube cover problem in BDD. A different variable ordering may lead to a different partitioning. Given a fixed ordering, the partitioning remains the same. Although finding the best variable ordering is NP-hard, there exist efficient heuristics to find the optimal ordering. In this work, we compute a partitioning under one variable ordering as provided by the CUDD package [37] and compute the smallest subsets of candidate controls that are valid temporary target control sets.
Algorithm 1 implements the idea in pseudo-code. It takes as inputs the Boolean network G = (X , F ) and the target attractor A t . It first initialises two vectors L and Ω to store valid controls and the checked controls, respectively. (We use Ω to avoid duplicate control validations.) Then, it computes input nodes I and the non-specified input nodes I ns , I ns ⊆ I (line 3). The strong basin SB and the weak basin WB of A t of TS are computed using the decomposition-based procedures Comp_SB and Comp_WB developed in [31,32] (lines 4-5). The weak basin WB is then partitioned into m mutually disjoint schemata with procedure Comp_schemata. Realisation of this procedure relies on the function to compute the largest cube provided by the CUDD package [37]. For each schema W i , the indices of its support variables computed by procedure Comp_support_variables form a candidate control C i (line 11). The essential control nodes C e i of C i consist of the non-specified input nodes and the non-input nodes in C i constitute a set C r i for further optimisation (line 12). We search for the minimal subsets of C r i starting from size k = 0 with an increment of 1 and verify whether the union of a subset C sub j of C r i and the essential nodes C e i , namely C j i = C sub j ∪ C e i , is a valid temporary target control using procedure Verify_TTC in Algorithm 2. If C j i is valid, save it to L. When all the subsets have been traversed or a valid control has been found, we proceed to the next schema W i+1 . In the end, all the verified temporary target controls are returned.
The most time-consuming part of our method lies in the verification process. As shown in Algorithm 2, for each candidate control C, we need to reconstruct the associated transition relations F | C and compute the strong basin of the remaining strong basin in TS| C , i.e. bas S TS | C (SB| C ) (lines 6 and 7 of Algorithm 2). Even though we have developed an efficient method for basin computation, the computational time of Algorithm 2 still increases when the network size grows. To improve the efficiency, we propose two heuristics: Algorithm 1 is easily adapted to solve target control problem with instantaneous perturbations by focusing on the schemata of the strong basin of A t . In this way, we don't need to use Algorithm 2 for additional verification and the indices of support variables of each schema form an instantaneous control.

EVALUATION
Our temporary target control method, described in Algorithms 1 and 2, is implemented in the tool ASSA-PBN [25] based on the model checker MCMAS [20] to encode Boolean networks into the efficient data structure BDD. All the experiments are performed on a high-performance computing (HPC) platform, which contains CPUs of Intel Xeon Gold 6132 @2.6 GHz.
As discussed in the introduction, both our method (TTC) and the stable motif-based control (SMC) [46] focus on temporary target control of asynchronous Boolean networks. We apply our method on several real-life biological networks and compare its performance with SMC. Here we give a brief description on the networks. An overview of the networks can be found in Table 1.
• The myeloid differentiation network is designed to model myeloid differentiation from common myeloid progenitors to four cell types, including megakaryocytes, erythrocytes, granulocytes and monocytes [17]. • The apoptosis network consists of necessary pro-apoptotic and anti-apoptotic pathways to capture decision-making on cell survival or apoptosis [41]. • The cardiac gene regulatory network integrates major genes that play important roles in early cardiac development and FHF/SHF determination [13]. • The ERBB receptor-regulated G1/S transition protein network combines ERBB signalling with G1/S transition of the mammalian cell cycle to identify new targets for breast cancer treatment [34]. • The HSPC-MSC network describes intercommunication pathways between hematopoietic stem and progenitor cells (HSPCs) and mesenchymal stromal cells (MSCs) in bone marrow (BM) [6]. • The PC12 cell network models temporal sequence of protein signalling, transcriptional responses and subsequent autocrine feedbacks [29]. • The network of hematopoietic cell specification covers major transcription factors and signalling pathways for lymphoid and myeloid development [4]. • The bladder cancer network allows us to identify deregulated pathways and their influence on bladder tumourigenesis [33]. • The MAPK network is constructed to study MAPK responses to different stimuli and their contributions to cell fates [12]. • The model of HGF-induced keratinocyte migration captures the onset and maintenance of hepatocyte growth factorinduced migration of primary human keratinocytes [35]. • The Th-cell differentiation network models regulatory elements and signalling pathways controlling Th-cell differentiation [28].  Table 1: An overview of the networks and a comparison of the two methods (TTC and SMC). Symbol '-' means that the method failed to finish the computation within five hours.
• The HIV-1 network models dynamic interactions between human immunodeficiency virus type 1 (HIV-1) proteins and human signal-transduction pathways that are essential for activation of CD4+ T lymphocytes [30].
Attractors of the networks. Before the computation of target control, attractors are identified with our decomposition-based attractor detection method [24] and SMC, respectively. Our method identifies all the exact attractors (the number of states and the structures for both singleton and cyclic attractors) introduced in Section 2.3, while SMC identifies exact singleton attractors and quasi-attractors, which correspond to cyclic attractors. A quasiattractor can be considered as a superset of an attractor: the values of oscillate nodes in the corresponding attractor are not specified in a quasi-attractor. Columns 4-5 and 6-7 of Table 1 show the number of attractors computed by the two methods. Most of the attractors identified by the two methods are the same except for the cyclic attractor of the apoptosis network (marked in bold in Table 1). SMC identifies its quasi-attractor, which consists of 64 states, while the corresponding cyclic attractor has 56 states. Columns 8 and 9 of Table 1 show the execution time for attractor detection. We can see that our attractor detection method is more efficient than SMC.
Effectiveness. We compute temporary target control for each attractor of the networks with TTC and SMC. Since neither of the methods guarantees the minimal control, they may find control sets of different sizes for one attractor. For comparison, we only consider the smallest control sets. Fig. 3 (a) shows the number of smallest control sets for the myeloid differentiation network. The blue bars and grey bars represent the control sets that only appear in the results of TTC and SMC, respectively. The green bars denote the intersection of the two methods. The equation above each bar |C | = k describes the size of control sets. For attractors A 1 , A 5 and A 6 , TTC identifies smaller control sets than SMC. Taking A 5 , as an example, the minimal number of perturbations required by TTC and SMC is 3 and 4, respectively. Since we only consider the smallest controls, SMC identifies no control sets of size 3, thus we can only see the blue bar for this case. For other attractors (A 2 , A 3 and A 4 ) of the myeloid network in Fig. 3 (a) and all the attractors of the cardiac network in Fig. 3 (b), two methods require the same number of perturbations, but our method has the potential to identify more solutions than SMC.
For the other networks listed in Table 1, we summarise the number of control sets for one of the attractors in Fig. 3 (c). It shows that our method is able to identify smaller control sets than SMC for the bladder cancer network (SMC failed to compute results for hematopoiesis, HGF and HIV-1 networks). Our method also has the capability to provide more solutions, which may give more flexibility for clinical applications. Another interesting observation is that even for large networks, the number of perturbations is relatively small. Now we use the myeloid differentiation network as an example to show the consistency of our results with biological conclusions in [17]. This network consists of six attractors, four of which correspond to erythrocytes, megakaryocytes, monocytes and granulocytes. To realise the conversion to granulocytes (A 5 in Fig. 3 (a)) from any initial state, TTC needs to perturb C/EBPα, PU.1, together with one of the nodes in {cJun, EgrNab, Gfi1}. It has been verified that coordinated overexpression of C/EBPα and PU.1 is required for the convergence to GM lineage (granulocytes and monocytes) [17]. One more control node in {cJun, EgrNab, Gfi1} helps to further distinguish granulocytes from monocytes. Efficiency. The last two columns of Table 1 summarise the execution time for computing temporary target control for all the attractors of the networks. We can see that our method is more efficient than SMC for most of the cases. SMC failed to finish the computation for three networks (hematopoiesis, HGF, and HIV-1) within five hours. For the hematopoiesis network, SMC failed in the identification of stable motifs, which has been pointed out to be the most time-consuming part of SMC [46]. The reason could be that the number of cycles and/or SCCs in its expanded network is computationally intractable. For the HGF-induced keratinocyte migration network, SMC is blocked in the optimisation of stable motifs due to that this network has 19 stable motifs and most of the stable motifs contain more than 16 nodes. SMC failed to construct the expanded network representation for the HIV-1 network because some of its Boolean functions depend on many parent nodes (K ≥ 10). Detailed discussion on the complexity of SMC can be found in [46]. The efficiency of our method is influenced by not only the network size, but also the number of attractors and the number of required perturbations. The results show that our method is quite efficient and scales well for large networks.

CONCLUSION
In this work, we have developed a temporary target control method for asynchronous Boolean networks to identify a set of nodes, the temporary perturbation of which can drive the network from any initial state to the desired target attractor. We have evaluated our method on various biological networks to demonstrate its efficacy and efficiency. We compared our method with SMC, a promising method to solve the same control problem. SMC explores both structures and Boolean functions of Boolean networks, and is potentially more scalable for large networks. In contrast, our method is essentially based on the dynamics of the networks, and it will suffer the state space explosion problem for networks of several hundreds of nodes. We believe that these two methods complement each other well. In the near future, we aim to find a way to combine the strengths of both methods by simultaneously exploring network structure and dynamics to achieve more efficient computational methods for the control of large biological networks.