Analysis of a dynamic model of guard cell signaling reveals the stability of signal propagation

Background Analyzing the long-term behaviors (attractors) of dynamic models of biological systems can provide valuable insight into biological phenotypes and their stability. In this paper we identify the allowed long-term behaviors of a multi-level, 70-node dynamic model of the stomatal opening process in plants. Results We start by reducing the model’s huge state space. We first reduce unregulated nodes and simple mediator nodes, then simplify the regulatory functions of selected nodes while keeping the model consistent with experimental observations. We perform attractor analysis on the resulting 32-node reduced model by two methods: 1. converting it into a Boolean model, then applying two attractor-finding algorithms; 2. theoretical analysis of the regulatory functions. We further demonstrate the robustness of signal propagation by showing that a large percentage of single-node knockouts does not affect the stomatal opening level. Conclusions Combining both methods with analysis of perturbation scenarios, we conclude that all nodes except two in the reduced model have a single attractor; and only two nodes can admit oscillations. The multistability or oscillations of these four nodes do not affect the stomatal opening level in any situation. This conclusion applies to the original model as well in all the biologically meaningful cases. In addition, the stomatal opening level is resilient against single-node knockouts. Thus, we conclude that the complex structure of this signal transduction network provides multiple information propagation pathways while not allowing extensive multistability or oscillations, resulting in robust signal propagation. Our innovative combination of methods offers a promising way to analyze multi-level models. Electronic supplementary material The online version of this article (doi:10.1186/s12918-016-0327-7) contains supplementary material, which is available to authorized users.


