COMPUTATIONAL METHODS FOR ASYNCHRONOUS BASINS

. For a Boolean network we consider asynchronous updates and deﬁne the exclusive asynchronous basin of attraction for any steady state or cyclic attractor. An algorithm based on commutative algebra is presented to compute the exclusive basin. Finally its use for targeting desirable attractors by selective intervention on network nodes is illustrated with two examples, one cell signalling network and one sensor network measuring human mobility.


1.
Introduction. There are three parts to network modeling in the life sciences: 1) build a dynamical model of interacting molecules or other agents using data, statistical learning, expertise, and literature, 2) classify the realistic attractors or steady state possibilities as uninformative, useful, or undesirable, according to the meaning in the model of the node values of the attractor, and 3) seek settings or interventions on accessible nodes to guarantee termination in a useful attractor or prevent termination in an undesirable configuration. This at least is one scenario for applied network modeling. This paper gives an example of a 15 node dynamic network of human mobility variables including cognitive impairment to illustrate step 1). The network is based on our particular study of household remote sensor data of motion variables related to human health, a type of data that is sometimes called "life kinetics" data. A detailed overview is given of the statistical methods that resulted in the model, but the complete analysis that would allow the reader to duplicate each step is not provided here. This model is given as a realistic computational example for the techniques of this paper and not as a conclusive model for the medical or biological sciences. We also give computational methods for problems 2) and 3) that are applied to both a cell signalling example and the life kinetics model.
Networks in biology and life sciences are modeled as dynamical systems in order to understand interactions, causal relationships, and evolution. Boolean or on/off models are the most tractable mathematically and have the most intuitive interpretation. Boolean models continue to be used widely ( [15], [22], [25], [26], [29], [38]). Cancer signalling is a common application, and a growing application area is aging and Alzheimer's disease ( [9], [18], [28]), which is related to the life kinetics example 3392 IAN H. DINWOODIE of Section 4.2. Methods of commutative algebra related to the ones in this paper have been used previously by other authors (for example [19], [20], [24], [40]).
To include variability in the timing of the interactions in the network, random perturbations of the dynamics are of interest ( [6], [30], [36], [42], [45]). The randomized dynamics are a special case of "probabilistic Boolean networks" [39].
The randomization techniques are sometimes called "asynchronous" dynamics to distinguish them from the plain deterministic dynamics where the coordinate maps are simultaneously or synchronously applied to the present state, mapping all coordinates or nodes together to the next state. The asynchronous systems are Markovian, and bring in probabilistic notions of basin of attraction that extend the traditional definition. Here we will focus on one particular scheme, known sometimes as "general asynchronous", in which a node j is selected uniformly at random for update using its transition or update function f j . This paper defines precisely the exclusive asynchronous basin in Section 2 as the set of points that can never leave a deterministic basin under asynchronous updates. Theorem 2.1 gives useful properties following from the definition, especially one that no point in the exclusive basin can reach any other attractor. In Section 3 we describe algorithms based on commutative algebra that compute and analyze the exclusive asynchronous basin. In Section 4 we apply the algorithms to examples from the life sciences of 13 and 15 nodes which demonstrate the effectiveness of the computational tools. The algebraic computations were done in Singular [8], but other algebra software is also suitable, such as Cocoa [1] and Macaulay2 [14].
2. Exclusive asynchronous basin. Consider a state space {0, 1} d , a d-fold product of logical or on/off or high/low symbols 1 and 0. These will be states or configurations on a network with d nodes. Let F = (f 1 , . . . , f d ) be a transition map or transition function or update function on {0, 1} d , where f j : {0, 1} d → {0, 1} is the rule for node j. This map is deterministic, it is the simplified algebraic or logical model of interactions from one time step to the next and F is often called the synchronous update.
Two common choices for randomizing the timing of the updates in F are "random node" (recommended by [36] and [45] where it is called "general asynchronous" or GA), and random order (implemented in [6]). In this work we will only deal with general asynchronous, which we will call simply asynchronous.
Asynchronous updates proceed by choosing a node or coordinate j ∈ {1, 2, . . . , d} with the uniform probability distribution, then changing current state The resulting process x 0 , X 1 , X 2 , X 3 , . . . is a Markov chain in state space {0, 1} d , starting at nonrandom initial configuration x 0 and with random variables X n for successive random states (here we follow the convention that capital letters denote random variables). The probability that X n is in a set S will be written: In theory the process has a probability transition matrix of size 2 d × 2 d but usually it is not practical to work with such matrices directly. Our examples have 13 and 15 nodes respectively, giving cardinalities of 8192 and 32768 for the state space size 2 d . We hope and expect that the algebraic methods in this paper can be applied on examples of up to 50 nodes in the future, as symbolic commutative algebra software gets more powerful.
where F n is the n-fold composition of F . The resulting sets, as x varies in the state space {0, 1} d , are the limiting sets or attractor sets of the system. These are fixed points or cycles, and every state enters one of them and stays there under deterministic dynamics F .
Let B A be the usual deterministic basin of attraction for an attractor A, the points that eventually hit A: Define the single coordinate updates on a state x = (x 1 , . . . , x d ) by the inverse of F j , and B c the complement of a set B, define the exclusive asynchronous basin of the attractor A by Probabilistically, this can be described as In words, the points in B ex,A are the states that can never leave the basin of attraction B A over time under perturbations in the form of asynchronous updates. The sets B i are the initial states x 0 that keep the initial segment x 0 , X 1 , X 2 , . . . , X i in B A , and (2) gets them recursively by removing the points from B i−1 that can move out of B i−1 in one initial step. At the end of this section we discuss the definition further. Note that B 0 ⊃ B 1 ⊃ B 2 · · · and the domain {0, 1} d is finite so the limit B ex,A exists, while possibly empty, and will be a fixed point for the set map T : Our definition of B ex,A starts with the deterministic basin B A , and this gives useful properties on keeping the asynchronous dynamics out of other attractors.
Theorem 2.1. If B ex,A is the exclusive asynchronous basin for attractor A then 1) asynchronous or synchronous dynamics starting at any point in B ex,A cannot reach any point of any attractor other than Proof. For 1), let x ∈ B ex,A ⊂ B A . Then F n (x) ∈ B A , n ≥ 0, and x can never reach any other attractor under deterministic dynamics F . For the asynchronous dynamics, (5) implies ) and x ∈ B ex,A means F j (x) ∈ B ex,A , j = 1, ..., d, in other words P x (X 1 ∈ B ex,A ) = 1. Thus x stays in B ex,A under single-node updates, and B ex,A is disjoint from any other attractor since B ex,A ⊂ B A . For 2), a fixed point for F is also fixed for the asynchronous updates (that is, if A = {x} and F (x) = x then F j (x) = x as well, hence P x (X 1 = x) = 1) and will never leave B A . For 3), which means we can stop the iteration as soon as it repeats, observe that if B i = B i−1 then B i = B i+1 = B i+2 = · · · and the limit is B i .
Note that B ex,A may be empty if A is a cycle. A simple example is useful to understand the definition and its consequences. Consider the state space diagram in Figure 1 for dynamics given in logical form by , and this is also the exclusive basin for the fixed point 00. Now the cycle attractor 01 → 10 → 11 → 01 has empty exclusive basin, because any of its states can transition to 00 with asynchronous updates. For example, 11 → 01 → 00 is possible by updating first coordinate 1 with map f 1 then coordinate 2 with map f 2 . So 11 can reach 00, and 00 is not in the basin of the cycle attractor as it is a fixed point, and 01, 10 are similar. Thus there are no states that cannot leave the basin {01, 10, 11}, and B ex,{01,10,11} = ∅.
This example illustrates also why the iteration (2) should start at B 0 = B A in the presence of cyclic attractors, rather than starting with the inclusive basin B in,A . Recall that the inclusive basin B in,A is the set of points that have positive probability of reaching the attractor A with asynchronous updates. In the above example, the inclusive basin for the fixed point 00 contains all states, and then the iteration (2) would give the entire state space if B 0 = B in,A , including the cyclic attractor, as the exclusive basin B ex,A . Thus, if we want to exclude all elements of other attractors from B ex,A we should start with B 0 = B A and not with B in,A at (2).
Starting with an initial state x 0 in B ex,A as defined in (3) may not guarantee reaching the attractor A under asynchronous dynamics (but it does using the original deterministic map F ). A simple example is a network in d = 2 where states 00, 11 map to 10 under F , 10 maps to 01, and 01 is fixed. Then the basin B 01 for 01 is all four states, as is the exclusive basin B ex,01 . But the state 10 cannot reach 01 under asynchronous updates, it can only stay within {00, 11, 10} with single node updates.
There are a variety of ways to define basins of interest with randomized versions of deterministic dynamics F , as explained and illustrated in [6] and [37]. Randomized versions can be either totally random or partially random, depending on how much of the network is randomized. Then updates can be randomized by choosing an order randomly or a particular node or subset to update. The totally randomized method with individual random node updates is recommended in [36] (and called general asynchronous or GA dynamics) for its realism combined with computational efficiency, and this is what we use.
The definition (3) for the exclusive asynchronous basin is slightly different than the notion in [36] for general asynchronous dynamics. The definition used in [36] is the collection of states in the inclusive basin of A that are not in any other inclusive basin. The two definitions work out the same on the main 6-node example of their Table 1. They are also the same in the Figure 1 example, where B ex,00 = {00} and B ex,{01,10,11} = ∅. But in the previous example (in d = 2 where states 00, 11 map to 10 under F , 10 maps to 01, and 01 is fixed), the two notions are different: definition (3) for B ex,01 gives the entire state space {00, 01, 10, 11} since there is no other attractor basin that the Markov chain could enter, whereas the notion in [36] would make the exclusive basin equal to {01}, since it is the only point in the inclusive basin. This example supports the value of (3) if synchronous dynamics are considered primary and the goal is to prevent the dynamics from entering unwanted attractors under perturbations. The definition in [36] is studied in [11], the algebraic methods work there as well.
The notion of exclusive basin made precise in (3) is computationally efficient, it has a clear probabilistic interpretation, and it has practical value in applications like Example 4.2. It gives a robust basin for the purpose of targeting useful attractors.
3. Algebraic computation. There are three computational problems when trying to target a given attractor under asynchronous dynamics with node interventions: 1) computing the basin of attraction B A , 2) computing the exclusive asynchronous basin of attraction B ex,A ⊂ B A , and 3) finding cylinder sets inside B ex,A . A cylinder set C is a set of the form C = {x : x j1 = a 1 , . . . , x jc = a c } for particular 0-1 values of a 1 , . . . , a c . In other words, certain coordinates are assigned particular values and the others are free. The fixed coordinates j 1 , . . . , j c are the targets for intervention -if C ⊂ B ex,A , then initializing the values of coordinates j 1 , . . . , j c at a 1 , . . . , a c for any starting state x will guarantee convergence to the attractor A under deterministic dynamics and will prevent entering any other attractor under asynchronous dynamics, by Theorem 2.1.
A convenient and effective way to do the three computations uses computational commutative algebra ( [7], [23], [34]) and its implementation in packages like [8]. Although the notion of exclusive asynchronous basin involves probability, the algebraic calculations are exact, unlike inferences about the basin from data generated by computer simulation of the asynchronous process. This technology works with the polynomials that vanish on a set of points (the ideal of the set) rather than directly with the points themselves (the variety). The ideal-variety correspondence [7] allows one to move between the two and find geometric properties of the set of points from algebraic properties of the ideal. To use commutative algebra for dynamics, one must write the dynamics as polynomial functions (see [24] and [40] for background). For example, the dynamics in Figure 1 would be written f 1 (s 1 , s 2 ) = s 1 + s 2 − 2s 1 · s 2 , f 2 (s 1 , s 2 ) = s 1 where the s 1 , s 2 are indeterminates for a polynomial ring C[s 1 , s 2 ], technically not the same as state variables x 1 , x 2 which are merely symbols for unspecified values of 0 or 1.
The remaining parts of this section are about doing the computations described in Section 2 using commutative rings in algebra software and will be of interest primarily to people implementing the methods in software.
Define the ideal I 1 by and define recursively a sequence of ideals I 2 , I 3 , I 4 , . . . by Stop the iteration when dim R/(I i + I 01 ) repeats in order to get the polynomials I B A that vanish on the basin of attraction B A . This solves the first computational problem of getting B A . This algorithm is in [13].

