Our model differs from previous models employed in the life sciences in the following fundamental aspects. First, we have reconstructed the signaling network from inferred indirect relationships and pathways as opposed to direct interactions; in graph theoretical terminology, we found the minimal network consistent with a set of reachability relationships. This network predicts the existence of numerous additional signal mediators (intermediary nodes), all of which could be targets of regulation. Second, the network obtained is significantly more complex than those usually modeled in a dynamic fashion. We bridge the incompleteness of regulatory knowledge and the absence of quantitative dose-response relationships for the vast majority of the interactions in the network by employing qualitative and stochastic dynamic modeling previously applied only in the context of gene regulatory networks [53].

Mathematical models of stomatal behavior in response to environmental change have been studied for decades [63,64]. However, no mathematical model has been formulated that integrates the multitude of recent experimental findings concerning the molecular signaling network of guard cells. Boolean modeling has been used to describe aspects of plant development such as specification of floral organs [65], and there are a handful of reports describing Boolean models of light and pathogen-, and light by carbon-regulated gene expression [66–68]. Use of a qualitative modeling framework for signaling networks is justified by the observation that signaling networks maintain their function even when faced with fluctuations in components and reaction rates [69]. Our model uses experimental evidence concerning the effects of gene knockouts and pharmacological interventions for inferring the downstream targets of the corresponding gene products and the sign of the regulatory effect on these targets. However, use of this information does not guarantee that the dynamic model will reproduce the dynamic outcome of the knockout or intervention. Indeed, all model ingredients (node states, transfer functions) refer to the node (component) level, and there is no explicit control over pathway-level effects. Moreover, the combinatorial transfer functions we employed are, to varying extents, conjectures, informed by the best available experimental information (see Text S1). Finally, in the absence of detailed knowledge of the timing of each process and of the baseline (resting) activity of each component, we deliberately sample timescales and initial conditions randomly. Thus, an agreement between experimental and theoretical results of node disruptions is not inherent, and would provide a validation of the model.

The accuracy of our model is indeed supported by its congruency with experimental observation at multiple levels. At the pathway level, our model captures, for example, the inhibition of ABA-induced ROS production in both ost1 mutants and atrboh mutants [12,19,21] and the block of ABA-induced stomatal closure in a dominant-positive atRAC1 mutant [22]. In our model, as in experiments, ABA-induced NO production is abolished in either nos single or nia12 double mutants [13,18]. Moreover, the model reproduces the outcome that ABA can induce cytosolic K⁺ decrease by K⁺ efflux through the alternative potassium channel KAP, even when ABA-induced NO production leads to the inhibition of the outwardly-rectifying (KOUT) channel [70]. At the level of whole stomatal physiology, our model captures the findings that anion efflux [35,71] and actin cytoskeleton reorganization [22] are essential to ABA-induced stomatal closure. The importance of other components such as PA, PLD, S1P, GPA1, KOUT, pH_c in stomatal closure control [8,20,31,43,58,72], and the ABA hypersensitivity conferred by elimination of signaling through ABI1 [28,29], are also reproduced. Our model is also consistent with the observation that transgenic plants with low PLC expression still display ABA sensitivity [73].

The fact that our model accords well with experimental results suggests that the inferences and assumptions made are correct overall, and enables us to use the model to make predictions about situations that have yet to be put to experimental test. For example, the model predicts that disruption of all Ca²⁺ ATPases will cause increased ABA sensitivity, a phenomenon difficult to address experimentally due to the large family of calcium ATPases expressed in Arabidopsis guard cells (unpublished data). Most of the multiple perturbation results presented in Figure 5 and Table 2 also represent predictions, as very few of them have been tested experimentally. Results from our model can now be used by experimentalists to prioritize which of the multitude of possible double and triple knockout combinations should be studied first in wet bench experiments.

Most importantly, our model makes novel predictions concerning the relative importance of certain regulatory elements. We predict three essential components whose elimination completely blocks ABA-induced stomatal closure: membrane depolarization, anion efflux, and actin cytoskeleton reorganization. Seven components are predicted to dramatically affect the extent and stability of ABA-induced stomatal closure: pH_c control, PLD, PA, SphK, S1P, G protein signaling (GPA1), and K⁺ efflux. Five additional components, namely increase of cytosolic Ca²⁺, Atrboh, ROS, the Ca²⁺ ATPase(s), and ABI1, are predicted to affect the speed of ABA-induced stomatal closure. Note that a change in stomatal response rate may have significant repercussions, as some stimuli to which guard cells respond fluctuate on the order of seconds [74,75]. Thus our model predicts two qualitatively different realizations of a partial response to ABA: fluctuations in individual responses (leading to a reduced steady-state sensitivity at the population level), and delayed response. These predictions provide targets on which further experimental analysis should focus.

Six of the 13 key positive regulators, namely increase of cytosolic Ca²⁺, depolarization, elevation of pH_c, ROS, anion efflux, and K⁺ efflux through outwardly rectifying K⁺ channels, can be considered as network hubs [45], as they are in the set of ten highest degree (most interactive) nodes. Other nodes whose disruption leads to reduced ABA sensitivity, namely SphK, S1P, GPA1, PLD, and PA, are part of the ABA → PA path. While they are not highly connected themselves, their disruption leads to upregulation of the inhibitor ABI1, thus decreasing the efficiency of ABA-induced stomatal closure. Similarly, the node representing actin reorganization has a low degree. Thus the intuitive prediction, suggested by studies in yeast gene knockouts [76,77], that there would be a consistent positive correlation between a node's degree and its dynamic importance, is not supported here, providing another example of how dynamic modeling can reveal insights difficult to achieve by less formal methods. This lack of correlation has also been found in the context of other complex networks [78].

Comparing Figure 3 and Figure 6C, one can notice a similar heterogeneity in the measured stomatal aperture size distributions and the theoretical distribution of the cumulative probability of closure in the case of multiple node disruptions. While apparently unconnected, there is a link between the two types of heterogeneity. Due to stochastic effects on gene and protein expression, it is possible that in a real environment not all components of the ABA signal transduction network are fully functional. Therefore, even genetically identical populations of guard cells may be heterogeneous at the regulatory and functional level, and may respond to ABA in slightly different ways. In this case, the heterogeneity in double and triple disruption simulations provides an explanation for the observed heterogeneity in the experimentally normal response: the latter is actually a mixture of responses from genetically highly similar but functionally nonidentical guard cells.