Background
Modeling offers a comprehensive way to understand biological processes by integrating the components involved in them and the interactions between components. Models can recapitulate and explain the emergent outcome(s) of the process [1,2]. Representing cellular processes that involve many proteins and small molecules by a signal transduction network can reveal indirect relationships between components and provide new insight [3][4][5]. Such network usually consists of nodes representing biological entities, and edges representing interactions. Once a network has been constructed, dynamic modeling, where each node in the network is associated with a variable representing its abundance or activity, can further describe the behavior of the network. Dynamic models can have continuous variables whose change is described by differential equations [6], discrete variables described by discrete (logical) regulatory functions [7,8], or a combination of continuous and discrete variables [9]. The major advantage of discrete dynamic and continuous-discrete hybrid models is that they use many fewer parameters than continuous models and thus need less parameter estimation [10][11][12]. Modeling allows one to analyze the biological system represented by the network in silico, when performing the relevant experiment is infeasible. It also helps identify general principles of biological systems [13,14].
The biological process of stomatal opening in plants is a good example of a complex system wherein modeling leads to significant gain in understanding [15,16]. Stomata are pores on leaf surfaces that allow the plant to exchange carbon dioxide (CO 2 ) and oxygen with the atmosphere. Stomata are formed by two guard cells that can change shape: swelling of guard cells leads to stomatal opening; their shrinking leads to stomatal closure. The shape of each guard cell is directly controlled by water flow through the membrane, which is in turn controlled by ion flow. Different signals can affect the guard cell, changing its ion concentration in direct and indirect ways, resulting in stomatal opening or closure [17][18][19]. These signals include light of different wavelengths, CO 2 concentration in the air, and plant hormones like abscisic acid (ABA). The regulation of stomatal opening is essential to plants, as it controls vital activities like the uptake of CO 2 for photosynthesis, and the unavoidable water loss through evaporation [20]. Through extensive experimentation over several decades, more than 70 proteins and small molecules have been identified to participate in this process.
Sun et al. [15] recently constructed a signal transduction network based on conclusions from more than 85 articles in the literature, describing how more than 70 nodes (proteins, small molecules, ions) interact with each other in the stomatal opening process. The network, reproduced as Fig. 1 [15], includes four source nodes that correspond to the signals red light, blue light, CO 2 , and ABA. The more than 150 edges are directed and signed, with arrowheads indicating activation and terminal black circles indicating inhibition.
Translating this network into a dynamic model, Sun et al. characterized each node with a discrete variable describing its activity and with a discrete (logical) regulatory function describing its regulation. Twenty-one out of the 70 nodes in the model are multi-level, the rest are Boolean (binary). The levels reflect relative and qualitative information: a level of 2 is a higher level than 1, but should not be interpreted as twice as high. A few Fig. 1 The signal transduction network responsible for stomatal opening, as reconstructed by Sun et al. [15]. The color of a node marks which signal regulates this node. Red nodes are regulated solely by red light. Blue nodes are regulated solely by blue light. Yellow nodes are regulated solely by ABA. Grey nodes are regulated by CO 2 . Purple nodes are regulated by both blue and red light. Green nodes are regulated by blue (and potentially, red) light and ABA. White nodes are source nodes not regulated by any of the four signals. To improve visualization, certain pairs of edges with the same starting or end nodes overlap. Nodes with multiple levels in the dynamic model are represented by red shadows; the others are Boolean. The full names of the network components denoted by abbreviated node names are given in Table 1. This figure and part of its caption is reproduced from Sun Z, Jin X, Albert  discrete values are not integers; e.g. stomatal opening is a weighted sum with non-integer weights. The dynamic model has~10 31 states. The logical regulatory functions, describing each node's future state based on the states of the node's regulators, use a combination of Boolean logic operators (And, Or, Not), algebraic operations, and input-output tables. For example, the regulatory function of PRSL1 is: Here for simplicity the node states are denoted by the node names; the asterisk in "PRSL1*" indicates that this will be the next state of the PRSL1. The "Or" Boolean operator expresses that either of the blue light receptors, i.e. the phot1 complex or phot2, can independently activate PRSL1.
The Sun et al. model starts from an initial condition representative of closed stomata. Then a combination of the four input signals is applied. Red light, blue light, and ABA are represented as binary variables, and external CO 2 is represented with three states: 0 (CO 2 free air), 1 (ambient CO 2 ) and 2 (high CO 2 ). The system's response is simulated through repetitive re-evaluation of each node's state until a stable value of stomatal opening is observed. The model successfully captures stomatal opening in response to combinations of the signals. It also successfully reproduces stomatal opening under most of the experimentally studied perturbation scenarios (i.e. genetic knockouts or external supply of components). In total, the model is consistent with 63 out of 66 experimental observations collected by Sun et al. [15]. The model predicts the outcome of a large number of scenarios that have not been explored experimentally so far. It also revealed a gap of knowledge regarding the cross-talk of red light and ABA signaling, and filled it with a newly predicted interaction.
Although the Sun et al. model recapitulates existing knowledge and offers new predictions, the model's full dynamic repertoire could not be characterized due to its large state space. Instead, Sun et al. focused on tracking the output node, stomatal opening, and a few selected internal nodes, in time. In this paper we apply multiple methods to analyze the model and aim to fully map all its potential long-term behaviors, or in other words, attractors.

Attractors of a dynamical system
An attractor is a set of states from which only states in the same set can be reached. Attractors that consist of a single state are called stable steady states or fixed points; attractors that contain multiple states are called complex attractors or oscillations [10]. In biological networks, attractors often have significant biological meaning. In a cell signaling network, attractors correspond to cell types, cell fates or behaviors [21]. For example, one attractor can represent a healthy differentiated cell, while another attractor can represent an abnormally motile cancer cell [22].

Update scheme of a discrete time model
In the Sun et al. model, as in most discrete dynamic models, time is an implicit variable. As there is very little information about the kinetics of the nodes in the stomatal opening network, the model incorporates an element of stochasticity in timing. The timing does not affect a system's fixed point attractors, but it can change the complex attractors and the possibility of reaching a given attractor from a given initial state [10]. In the Sun et al. model, a random-order asynchronous update is used. Specifically, at each time step, a random order of nodes (excluding the four input nodes and the output node stomatal opening) is generated, and each node's state is reevaluated in this order; stomatal opening is always updated last. In the next time step a different order is generated randomly. In this paper, we use a different type of stochastic update, called general asynchronous update, wherein a randomly selected node is updated at each time step. This is required by the network reduction method we use. Although this theoretically could cause a difference in complex attractors, we will show that in this specific model the two update methods yield the same attractors.