3.2.
Computing B ex,A . To get the exclusive asynchronous subset B ex,A , we remove points algebraically from B A following (2). Define ideals in C[s, t] for single coordinate updates of asynchronous dynamics by Also we will use the notation I s→t for substitution of indeterminates t for s in an ideal I ⊂ C[s].
The way to compute the ideal of B ex,A is In words, we get J 1 by looking first at the complement of B A with the colon ideal R : I s→t 0 in next-step variables t, and by adding I Fj we add equations that constrain the starting variables s to hit the complement of B A with the coordinate update f j . Adding the ideal I 01 is for technical reasons related to variable elimination and extending solutions (see the Extension Theorem in [7]). Then I 1 is the complement of J 1 , the points in s variables that do not exit the basin. Stop the iteration when the ideals I i repeat, using property 3) of Theorem 2.1, which can be detected by measuring their size with dim R/(I i + I 01 ) the vector-space dimension counting the number of points in the 0-dimensional ideal. If the size of I i stabilizes at i then 3.3. Finding control variables. One purpose of a signalling network model is to understand causal relationships among the nodes that can be used to target certain desirable limiting configurations. Examples are in cancer therapy ( [26], [36], [35]) and Alzheimer's disease ( [9], [18]). The effect of node settings is a purely mathematical one in this section. As discussed further in Example 4.2, the causal effect of network node interventions in the larger scheme of gathering data and fitting statistical models to biological data is more complicated.
Ways to find control variables, or nodes for intervention, are developed in [3]. One way is to find the prime decomposition (see [7]) of the ideal I B ex,A which is the algebraic way to find simple cylindrical subsets of the variety B ex,A . To summarize Theorem 2.2 of that work, we do three things to I B ex,A : 1) reduce or divide a standard basis for I B ex,A by I 01 to getĪ, 2) compute the radical ofĪ, say