Network reduction
To reduce the Sun et al. model's state space, we apply a network reduction method developed by Saadatpour et al. [23] that is proven to preserve the attractors of a Boolean model. Two types of nodes can be reduced (eliminated or merged): source nodes with no incoming edges, and simple mediator nodes that have one incoming and/or one outgoing edge. In the reduction, the source node's state is directly plugged into the regulatory function of all of its direct successor nodes; then the source node is eliminated. For a simple mediator node with one predecessor (regulator) and one successor (target), its regulator is connected to its target and the mediator node is merged into the regulator. If there is one regulator and several targets of the mediator node, but no direct edges between the regulator and any of the targets, the mediator node is merged into the regulator. Conversely, if there are several regulators and one target of the mediator node, but no direct edges among any of the regulators and the target, the mediator node is merged into its target. Although this method is not proven in the multi-level case, we conjecture that attractors are also conserved for a multi-level model, and will show from the results that in the Sun et al. model this reduction method preserved all attractors.

Elimination of redundant edges
During the process of creating a discrete dynamic model from biological data, when an influence is weaker than other influences, the modeler may choose to omit this influence or, alternatively, include it a redundant way. The latter choice was made by Sun et al. in four cases, leading to four regulatory functions that contain an input that does not affect the outcome of the regulatory function. One of these is The italicized words "And", "Or" and "Not" are Boolean logic operators; the non-italicized words represent node names. In this regulatory function every node is Boolean (binary). The first clause "NADPH And Atr-bohD/F" and the second "NADPH And AtrbohD/F And CDPK" are connected with an "Or" rule, with the result that the node "CDPK" does not have any influence on the outcome. Therefore, we can prune the edge from CDPK to ROS without changing the model's dynamics. We similarly prune three additional redundant edges.

Converting a multi-level model to Boolean
There are several possibilities to convert a multi-level model to Boolean [24]. The standard method used in the case of logical models of regulatory networks is the Van Ham mapping [25,26]. It preserves the dynamics of the original model if the variables in the original model can be represented by integers and if the original model only allows state transitions in which one node changes its state by one level [26]. The Sun et al. model does not satisfy these criteria. However there still is a conclusion that we can use: All types of conversions maintain the fixed points and the reachability of states (i.e. if there is a sequence of state transitions from state A to state B before conversion, there must be a sequence of state transitions from the corresponding state A' to state B' after the conversion) [26]. So the worst distortion of attractors due to the conversion is the merging of two complex attractors into one. In this light we choose to use an economic mapping of each multi-level node into as many Boolean nodes as necessary for the binary representation of the corresponding integer. We will show that in this specific model, the conversion did not change the attractors. Table 1 summarizes the full names of the network components denoted by abbreviated node names in Fig. 1.

Network reduction
The Sun et al. model has a huge state space of~10 31 states, making its analysis difficult. To obtain a smaller state space, we reduce the size of the network by applying a network reduction technique developed by Saadatpour et al. [23] that is proven to preserve the attractors of Boolean models (see Methods). All source nodes other than the four signals (blue light, red light, CO 2 , and abscisic acid) and all simple mediator nodes are identified and reduced. This process is done iteratively until it cannot be done any more. A total of 7 source nodes (14-3-3 protein-phot1 , PIP2 C , AtNOA1, Nitrate, PP1 cn , mitochondria, and CHL1), and 19 simple mediator nodes (phot1, phot2, NIA1, H + -ATPase, LPL, ATP, acid. of apoplast, PA, ABA receptors, OST1, PRSL1, PIP2 PM , AtrbohD/F, Nitrite, and phot1 complex ) are eliminated. Several of the simple mediator nodes form linear paths (e.g. phot1, OST1) thus their iterative reduction shortens the linear paths in the network. In addition, 16 of the 19 reduced mediators have a regulatory function of the form "B* = A". It is intuitive that reduction of this node type preserves the attractors.
We do not eliminate the four signal nodes because we want to simultaneously explore all the combinations of input signals. We also choose to not reduce the five nodes (K in , K out , K c , Ca 2+ -ATPase, mesophyll cell photosynthesis) whose merging with their sole regulator would result in a self-loop (self-regulation), because such selfloops may be difficult to interpret. Two additional nodes with significant biological meaning to the network (sucrose, stomatal opening), are not reduced either.
Another form of network reduction is the elimination of redundant edges (see Methods). After removal of redundant edges, the node CDPK becomes a sink node, thus it can also be eliminated. The reduction of the above-described nodes and redundant edges simplifies the network from 70 nodes to 42 nodes, with an estimated state space of~10 22 states.

Simplification of regulatory functions
In order to further reduce the state space from~10 22 to a manageable size, we grouped state values so that nodes are represented with fewer states. This grouping was guided by the 66 experimental observations summarized in Sun et al.; we aimed to maintain the reduced model's results consistent with these experimental observations.
For example, in the Sun et al. model [15] the regulatory function of Stomatal Opening is a weighted sum of different ions and sucrose: The weights of the anion contributions to the osmotic potential were chosen based on the literature. Also, the anion contributions must not exceed a proportion of [K + ] v due to charge balance. The anion contributions are [malate 2 The primary contributions come from [K + ] v and sucrose. We grouped the stomatal opening values into 6 groups with different [K + ] v and sucrose values (see Table 2 and Additional file 1).
The  Similarly to the original model, the simplified states represent qualitative, relative categories. For example, a stomatal opening level of 2 is not twice as high as level 1. We choose the simplified stomatal opening values so that there is no state "4", to better reflect an experimentally observed synergistic effect between blue and red light [18,19,27]. Simulation results with the simplified regulatory function are that under monochromatic red light stomatal opening =1; under monochromatic blue light stomatal opening =3; under dual beam the stomatal opening =5, which is larger than the sum "1 + 3". This qualitatively reproduces the experimental observation that under dual beam illumination stomata open to a size much larger than the sum of opening under monochromatic blue or red light.
We find by simulation of the reduced model, using the same initial condition as the Sun et al. model, that the simplification of the stomatal opening regulatory function results in only 3 additional cases of inconsistency with experimental observations out of a total of 66 experimentally studied scenarios. Additional file 2 lists all experimental observations and compares them to the relevant simulation results. Ignoring the contribution of malate 2− , NO 3 − , and RIC7 to stomatal opening each causes one additional discrepancy; ignoring Cl − does not cause any additional discrepancy. Ignoring these nodes trades a decrease in accuracy for a significant increase in simplicity.
The simplification of the stomatal opening regulatory function eliminates the effect of vacuolar anions and of RIC7 − ] a , ROP2, RIC7, ABC, and PEPC. The only edge from these nodes to other nodes is [malate 2− ] a → AnionCh. In section 3 of Additional file 3 we show that eliminating this edge does not change the system's longterm behavior, i.e. attractors. Also, the regulatory function describing the cytosolic K + concentration, [K + ] c , can be simplified without loss, as described in section 3 of Additional file 3. After this simplification we have a network of 32 nodes, 81 edges, indicated on Fig. 2. We will refer to this model as the "reduced model". A list of nodes and their regulatory functions is provided in Additional file 1.
Identifying strongly connected components (SCCs) is important for attractor analysis, as complex dynamic behavior such as oscillations or multi-stability requires feedback loops [7]. There are three SCCs in the network of the reduced model, as marked in Fig. 2. The NO cycle contains three nodes and three positive edges. The C i SCC contains three nodes, which form two negative feedback loops. The Ion SCC is the most complex, containing 13 nodes and 26 edges, 7 of which are negative.
Next we perform attractor analysis using two methods: 1. by converting the reduced model to Boolean and applying two analysis tools; 2. by analyzing the regulatory functions theoretically. The former method finds all stable steady states and candidate oscillations; the latter confirms the results of the first method and gives insight about perturbation scenarios.