√Ī
, 3) compute the primary decomposition of √Ī and look for components using variables that can be controlled.
The idea behind why this method works is as follows. Suppose the cylinder The ideal for the irreducible set {x ∈ C d : x j1 = a 1 , . . . , x jc = a c } is prime, say P , and the ideal for {0, 1} d is I 01 defined previously in Section 3.1. So we have inclusions using that the ideal of the intersection of these sets (varieties more precisely) is the sum of their ideals. Then reducing the generators for I B ex,A to getĪ leads toĪ ⊂ P and since P is prime we get √Ī ⊂ P . Then the prime ideal P = s j1 −a 1 , . . . , s jc −a c in the ring C[s] will show up in the prime or minimal decomposition ∩ j P j of the radical ideal

√Ī
. The details are in Theorem 2.2 of [3]. For the simple example in Figure 1, the three calculations would go as follows. Starting with basin B = {00}, its ideal I B = s 1 , s 2 , the polynomials whose roots define x 1 = 0, x 2 = 0. Then I ex,B = I B since 00 is a fixed point for both deterministic and asynchronous dynamics, and I ex,B = s 1 , s 2 . Finally, reducing I ex,B by I 01 does not change anything because s 1 , s 2 are not divisible by s 2 j − s j . Finally, the prime decomposition is nearly trivial, since the ideal s 1 , s 2 is already prime and it corresponds to the cylinder x 1 = 0, x 2 = 0, just a single point.
A more illustrative example of the algebra is one where a hypothetical exclusive basin B ex,A = {000, 001, 010} is a union of two cylinder subsets 00* and 0*0, with intersection the single point 000. This set has ideal given by I B ex,A = s 2 3 − s 3 , s 2 s 3 , s 2 2 + s 2 s 3 − s 2 , s 1 . Then reducing the generating elements of I B ex,A by I 01 and computing the radical gives ideal s 1 , s 2 s 3 , which has prime decomposition P 1 = s 1 , s 3 , P 2 = s 1 , s 2 , corresponding to 0*0 and 00* respectively. 4. Examples. Our two applications for illustrating the methods come from the life sciences. The first cell signalling example starts with dynamics from published sources, whereas the second involves a first step of statistical modeling using modern learning techniques described for example in [16]. We used the software BoolNet [31] to draw wiring diagrams, where green (light) edges mean activation and red (dark) edges mean inhibition, and to draw the state space diagram of Figure 1. We also used it to check calculations on attractors. BoolNet uses a concise and intuitive logical notation for coding functions: & codes for "and", | codes for "or" and ! codes for "not." In Figure 2, the first line "NOS, Ca2" means that the first node is called NOS, and its update has a value given by the current value of the calcium node Ca2. Then the function for coordinate 1 is f 1 (x) = x 11 , since Ca2 is node 11. For node 11, we have f 11 (x) = x 9 · (1 − x 10 ).