Conversion of nodes from multi-level to Boolean states and attractor analysis
We perform the conversion to Boolean to enable attractor analysis by existing software tools. Zañudo et al. [28] proposed an algorithm to find the attractors of a Boolean network based on the concept of "stable motif", a strongly-connected group of nodes that can stabilize regardless of their inputs. The algorithm finds all stable motifs, which determine the part of the network that stabilizes in an attractor. After a stable motif is found, one can plug in its stabilized state into the network, and obtain a smaller remaining network. After repeating this, eventually the remaining part is either nothing (indicating a fixed point/stable steady state) or a candidate oscillating sub-network. Compared with other software tools [29,30], the major advantage of this algorithm is that it finds all the attractors of Boolean networks with hundreds of nodes [28]. Application of this powerful method requires a Boolean model, so we convert the multi-level model into Boolean first (see Methods). An example of conversion is given in Table 3. More detailed examples of the conversion of the states and regulatory function of specific nodes are given in the Additional file 4. We will refer to the reduced model after conversion to Boolean variables as the "Boolean-  represent the same entity (the same multi-level node) are updated simultaneously. In this way the state transitions of the reduced model will be kept the same in the Boolean-converted reduced model, and therefore the Boolean conversion will not cause additional discrepancies from experimental observations. We apply the stable motif algorithm's implementation, downloaded from http://github.com/jgtz/StableMotifs/ [28], to the Boolean-converted reduced model. The algorithm uses the Boolean regulatory functions of the converted model (given in Additional files 5 and 6) as input. We consider every combination of sustained states of the five signal nodes (blue light, red light, ABA, CO 2 , CO 2 _high). We find two possible stable motifs, corresponding to the self-regulatory node PMV_pos (one of the two Boolean nodes associated with the multi-level node PMV, see Additional files 4 and 5), in conditions where the H + -ATPase complex is inactive. These two stable motifs indicate the bistability of PMV. Under its influence, another node, K out , will also be bistable. The algorithm also indicates that for any signal combination, every node, except [Ca 2+ ] c and Ca 2+ -ATPase, will stabilize in a fixed state. [Ca 2+ ] c has three states, and in the Boolean-converted model it is represented by two nodes, Cac and Cac_high. Cac_high, which represents the higher level of [Ca 2+ ] c , stabilizes at zero in all situations. Cac and Ca 2+ -ATPase may oscillate in conditions where blue light is present and ABA is absent (a total of six cases, two of which allow PMV bistability). Table 4 summarizes key features of the attractors found by the stable motif algorithm for all 24 input combinations. Attractors where Ca 2+ oscillation is not possible are fixed points (stable steady states).
We verified the obtained attractors with GINsim [12], a software suite capable of model construction, simulation, and analysis. GINsim can compute all stable steady states (called stable states in GINsim), or determine complex attractors by mapping the state transitions. The stable steady states found by GINsim are identical to those found by the stable motif algorithm. To verify and further explore the complex attractors, we use the simulation function of GINsim, starting from a state in the complex attractor. The result that the system oscillates between four states, where only the state of Cac and Ca 2 + -ATPase changes, agrees with the findings of the stable motif algorithm. We summarize the GINsim computation/simulation results in Additional file 7. Additional file 8 indicates the Boolean-converted reduced model in SBML-qual format [29], a general format for biological model to be analyzed using various tools including GINsim.
We can also connect the stable motif analysis results to network reduction. We have previously decided to not reduce the four nodes that correspond to input signals. If we do consider a specific input combination when using network reduction, e.g. blue light and red light with normal CO 2 without ABA, we can reduce much more of the network: two of the three SCCs, namely the NO cycle and the C i SCC, will stabilize and can be eliminated. Only the Ion SCC and its sole output stomatal opening remain, indicating that this SCC is not driven solely by the external signals and has the capacity for oscillations or multi-stability. This is consistent with the results found by stable motif analysis, according to which the NO cycle and the C i SCC attain a steady state and the Ion SCC admits a [Ca 2+ ] c -Ca 2+ -ATPase oscillation and PMV bistability. This consistency supports the appropriateness of the network reduction method and of the Boolean conversion.

Theoretical analysis of the reduced model
To gain additional insight into the attractors of the reduced model and their potential changes due to node perturbations, we analyze the reduced model theoretically. Specifically, we aim to answer the question: Can there be other types of oscillation, or can there be additional multi-stability, if a node is knocked out (fixed in the OFF state) or is constitutive active (fixed in the highest state)?
We first test whether the network and regulatory rules allow multi-stability or oscillations. This analysis is based on R. Thomas's conjectures [7]: The presence of a positive (negative) feedback loop -a cycle with an even (odd) number of inhibitory edges -in the network is a necessary but not sufficient condition for the occurrence of multiple steady states (oscillations). The conjectures have been proven in the case of discrete dynamic systems [31][32][33][34]. Since only feedback loops are candidates for potential multi-stability or oscillations, we analyze the regulatory functions of each strongly connected component of the network. For each feedback loop, we identify a sufficient condition for the nodes to stabilize in a specific state. The violation of this condition becomes a further necessary condition of multi-stability or oscillation. Here we describe the main steps and results of the analysis; the detailed analysis is in Additional file 3. The multi-level node shown in the 1st column is mapped into two Boolean nodes, shown in the 2nd and 3rd columns, using the binary representation of the corresponding integer.
The NO cycle is composed of the nodes PLD, ROS, NO, and the three positive edges between them. It does not have any negative edges, so it cannot oscillate. A fixed ABA value is sufficient to stabilize each node of the cycle in a specific state, thus the cycle does not admit multi-stability under any perturbation.
The C i SCC has three nodes, C i , mesophyll cell photosynthesis (MCPS), carbon fixation, and four edges that form two negative feedback loops, one between carbon fixation and C i , and the other between C i and MCPS. Despite the existence of negative feedback, this cycle will stabilize if given a fixed CO 2 value. From this we know that this cycle cannot oscillate or admit multi-stability under any perturbation.
The Ion SCC has 13 nodes. To reduce its complexity we show that the key node [Ca 2+ ] c , which has states 0,1, and 2, cannot enter state 2 in the long term under any perturbation. Since most nodes respond to [Ca 2+ ] c only if [Ca 2+ ] c =2, we can eliminate all edges that depend only on "[Ca 2+ ] c =2", and obtain a simplified Ion SCC, as shown in Fig. 3. The Ca 2+ SCC ([Ca 2+ ] c , Ca 2+ ATPase, PLC, CaR) now becomes a sink SCC. The only negative edge in this sub-network is from Ca 2+ -ATPase to [Ca 2+ ] c . These two nodes are known to oscillate. The positive feedback loop formed by [Ca 2+ ] c , PLC, and CaR will stabilize if given fixed inputs. So there cannot be multi-stability. For the nodes outside of the Ca 2+ feedback loops, we show that the edges from KEV and [K + ] v are redundant in the long term, so there are no feedback loops except the PMV selfloop. PMV is not capable of having oscillations, but can have bistability (as also indicated by the stable motif analysis). The bistability can affect at most one other node, K out , under any perturbation. This means that the bistability has very limited effect on the attractor of the reduced model. Now we can summarize our conclusions and return to the question we sought to answer: there is no oscillation except in the calcium nodes; there is no multi-stability except in the nodes PMV and K out . These statements are true under any perturbation. Moreover, for the calcium oscillation, [Ca 2+ ] c cannot enter the state 2, so the subnetwork between [Ca 2+ ] c and Ca 2+ -ATPase is a negative feedback loop between two Boolean nodes, with the regulatory functions Ca 2+ ATPase * = [Ca 2+ ] c ; [Ca 2+ ] c * = not Ca 2+ ATPase. It results in the simplest type of oscillation, as also found by GINsim simulation. For the PMV bistability, even if the bistability exists, most nodes, especially the output node stomatal opening, still have a unique value. Thus the theoretical analysis, in agreement with the computational analysis, leads to very strong conclusions about the reduced model's dynamic repertoire.
We can also show that the reduction or Boolean conversion did not change the attractors of the Sun et al. model. Although the reduction we used is only proven in the Boolean case, Naldi et al. showed that for multivalued models, removal of non-autoregulated nodes, like in our reduction, preserves crucial dynamical properties [35], including fixed point attractors and the two-node simple oscillation we found. So our reduction is valid in this specific model. To confirm that the Boolean conversion preserved attractors, we note that in the Boolean- The first five columns indicate the input signal combination. The setting CO 2 _high = 1 and CO 2 = 0 is not included because it is not biologically meaningful. The "SO (Bool)" column indicates the state of the Boolean node combination representing stomatal opening. The "SO" column is the state of stomatal opening when converted back to an integer. Note that the stomatal opening level of four is not defined, and no attractors have a stomatal opening level of two. The next column indicates whether Ca 2+ oscillation can possibly happen under the given signal combination. The last column indicates whether bistability of PMV_pos can be observed under this setting. In those cases, two stable steady states with (PMV_pos = 0, K out = 0) and (PMV_pos = 1, K out = 1) can be observed. The rest of the nodes are unaffected by this two-node bistability converted reduced model we found fixed point attractors and a complex attractor in which only two nodes oscillate. Because the only potential change to attractors as a consequence of the conversion is merging of complex attractors [26], it is straightforward that the attractors have been conserved during the conversion, as the twonode oscillation found is the simplest type of complex attractor and cannot be a result of attractor merging. In addition, using general asynchronous update instead of random order asynchronous update does not cause any changes to the attractor, because the update schemes do not affect fixed points or the two-node simple oscillation we found.