4.1.
Cell signalling network. The first example is a signalling network from Table 1 of [36]. The network models the agents (enzymes, proteins, and small molecules) that are involved in signalling stomatal closure in a cell. The closure is a function of the calcium node Ca2 and is not actually present in the coding of Figure 2. The wiring diagram is in Figure 3. The meanings of the names of the nodes in Figure 2 are given in Figure 1 of [37].  The Boolean network has one fixed point 0000000000010 in which node 11 Ca2 is off (implying no stomatal closure) and two cycles of length 4. The fixed point has deterministic basin of size 108, and an exclusive asynchronous basin of size 56, composed of two cylinders x 1 = x 2 = x 3 = x 6 = x 7 = x 8 = x 9 = x 11 = 0 and x 1 = x 2 = x 4 = x 5 = x 7 = x 8 = x 9 = x 11 = 0 each of size 32, and overlapping in 8 points. The computer time to get the synchronous basin was .060 seconds, whereas the asynchronous basin required 74.420 seconds (Intel x86-64 processor i7-3770 at 3.40GHz) and memory use was under 2 Gbytes always.
Any of the 8136 = 2 13 − 56 points outside of its exclusive basin has a positive probability of entering the deterministic basins of the two cycles under asynchronous dynamics. Both cycles show oscillatory behavior in the Ca2 node, and oscillation is of special interest to biologists [37].
One cycle shows oscillation in all 13 nodes: on) will move to this attractor. Note that the values on these nodes will change over time, as they must since the attractor shows oscillation on all nodes. The exclusive asynchronous basin is empty, meaning that changing the dynamics with random node updates will allow the system to move into other attractor basins with positive probability. A second cycle is fixed in the CIS node 9: 1100001111001 0111000110010 0011110010110 1000111011101 Here the deterministic dynamics give 16 subcylinders of the basin, with some two-node sets like x 7 = 1, x 11 = 1. Again, the exclusive asynchronous basin is empty.
The empty asynchronous basins show the network is fragile in the sense that randomized dynamics can change the limiting behaviour significantly with positive probability, a property also noted in [6].