Stability of guard cell signal transduction
Our previous results indicate the stability of the system in the sense that all the initial conditions lead to the same attractor except for up to four nodes. We also examine another facet of the system's stability: the robustness of the stomatal opening in response to node perturbations that render them non-functional. We perform a systematic analysis of single-node knockouts of every non-signal node in the reduced model, under all combinations of light, CO 2 and ABA conditions. For each signal combination, we set the perturbed node's initial state and regulatory function to 0, initialize the rest of the nodes in the condition representative of closed stomata, and then simulate the reduced model until it reaches its attractor. In the absence of ABA under each light and CO 2 condition, 60-90 % perturbation scenarios produce the same stomatal opening value as the unperturbed system (Table 5). These results are similar to those reported by Sun et al. for the original model [15] (see Additional file 9). In the presence of ABA 50-90 % perturbation scenarios produce the same stomatal opening value as the unperturbed system, and 4-16 % knockouts lead to a higher stomatal opening value. Perturbations in the ABA = 1 case were not studied by Sun et al., but our simulations of the original model give the same qualitative results as the reduced model. These results indicate the closeness of the perturbed attractor (at least in terms of the stomatal opening value) to the unperturbed attractor in more than 50 % of single node perturbations. They also suggest the resilience of the stomatal opening process against internal failures and perturbations.
Extending the conclusions to the original model We found that in the reduced model there is no oscillation except in the calcium nodes; there is no multistability except in the nodes PMV and K out . Because the reduction we used has been shown to conserve attractors [23,35], we know that our attractor conclusions can be immediately extended to all nodes in the original model except the reduced nodes and stomatal opening. Next we extend the attractor analysis to include the reduced nodes as well.
First we consider the nodes reduced during the first step of network reduction, i.e. non-signal source nodes and simple mediator nodes. These nodes are trivially incapable of having multi-stability and oscillations themselves, so we need only to consider their perturbations. Perturbation of a simple mediator node can always be replaced by a corresponding (set of ) perturbation(s) in the mediator node's direct successor(s), so these perturbations have already been considered. Perturbing a nonsignal source node may theoretically cause a difference, however the nodes in this category in the Sun et al. model represent molecules that are abundant in the cell or cell environment, thus their perturbation is not biologically relevant or practical.
Next we consider the anion nodes reduced due to the simplified stomatal opening rule. Recall that these nodes do not affect other nodes except stomatal opening in the long term. There cannot be multi-stability in anion nodes unless the assumptions of sufficient initial [NO 3 − ] a and starch concentration, and sufficient initial mitochondrial TCA cycle activity are violated (details are provided in Additional file 3, section 5 and 6). Since there is no support for interventions that would lead to the violation of these assumptions, it is reasonable to conclude that no multi-stability can be found in the reduced nodes under biologically relevant situations. We also found that there can be an additional oscillation in the RIC7 path (involving the nodes ROP2, RIC7 and SO) when a special set of perturbations is applied. Under that case, the nodes RIC7 and SO will oscillate. Since the effect of this behavior is small (within 5 % of the unperturbed SO value in the Sun et al. model [15]), it has little biological significance. There are no more possible oscillations as there are no more negative feedback loops. To conclude, the original Sun et al. model has oscillations only in cytosolic Ca 2+ ([Ca 2+ ] c ) and Ca 2+ ATPase, and has multi-stability only in PMV and K out , under situations that are biologically meaningful.