4.2.
Life kinetics network. Our second example involves longitudinal data from in-home motion sensor recordings. This data was generated by the Intelligent Systems for Assessing Aging Changes (ISAAC) study [21], an aging study conducted by the Oregon Center for Aging & Technology (ORCATECH) at OHSU. The particular data set involved 153 people, 197 homes, and spanned the time period May 2007 to February 2014. Our data file recorded summary data for each day, for each residence. The variables to be measured and computed were determined by ORCATECH based on expertise and current research including [2], [17], [32], where variables like the number of transitions from one room to another, or the time spent out of house, were shown to be relevant. It is known that variability of quantities like walking speed is important, so we also had measures of dispersion like the standard deviation of walking speed over each day, and the coefficient of variation of walking speed (standard deviation divided by mean) on each day to normalize for subject differences. The data is observational and not controlled experimental data, as variables like "walking speed" were not applied in a randomized way to subjects. So dynamical patterns in the data found with statistical learning methods are comparable to correlations, and causal effects can only be determined with further study.
A variable of interest is the binary variable MCI (mci) which codes for mild cognitive impairment. While MCI is not the same as Alzheimer's disease, the two are related [33]. One goal is to understand the predictive roles of other behavioural variables. In particular, it would be useful to know if interventions may be possible on some variables to reduce the risk of future MCI.
Let us say why a Boolean network model may be useful for this data. A Bayes network can be used to model dependencies in multivariate data, but the model is static and it does not help to understand the predictive role of variables related to MCI months in advance. In the work [5], the predictive role of walking speed years in advance on MCI was established using change point methods, statistical methods which are related to on/off signalling. But [4] emphasizes that monotherapeutics (therapies focusing on one variable) have been ineffective, and presents a new therapeutic system involving over twenty actions designed to modify behaviour variables like sleep and exercise, variables which also appear in some form here. The Boolean network method can give an enhanced understanding on how combinations of interventions have signalling effects over time. It is one of the simplest ways to understand multivariate signalling over time, our main goal. The methods of survival analysis are also applicable, but they are not able to capture cyclic behaviour of MCI in some individuals. It is possible to extend the techniques of this paper to predictive variables with more levels than two [10]. Finally, many statistical methods are useful for understanding multivariate longitudinal data and the Boolean model that we use here is just one tool.
A large scale evolving study like ISAAC will have challenging statistical issues. About 16% of the entries are missing in our working data file of 43 Megabytes, and outliers are present due to imperfect sensor readings. The reasons why we consider our fitted model in Figure 4 only a first step in understanding the dynamics of the data are: 1) our missing data procedure is crude, 2) tree classification methods are known to be sensitive to small changes in the data, and 3) the best time lag for regression is not established. Now we describe the statistical procedure that gave the estimated dynamics of Figure 4. First, the data was revised by removing the few records with missing MCI value or missing subject identifying number. Then other missing values were imputed with the function na.roughfix of the randomForest package [27]. This is a very crude tool for missing data, but the more sophisticated function rfImpute would not handle so much data with only 32 Gigabytes of computer memory. A more refined analysis in the future should improve this step. Then, after imputing to get complete data, we used the variable importance feature of the randomForest package to rank the variables for inclusion in the final model. This method essentially looks at how much worse prediction with cross validation becomes when the variable is removed from the model (see p. 593 of [16] p. 593 for a detailed description). From approximately 30 variables, we chose 15 including the MCI value. Many of the variables are continuous, so it is necessary to discretize into high/low or on/off values for a Boolean network. Rather than using quantiles or other methods that are unrelated to signalling or prediction, we used a new predictive method that is based on random forests. For this, we get 500 classification trees from the randomForest command predicting MCI on one day from the predictor values on the previous day. Then we examined each tree for the value of the first split on each predictor variable. The cutoff or threshold for on/off was the mean of the 500 values of the first split (each of the 500 trees split every predictor variable in its partition of the domain for building a function y = f (x 1 , . . . , x 14 ), where y is a binary code for MCI). The variable definitions are in Table 1.
Then finally tree classification was implemented on the resulting binary variables. For the tree classification, we used the R package rpart [41]. This procedure is described in [16] and is different than the interpolation method of [43]. Converting the trees to algebra was done as in [12]. The wiring diagram in Figure 5 is a graphical representation of Figure 4, where green arrows indicate an effect of switching on or up (activation) and red arrows represent switching off or down (inhibition). The wiring diagrams are slightly ambiguous because of interaction effects. The arrow shows the direction of the effect in time. The Boolean dynamical system from statistical learning is in Figure 4. This model is one that is consistent with many observed interactions among variables related to aging, but should not be taken as a final analysis. It has 52 attractors, including 7 steady states and 45 two element cycles. Six of the steady states have MCI off, so these are interesting and useful as targets, as are the 33 cycles in which MCI is off.
One of the interesting steady states is 01101 10000 00000, with coordinate labels in the order of Figure 4. Thus MCI has value 0, which is a desirable configuration. The deterministic basin has size 3456, whereas the exclusive asynchronous basin has size 960. Such a significant reduction in size is expected, as this was observed in   Figure 5. Wiring Diagram for Life Kinetics and Cognitive Impairment [36] and [37]. The computational time in Singular [8] (version 3-1-6) is considerably higher for the asynchronous basin. On a desktop computer (32 Gbytes of memory, Intel x86-64 processor i7-3770 at 3.40GHz), the synchronous basin B A took 2.320 seconds of computing time, whereas the exclusive asynchronous basin B ex,A took 1090.100 seconds. The computer memory use for the asynchronous basin calculation was less than 2 Gbytes throughout the computation. The basin B ex,A is the union of five intersecting cylinder subsets: C 1 = {x : x 2 = 1, x 3 = 1, x 8 = 0, x 9 = 0, x 10 = 0, x 11 = 0} C 2 = {x : x 2 = 1, x 3 = 1, x 6 = 0, x 9 = 0, x 10 = 0, x 11 = 0} C 3 = {x : x 1 = 0, x 2 = 1, x 3 = 1, x 5 = 1, x 6 = 1, x 8 = 0, x 9 = 0, x 11 = 0} C 4 = {x : x 2 = 1, x 3 = 1, x 5 = 1, x 6 = 1, x 8 = 1, x 9 = 0, x 10 = 1, x 11 = 0} C 5 = {x : x 1 = 0, x 2 = 1, x 3 = 1, x 5 = 1, x 6 = 1, x 9 = 0, x 10 = 1, x 11 = 0} Cylinders C 3 and C 5 use node 1 (MCI) so these are of no practical use for interventions. However C 1 may be useful: it says if we initialize numtrans at high, numfir at high, wscv low, ttib low, wssigma low, and timeasleep low, then the system will lead to this steady state deterministically and will not lead to any attractor with MCI on under asynchronous dynamics. The six nodes for this intervention are known to be related to MCI individually, but this new analysis gives the precise combination of settings that will turn MCI off over time. The results are consistent with current knowledge. One can study all steady states to seek the most practical nodes for intervention. All of the 33 cycle attractors have empty asynchronous basins.

Conclusions.
We have defined an exclusive asynchronous basin in order to eliminate initial states from a deterministic basin of attraction that could move to unwanted attractors under uncertainty in the timing of the node updates. Variations on the definition are possible by changing the initial set B 0 in the iteration. Algebraic computational methods are given. Altogether, for Boolean networks of up to around 15 nodes, the methods give a way to find node interventions that will set an initial state moving towards a desired limiting state even under random perturbations of a model.
The examples show that analyzing asynchronous dynamics is more demanding computationally than a deterministic analysis. Still, existing symbolic computational algebra software is sufficient for precise geometric calculations needed to understand network control variables on real examples. While other computational methods exist for finding deterministic and asynchronous basins, the algebraic geometry is especially good for finding combinations of node settings that target a given attractor in the presence of randomized dynamics. Calculations on a 15 node network use a space of polynomials with 30 indeterminates and this works well in practice. An interface between a package like BoolNet [31] directly to Singular [8] would be convenient and not difficult.