Discussion
The conclusions we obtained can tell us how to control this network model. Generally in engineering applications, control means to drive a system into an arbitrary state [36,37]. However in biological systems, it is more meaningful to drive the system into one of its natural attractors rather than into an arbitrary state, as the attractors correspond to stable phenotypes [38]. To control the attractor of a Boolean system, one needs to control only its input nodes and a subset of nodes in each stable motif [39]. Our integrated analysis, involving Boolean conversion, indicates that to control the attractor that the stomatal opening network evolves into, one only needs to control the input signals and PMV, even in case of perturbations. In particular, to control the stomatal opening value, one only needs to control the input signals, under any perturbation.
The reduced model provides new biological insights. Normally, when ABA is present, stomata will close. However in some knockout mutants stomata can open to a certain extent in the presence of ABA, although the opening level is not as much as in the case without ABA [15]. Such partial reversals of the effect of ABA are important for understanding the mechanism of stomatal opening. For example, Sun et al. reported that OST1 knockout (OST1 is kept 0) and inhibition of the NADPH oxidase (AtrbohD/F is kept 0) yielded partially restored SO level in simulations, in agreement with experimental observations (see Additional file 2 for the comparison of the equivalent simulations in the reduced model with experiments). Simplification of the Sun et al. model allows easier simulation of more perturbation scenarios, e.g. the systematic identification of possible partial reversals. Table 6 indicates all the partial reversals due to single node knockouts in the reduced model.
Our results reproduce the observation that knockout of nodes in the ABA pathway (PLD, NO, ROS) can cause partial reversals of ABA's effect. We find that AnionCh knockout can partially restore stomatal opening inhibited by ABA, a result not reported by Sun et al., but which is supported by experimental evidence [40]. In addition, Table 6 offers a new biological prediction: low CO 2 concentration can partially restore stomatal opening when ABA is present. This is consistent with the knowledge that CO 2 -free air promotes stomatal opening in the absence of ABA [41]. This CO 2 effect suggests a mechanism of cross-talk between CO 2 and ABA. Importantly, apart from the five nodes listed in Table 6, no other node's knockout can reverse ABA's inhibition of stomatal opening. The perturbation results of Table 5 offer many more new predictions.
Our combination of techniques offers a powerful framework for determining the dynamic repertoire of a multi-level dynamic model. Multi-level models are more accurate than Boolean models in describing the quantitative characteristics of dynamic systems, but there are few general methods to analyze multi-level models [10,12]. By combining different existing methods, we were able to overcome the limitations of each method. Our successful combination of existing methods offers a promising way to analyze multi-level models, and might point towards a general strategy to analyze the attractors of multi-level models, biological or non-biological.
A notable future direction for this work is to develop an alternative way to determine the attractors of multilevel models by extending the concept of stable motifs. Compared with conversion to a Boolean model, then applying Boolean stable motif algorithm, extending the stable motif algorithm to multi-level models can avoid potential attractor change issues. Development of such a technique will allow easy and powerful attractor analysis for multi-level models.

Conclusions
We obtained a very strong conclusion about the attractors of the Sun et al. stomatal opening model: under any combination of sustained signals, all nodes in the model converge into steady states, with the potential exception of the cytosolic Ca 2+ ([Ca 2+ ] c ) and Ca 2+ ATPase. Variations in the initial condition of non-source nodes or in process timing (node update sequence) can drive at most two nodes, PMV and K out , into a different attractor. This high degree of attractor similarity is somewhat unexpected, as the network has a large strongly connected component and several feedback loops. Thus, despite the decidedly non-linear structure of the network, most parts of the system behave in the consistent manner of a linear pathway. This is a distinct feature of the stomatal opening model: many dynamic models of biological systems have multiple, diverse attractors [22,42]. The models of these systems will evolve into drastically different attractors when starting from different initial conditions, sometimes even when starting from the same initial condition, demonstrating different biological trajectories. In the stomatal opening model, however, the uniqueness of the steady state stomatal opening level suggests that the final extent of the stomatal opening response is robust and resilient Red Light 0 3 1 The first set of columns, with the header 'Light, CO 2 and ABA condition' , indicate the input signal combinations. The 2nd column is the stomatal opening without perturbations. The 3rd column set indicates the nodes whose knockout would yield a stomatal opening level that is higher than the unperturbed value of 0. CO 2 knockout means CO 2 being set to zero (CO 2 free air). No entry means the setting does not cause partial reversal