Model-driven discovery of calcium-related protein-phosphatase inhibition in plant guard cell signaling

The plant hormone abscisic acid (ABA) promotes stomatal closure via multifarious cellular signaling cascades. Our previous comprehensive reconstruction of the stomatal closure network resulted in an 81-node network with 153 edges. Discrete dynamic modeling utilizing this network reproduced over 75% of experimental observations but a few experimentally supported results were not recapitulated. Here we identify predictions that improve the agreement between model and experiment. We performed dynamics-preserving network reduction, resulting in a condensed 49 node and 113 edge stomatal closure network that preserved all dynamics-determining network motifs and reproduced the predictions of the original model. We then utilized the reduced network to explore cases in which experimental activation of internal nodes in the absence of ABA elicited stomatal closure in wet bench experiments, but not in our in silico model. Our simulations revealed that addition of a single edge, which allows indirect inhibition of any one of three PP2C protein phosphatases (ABI2, PP2CA, HAB1) by cytosolic Ca2+ elevation, resolves the majority of the discrepancies. Consistent with this hypothesis, we experimentally show that Ca2+ application to cellular lysates at physiological concentrations inhibits PP2C activity. The model augmented with this new edge provides new insights into the role of cytosolic Ca2+ oscillations in stomatal closure, revealing a mutual reinforcement between repeated increases in cytosolic Ca2+ concentration and a self-sustaining feedback circuit inside the signaling network. These results illustrate how iteration between model and experiment can improve predictions of highly complex cellular dynamics.


Introduction
Stomatal pores on the surfaces of aerial plant parts mediate both uptake of carbon dioxide (CO 2 ), the substrate for photosynthesis, and transpirational water vapor loss. The aperture of each stomatal pore is increased or decreased by dynamic swelling or shrinking, respectively, of the surrounding guard cell pair. Given the central role that stomatal apertures play in control of both carbon assimilation and plant water status, guard cells have evolved exquisite sensitivity to a plethora of environmental and endogenous signals. In particular, drought and other dehydrating stresses result in abscisic acid (ABA) synthesis and redistribution that comprises a strong stimulus for both inhibition of stomatal opening and promotion of stomatal closure.
Guard cell researchers can directly observe responses to ABA in the form of electrophysiological measurements of the ionic fluxes that drive guard cell volume changes, microscopybased direct observations of stomatal apertures, and whole-leaf as well as whole-canopy measurements of alterations in CO 2 and water vapor exchange with the atmosphere. Owing to these advantages, as well as to the vital role of this specialized cell type in the control of plant water use efficiency, guard cells have become a leading system for elucidation of cellular signaling processes [1][2][3]. A core ABA response system in guard cells consists of ABA binding to the soluble pyrabactin resistance 1/pyrabactin resistance 1-like/regulatory component of ABA receptor (RCAR/PYR1/PYL) receptors that inhibit protein phosphatase 2Cs (PP2Cs), thus activating the Ca 2+ -independent protein kinase OPEN STOMATA 1 (OST1). Resultant activation of nicotinamide adenine dinucleotide phosphate (NADPH) oxidases results in reactive oxygen species (ROS) production that enhances Ca 2+ uptake. Dozens of downstream processes, both Ca 2+ -dependent and Ca 2+ -independent, ultimately activate anion channels; consequent anion loss causes membrane depolarization that drives K + efflux, and the net loss of ions drives water efflux, guard cell shrinkage and stomatal closure.
In total, over 80 distinct signaling entities, including enzymes, signaling proteins, metabolites, cytosolic pH (pH c ), cytosolic Ca 2+ (Ca 2+ c ), and cellular redox status are implicated by experimental analyses as components in ABA-induced stomatal closure. This is a massive amount of information, and computational approaches can aid in its comprehension and synthesis. We recently published an updated network of ABA-induced stomatal closure based on synthesis of information from over 100 scientific publications. The network contains 81 components (nodes) and 153 connections (edges) [4]. While networks themselves are static entities, they can form the basis of dynamic modeling approaches that can provide both analysis and predictions. Appropriate choice of modeling approaches requires recognition of the types of experimental data and their limitations. For example, in the case of ionic fluxes, detailed quantitative information in real time is available, making quantitative kinetic modeling approaches such as those invoked by the OnGuard software appropriate [5][6][7]. However, in the case of guard cell signaling networks, much of the information is qualitative rather than quantitative. In such instances, discrete models, in which components are designated to exist in a discrete, rather than continuous, number of states, are more appropriate, of which the simplest is the Boolean ON/OFF formalism. In this formalism the ON state (alternatively denoted as 1) is interpreted as a higher-than-threshold abundance or activity, and the OFF or 0 state is interpreted as lower-than-threshold abundance and activity.
In our previous model [4], we synthesized published knowledge concerning the mechanisms by which guard cell components regulate, and are regulated by, each other during ABA signaling, resulting in 58 Boolean regulatory functions that describe this regulation with the Boolean operators AND, OR, and NOT. We established initial (resting) node states where known, randomized states of the remaining nodes, and performed dynamic modeling in both the presence and absence of ABA. Because little is known concerning the relative temporal relationships of activation/deactivation among the majority of the 80+ signaling components, which function in highly interconnected signaling pathways, our modeling employed randomized asynchronous update, in which nodes are updated sequentially in a randomly selected order. We performed a large number of simulations (considering a large number of in silico stomata) and determined the percentage of simulations wherein the final state (attractor) representing stomatal closure was reached. Taking advantage of the ease by which components can be manipulated in silico, we also predicted the impact of knock out or constitutive activation (as can occur via pharmacological intervention or overexpression of an active protein) of every node in the network. Our previous model had 2 81 possible system states (i.e., more than 10 23 ) and thus was far from determined by the available wet bench data. Nevertheless, we observed that modeling predictions were congruent with experimental evidence in 85 out of 112 cases (76%), supporting the validity of our model and the usefulness of discrete dynamic modeling for complex non-linear systems, including ABA-induced stomatal closure. From the above model, 9 of the predictions that were inconsistent with experimental data were simulations of stomatal behavior upon constitutive activation of an internal node in the absence of ABA. For example, while the model correctly predicted that elevation of ROS would lead to stomatal closure even in the absence of ABA, which is well-documented experimentally [8][9][10], the model failed to correctly predict that elevation of Ca 2+ c would lead to stomatal closure in the absence of ABA, even though manipulations that elevate Ca 2+ c are confirmed strong drivers of stomatal closure [11][12][13][14][15]. We determined that for 8 of these 9 cases, the discrepancy with experimental results would be eliminated if we made just one additional assumption, namely that the four PP2Cs with a documented role in guard cell ABA signaling are inhibited by elevation of Ca 2+ c levels. These four PP2Cs are ABSCISIC ACID INSENSITIVE 1 (ABI1), ABSCISIC ACID INSENSITIVE 2 (ABI2), HYPERSENSITIVE TO ABA 1 (HAB1) and PROTEIN PHOSPHATASE 2CA (PP2CA). In the present report, we reduce our original model into a smaller network that preserves all the relationships of the complete network while being intuitively more accessible. We then evaluate, both computationally and experimentally, the previously posed hypothesis that Ca 2+ c either directly or indirectly inhibits guard cell PP2Cs. Our analysis also highlights the network motif whose activation underlies the efficacy of the observed closure stimuli. Deeper analysis of the relationship of this motif with Ca 2+ c oscillations provides insight into how various regular or irregular Ca 2+ c patterns can affect stomatal closure.

Model reduction
Our published network model of ABA-induced stomatal closure contains 81 nodes, 153 edges, and 58 regulatory functions (23 nodes do not have regulators) [4]. The full biological names corresponding to the abbreviated node names of this network are reproduced in S1 Table. To ease comprehension of our new analyses, we employ dynamics-preserving methods to reduce this network [16][17][18]. These manipulations are possible given the inherent logic structure of Boolean networks, can be implemented for any Boolean system, and have no impact on the outcomes of dynamic modeling. All network reduction processes that were implemented are described in the Methods section and in S2-S5 Tables. To give an example of one type of reduction, two nodes related by an equivalence relationship can be merged. That is, given an edge A ! B such that A is the sole regulator of B and B is the only target of A, the state of B will equal the state of A in every long-term state (steady state) of the system. Because of this strong determination between node A and node B, they can be merged into a single equivalent node AB. The steady state of this merged node reflects the shared steady state of node A and B, and it also equivalently expresses interventions into nodes A or B. For example, since knockout of node A causes the OFF state of node B, the knockout of node AB in the reduced network is equivalent to the knockout of A, and also to the knockout of B, in the original system. The instances of node merging employed in the network reduction are supported by 40 logic equivalence relationships (see Table 1), and by the vast majority of the experimental observations of the effect of knockout or constitutive activation of a node on ABA induced closure (see S2-S5 Tables). This node merging was not done in our original model so as to fully explicate and acknowledge all the biological information in the literature; however, such network reduction simplifies and increases the comprehensibility of the model. Ultimately, network reduction resulted in the equivalent network depicted in Fig 1, with 49 nodes (of which 14 are merged nodes) and 113 edges. The regulatory functions of the reduced model are indicated in S1 Text. For convenience, in Fig 1 and in the following descriptions we use a shorthand notation for the merged nodes; Table 1 indicates the full notation of each merged node.
Similar to the complete ABA signaling network, the reduced network has a large feedbackrich center, which in graph theory is called strongly connected component (SCC, see Methods). As shown in Fig 2, the SCC contains 27 nodes and 50 edges, in other words almost half of the network. Four nodes, including ABA, can reach the SCC through paths (i.e. they are in the in-component) and 16 nodes, including Closure, can be reached from the SCC (i.e. they are in the out-component); these nodes are listed in S6 Table. The nodes of the SCC that are directly regulated by nodes in the in-component are shown in fuchsia background in Fig 2A while the nodes that directly regulate nodes in the out-component are shown in yellow background. Four nodes, namely, ABI1, ABI2, PP2CA, and Vacuolar Acidification, are regulated by a node in the in-component and also directly regulate nodes in the out-component. Three of these nodes are PP2Cs, for which both the edges from the in-component and the edges to the outcomponent are inhibitory. These nodes need to be inactivated either by nodes of the in-component (in response to ABA) or by nodes of the SCC such as ROS, phosphatidic acid (PA), or as proposed here, Ca 2+ c , in order to achieve closure. One of the important organizational features of the ABA-induced stomatal closure network is its stable motifs. Stable motifs are generalized positive feedback loops that maintain a fixed state of their constituent nodes once stabilized [19]. The effects of this stabilization propagate through the network, and channel the dynamics of the system into a subset of its full state space. Ultimately, a stable motif or a combination of stable motifs determine the attractor(s) (long-term dynamic behaviors) of the system. We verified that the reduced model preserved the original stable motifs (compare S1 Fig with Fig 2B-2E of [4]). Multiple stable motifs are associated with ABA-induced stomatal closure (S1A Fig), and a family of related motifs describe the lack of closure in the absence of ABA (S1B and S1D Fig). The reduced model also retained a stable motif associated with closure in the absence of ABA (S1C Fig and orange nodes in Fig 2B), which, consistent with biological reality, does not stabilize in any of the trajectories of the simulated wild type system. A key reason for the lack of realization of the stable motif is that in the absence of ABA the PP2Cs are active in all trajectories.
The reduced model's attractors (determined as described in the Methods section and summarized in S7 Table) are equivalent to the attractors of the original model, as expected. In the attractor associated with closure in the presence of ABA, 30 nodes stabilize in the ON state, 9 nodes stabilize in the OFF state, and 9 nodes oscillate. The oscillating nodes are driven by the negative feedback loop formed by Ca 2+ c increase and Ca 2+ ATPase (the node with a function to decrease Ca 2+ c ). These nodes are shown in purple background in Fig 2B. Indeed, the existence of Ca 2+ c oscillations in ABA-induced stomatal closure is well-established [11,20]. In the simulations, these 9 nodes switch between ON and OFF periods of one to three timesteps, with a variability in duration caused by the stochasticity of the model (see Methods). The average length of the ON and OFF periods of these nodes (see Methods) show clear regularities, as described in S2 Text. Ca 2+ c and Ca 2+ ATPase have average ON and OFF periods of 1.33 time steps. This is the same value as the average ON/OFF period of the nodes of a two-node S1P GCR1 = 0, GPA1 = 1, SPHK1/2 = 1 SPP1 = 1, Sph = 0, GPA1 = 0, SPHK1/2 = 0 The first column lists the coded node names of the nodes that are reduced or merged in the reduction process, following the naming procedures described in the Methods. Explication of the meaning of the coded node names is given in S2-S5 Tables. The second column lists the shorthand notation for the nodes, used for convenience in the figures and tables. The third column lists the settings (in terms of a sustained state of a node that was eliminated during the network reduction process) that yield the sustained ON (1) state of the merged node. Therefore, simulated constitutive activation (i.e., sustained ON state) of the merged node incorporates all of these individual settings. The fourth column lists the settings that yield the sustained OFF state of the merged node; simulated knockout of the merged node incorporates all these settings. The settings that were experimentally studied are indicated in boldface in the last two columns (see S2-S5 Tables for more details [4]. Edges ending in an arrowhead denote an activating process or regulation while those ending in filled circles denote an inhibitory regulation. Nodes that are enzymes are shown in red color, signaling proteins in green, membrane-transport related nodes in blue and secondary messengers and small molecules in orange. Symbols with a bold boundary (and in some cases, two colors) mark nodes that merge multiple nodes of the original model; these are listed in detail in Table 1. https://doi.org/10.1371/journal.pcbi.1007429.g001 Calcium-related PP2C inhibition in stomatal signaling trisphosphate and hexakiphosphate (InsP3/6), as well as diacylglycerol (DAG), both of which are at distance 2 from Ca 2+ c , have a longer average ON/OFF period than the nodes directly regulated by Ca 2+ c . The increase in the average ON/OFF period with increasing distance from Ca 2+ c reflects the increased chance that the effect of a transient change in Ca 2+ c does not translate into activation of these downstream nodes (see S2 Text). In summary, the attractor corresponding to closure in the presence of ABA features a constant state for 39 nodes and three classes of oscillations for 9 nodes. This attractor signature serves as a comparison for all model versions considered subsequently.
There are multiple attractors associated with lack of closure in the absence of ABA, but in all of them, the same 8 nodes stabilize in the ON state and the same 31 nodes stabilize in the OFF state (see S7 Table). Six nodes stabilize in either the OFF or ON state. These nodes are directly or indirectly affected by positive self-regulation thus if they are activated during the dynamics of the system, they will maintain that activity. Finally, in the absence of ABA, three nodes affected by a negative feedback loop between depolarization and K + efflux through outwardly-rectifying channels (KOUT) either oscillate or are stabilized in the OFF state.
We performed in silico simulations with the original 81-node and reduced 49-node model (see Methods and S1 Text) and determined the percentage of closure in the presence and absence of ABA and in settings where any node (except ABA and Closure) is knocked out or constitutively activated. Compilation of the results for simulated knockout and constitutive  Table) and in the attractor corresponding to closure in the absence of ABA (shown in S13 Table). The two nodes with grey background have different states in these two attractors. In both panels, edges ending in an arrowhead denote an activating process or regulation while those ending in filled circles denote an inhibitory regulation.  Table) supports the validity of network reduction. S8 Table marks the groups of nodes that are merged into a single node in the reduced model. For the majority of groups all interventions in the same group yield the same ABA response. The exceptions consist of tautological interventions such as external supply of an already-abundant molecule (the elimination of which is a minimal loss of information) and a small number of interventions into source nodes that have multiple targets or whose targets have multiple regulators. The reduced model reproduces the results of the original model [4] with 10 minor discrepancies (all involving the speed of reaching 100% closure) out of the 188 possible knockout or constitutive activity settings, see S9, S10 and S11 Tables. A characteristic example of discrepancy between the full and reduced models is the knockout of the nitrate reductases NIA1/2 in the presence of ABA. Knockout of nitrate reductase activity is predicted to lead to a close to wild type response in the reduced model, consistent with experimental evidence [21], while the full model predicted ABA hypersensitivity. Collectively, these results verify that the reduction methods preserve the dynamics of the system.

Model predictions concerning Ca 2+ c inhibiting one or more PP2Cs
Our original model successfully recapitulated 13 instances of experimentally observed stomatal behavior upon constitutive activation of an internal node in the absence of ABA and did not recapitulate 9 experimental findings of closure [4]. We made two key observations in our previous analysis: (i) the 9 nodes whose constitutive activation can lead to closure according to experiments but not in the model participate in paths that converge to Ca 2+ release from stores and (ii) even transient inactivation of the PP2C phosphatases is able to resolve 8 discrepancies.
We concluded that analysis with the observation that predicted Ca 2+ c inhibition of all four PP2C protein phosphatases would resolve all but one case of discrepancy between experiment and model prediction (see S3 Text for a more detailed description of the argument). As illustrated in the second column of Table 2, because of node merging in the reduction process, the total number of nodes for which there are experimental data on closure in the absence of ABA in the reduced model is 17. The number of cases of agreement with experiments in the absence of ABA in the reduced model is 10 (2 of which are merged nodes) and the number of instances of disagreement with experiment is 7 (5 of which are merged nodes). As shown in the second column of Table 2, in the majority of cases of constitutive activation of a node the percentage of closure was close to the baseline (i.e., close to the case of no constitutive activation); 0% in this case. Only the constitutive activation of ROS yielded a significantly increased percentage of closure.
In this report, we newly study the relationship between Ca 2+ c and the four PP2C protein phosphatases in detail. We separately analyze various scenarios of Ca 2+ c inhibiting a single PP2C (which would be the most parsimonious assumption, see S3 Text) or multiple ones and compare model simulations with the available experimental results (Tables 2 and S12). We first performed four simulations, adding an inhibitory edge from Ca 2+ c to each one of the PP2Cs (ABI1, ABI2, HAB1 and PP2CA), with no other changes. Simulations of direct Ca 2+ c inhibition of ABI1 yielded results identical to those obtained in the absence of any additional inhibitory edge from Ca 2+ c (third column of Table 2). We uncovered that this is because an indirect inhibitory effect between Ca 2+ c and ABI1 was already incorporated by multiple pathways in the model, such that the newly posited edge, in the case of ABI1, is completely redundant. One such indirect pathway is Ca 2+ c ! PLC ! DAG ! PA-• ABI1; this pathway is discussed in more detail in the next section. In contrast, when Ca 2+ c is assumed to inhibit any one of ABI2, HAB1 or PP2CA we observed that the constitutive activation (CA) of Ca 2+ c can lead to complete stomatal closure (fourth column of Table 2). This computational outcome agrees with multiple observations that experimental elevation of Ca 2+ c suffices to elicit stomatal closure [11,13,15]. Moreover, addition of an inhibitory edge from Ca 2+ c to any one of ABI2, HAB1, or PP2CA also successfully captures all 9 cases of experimentally observed closure due to constitutive activation of a node, including pH c , nitric oxide (NO), cyclic ADP-ribose (cADPR), PA, and InsP3 [21,25,27,29,36]. Fig 3 compares the percentage of closure for constitutive activation of the node that merges sphingosine, sphingosine kinases and sphingosine-1-phosphate (S1P), the node that merges 8-nitro-cyclic GMP, ADP ribosyl cyclase and cADPR (denoted cADPR), ROS (which is a merged node representing reactive oxygen species and the enzyme RBOH), PA, Ca 2+ c in the reduced model (panel A) and the model version that assumes that Ca 2+ c inhibits ABI2 (panel B). Constitutive activation of ROS yields 100% closure in both model versions. The model version that assumes that Ca 2+ c inhibits ABI2 outperforms the original reduced model by yielding an increased closure percentage when S1P or PA are constitutively activated, and 100% closure for constitutive activation of cADPR, InsP3/6, Ca 2+ c . All these results are consistent with experimental knowledge [15,21,25,27,[29][30][31]36]. An explanation of the two categories of increase in the percentage of closure is provided in a later section.

Calcium-related PP2C inhibition in stomatal signaling
This attractor is reachable due to the possibility of activating the stable motif corresponding to closure in the absence of ABA ( Fig 2B); thus, trajectories that reach this attractor yield a nonzero baseline percentage of closure. In contrast, this attractor is not reachable in the reduced model, nor in the reduced model version wherein we add a direct inhibitory edge from Ca 2+ c to ABI1, thus the baseline percentage of closure is 0 in these models. The mechanisms of activating the stable motif underlying closure in the absence of ABA are explained in detail in the next section.
The response categories of the constitutive activation simulation results are the same whether any two or more than two PP2Cs are inhibited by Ca 2+ c . The baseline probabilities of closure are slightly different in each model version, see S12 Table. A main contributor to the similarity of the model outcomes is the symmetry of the regulatory function of OST1, which incorporates the assumption that the inactivity of any two PP2Cs is sufficient to activate OST1 (see S2 Text of [4]). As detailed in the next section and in S4 Text, OST1 is a part of a stable motif that drives the system to closure, thus the equivalence of two settings in activating OST1 translates to their being equally effective in activating the stable motif and yielding closure.
While the motivation underlying our hypothesized inhibition of the PP2Cs is to accurately predict experimental results of closure in the absence of ABA upon constitutive activation of key internal nodes, it is important to check that the thus-augmented model maintains agreement with experimental results in the presence of ABA. As summarized in S14 Table, the assumed inhibition introduces a single discrepancy: the knockout of the node RCARs yields a milder defect in the model than in experimental observations (see the caption of S14 Table for a discussion of this discrepancy).

Theoretical analysis of the stable motif associated with closure in the absence of ABA
For a deeper understanding of how the inhibition of PP2Cs by Ca 2+ c can yield ABA-independent closure, we used causal logic analysis [16] to determine the drivers of the stable motif associated with closure in the absence of ABA (see Methods). Stabilization of a stable motif (or successive stabilization of a group of stable motifs) can drive an entire network to a particular attractor. These motifs can be identified by finding strongly connected components (SCCs) in an expanded network that expresses the regulatory function of each node [19]. Two different stable motifs were identified in our original model [4] in the absence of ABA (and reproduced by the reduced model, see S1 Fig), one of which is associated with closure in the absence of ABA. This stable motif, which can be stabilized only after the stabilization of the Vacuolar Acidification motif, is shown in reduced form in Fig 4A. The stable motif primarily contains the nodes ROS, PA, PLDδ (Phospholipase D δ), pH c , OST1, S1P, and the PP2Cs ABI1 and ABI2. We found that the node ROS is a driver node of this motif, which means that if ROS is set to the ON state, the entire motif stabilizes and eventually leads to the closure attractor of the network. There also exist three-node combinations whose stabilization can lead to the stabilization of the motif. For example, PLDδ, pH c and ABI2 make up a three-node driver set, that is, if PLDδ and pH c are fixed to ON and ABI2 is fixed to OFF, then the entire motif stabilizes, leading to closure. In the reduced model, ROS is the only node in the network whose constitutive activation can yield the stabilization of the motif and closure in the absence of ABA, consistent with the results in Fig 3A and S11 Table. We next studied the model version wherein Ca 2+ c directly inhibits ABI2. In this version, Ca 2+ c becomes an external driver node of the motif associated with closure in the absence of ABA (Fig 4B). Moreover, as described in S4 Text, we find that a Ca 2+ c pulse (sustained ON state) of four or more time steps, or multiple repeated Ca 2+ c spikes (where one spike is equated We also found that when InsP3/6 or cADPR (two of the nodes that are upstream of Ca 2+ c and whose constitutive activation, computationally or experimentally, leads to closure) are in their ON state, this eventually leads to Ca 2+ c exhibiting a sustained oscillation pattern. Our model predicts that the constitutive activation of InsP3/6 or cADPR in the absence of ABA leads to sustained Ca 2+ influx from internal stores to the cytosol (CIS), in agreement with experiments [25,36]. The model also predicts that Ca 2+ c then follows the oscillatory pattern described in S2 Text, which yields the stomatal closure attractor in 100% of our simulations. In total there are eight nodes whose constitutive activation can elicit cytosolic Ca 2+ oscillations, which then lead to the stabilization of the stable motif associated with closure, and ultimately yield stomatal closure: Actin reorganization, CaIM, CIS, cADPR, GUARD CELL HYDRO-GEN PEROXIDE RESISTANT 1 (GHR1), InsP3/6, PLC, ROS. Indeed, these eight nodes and Ca 2+ c populate the category "significantly increased closure" in the absence of ABA (see S12 Table). These results support the view that Ca 2+ c is a key driver of the process of stomatal closure in the absence of ABA [11,13,15].

Experimental assessment of PP2C inhibition by Ca 2+
Our simulation results and their agreement with reported experimental results gave us a strong basis for predicting that Ca 2+ c inhibits the PP2Cs. This prediction was congruent with the observation that ABI1 has a candidate EF-hand Ca 2+ -binding motif [37,38]. Indeed, the presence of this motif previously inspired three independent studies of potential direct Ca 2+ regulation of ABI1, all of which assayed the phosphatase activity of recombinant ABI1 in vitro. Giradaut and co-workers observed essentially the same extent of Ca 2+ inhibition of phosphatase activity for recombinant full-length ABI1 and truncated ABI1 interrupted in the EF hand motif [39]; in both instances, inhibition occurred in the mM Ca 2+ range, leading to the conclusion that Ca 2+ was inhibiting the catalytic domain, rather than binding the EF-hand. Grill and colleagues similarly observed Ca 2+ inhibition of recombinant ABI1 only at unphysiologically high Ca 2+ concentrations, with a half-inhibitory concentration of 8 mM; they suggested that Ca 2+ may be hindering the Mg 2+ binding required for phosphatase activity [40]. Finally, Sheen stabilization of the Vacuolar Acidification motif (and the corresponding ON state of Vacuolar Acidification) is a required precursor of this stable motif in all versions of the model. (A) The stable motif in the absence of any additional inhibitory edge to any of the PP2Cs. ROS (shown in orange color) is a single-node driver of this motif, i.e., if ROS has a sustained ON state, it can stabilize the entire motif, hence eventually leading the network to the Closure = ON attractor. PA, pH c and ABI2 (purple nodes) form a three-node driver set, that is, if PA is ON, pH c is ON and ABI2 is OFF, then the entire motif stabilizes, leading to closure. Note that not every combination of three nodes belonging to the motif shown in (A) can stabilize the entire motif; see S5 Text for more information. (B) The corresponding stable motif in the model version where Ca 2+ c is assumed to inhibit ABI2. In this case, if Ca 2+ c is set to ON, it stabilizes PA = ON, pH c = ON, and ABI2 = OFF (through its edges incident on nodes of the motif), which forms a three-node driver, stabilizing the entire motif, hence leading to closure. (C) The stable motif in the model version where PA is assumed to inhibit ABI2. PA and Vacuolar Acidification (shown in green color) make up a two-node driver set of the motif, as PA inhibits ABI1 and ABI2, which together with Vacuolar Acidification activate pH c , making this combination equivalent to the PA = ON, pH c = ON, ABI2 = OFF three-node driver of the motif. The Ca 2+ c ! PLC ! DAG ! PA path is sufficient for the activation of PA, thus Ca 2+ c remains an external driver of the stable motif associated with closure. https://doi.org/10.1371/journal.pcbi.1007429.g004 Calcium-related PP2C inhibition in stomatal signaling mutated the aspartate that is required for Ca 2+ binding in the candidate EF-hand and observed no impact on phosphatase activity [41]. In summary, these studies concluded that no direct Ca 2+ regulation of ABI1 could be observed within the physiological range prevailing in the cytosol (approximately 100 nM to low μM concentrations). The fact that these experiments were conducted in vitro under carefully controlled conditions provides no support for our hypothesis of direct Ca 2+ c inhibition of PP2Cs. We therefore considered whether the regulation of PP2Cs by Ca 2+ c could be via an indirect mechanism. Indeed, by the causal logic theory [16], almost the same system dynamics, and the same attractors, could be obtained under the assumption that the PP2C inhibition by Ca 2+ c occurs indirectly, via a pathway that has the same logical implication as a direct edge from Ca 2+ c to the PP2Cs. Accordingly, we decided to assay whether physiological levels of Ca 2+ in the presence of cellular (protoplast) lysates without the influence of high apoplastic Ca 2+ would inhibit PP2C activity. As shown in Fig 5, we did indeed observe inhibition of PP2C activity in the lysates by low μM levels of Ca 2+ , corresponding to levels that prevail in the cytosol. These results, in conjunction with the prior in vitro studies [39][40][41], lend support to the idea that inhibition of PP2C activity by Ca 2+ c does occur in vivo, but that the mechanism of inhibition is indirect.

PP2C inhibition mediated by PA gives similar results
We then turned to our model to elucidate possible pathways for indirect Ca 2+ c inhibition of PP2Cs. We already know about one such indirect pathway: Ca 2+ c ! PA-• ABI1, wherein Ca 2+ c contributes to PA production by upregulating PLDα and PLC (Fig 6A) [42][43][44][45][46], and PA is known to directly bind to and inhibit one of the PP2Cs (ABI1) [47]. Thus, we next considered the alternative hypothesis that PA binds to and inhibits additional PP2Cs (see S3 Text) and evaluated this hypothesis by simulations and stable motif analysis. We found that this assumption led to almost identical results as assuming that Ca 2+ c directly inhibits the PP2Cs. The attractors of the model in the presence or absence of ABA in the case of inhibition of the PP2Cs by PA are identical to those in case of inhibition of the PP2Cs by Ca 2+ c , shown in S13 Table. The stable motif corresponding to closure in the absence of ABA is identical except for the replacement of the Ca 2+ c−• ABI2 edge with a PA-• ABI2 edge (Fig 4C). ROS and Ca 2+ c remain internal and external drivers, respectively, of this stable motif. The only difference is that PA becomes an internal driver of the stable motif if Vacuolar acidification is already stabilized at ON. Moreover, we found the same categorization of the constitutive activation of nodes as in the case where Ca 2+ c directly inhibits the PP2Cs (fifth column of Tables 2 and S14). This means that the two model versions are equally supported by experimental results. As a further illustration of this point, Fig 6B compares Fig 6B) has the same effect in both versions of the model. This is because constitutive activation of Ca 2+ c yields constitutive activation of PA. The effect of constitutive activation of multiple other nodes is very similar as well. Constitutive activation of PA yields a stronger effect in the model version where PA inhibits the PP2Cs. This is because constitutive activation of PA does not imply constitutive activation (sustained high concentration) of Ca 2+ c , nor of CIS or CaIM. Thus, constitutive activation of PA yields a lower probability of activating the stable motif associated with closure, and thus a lower percentage of closure, in the model version where Ca 2+ c directly inhibits the PP2Cs than in the model version where PA directly inhibits the PP2Cs.
In summary, our modeling outcomes predict that the known inhibition of ABI1 by PA, plus putative inhibition of one or more of the other three PP2Cs by PA, would suffice to equal Calcium-related PP2C inhibition in stomatal signaling the outcome of positing direct inhibition of the same PP2Cs by Ca 2+ c in terms of system dynamics. Calcium-induced PLC activation, which, according to our model, is upstream of PA production, occurs in both guard cells and mesophyll cells [46], thus the Ca 2+ c ! PA-• PP2C mechanism is consistent with, and a likely contributor to our experiments in cellular lysates. For ABI1, lipid binding experiments with truncated ABI1 constructs have demonstrated that PA binding occurs within the first 104 amino acids of the N-terminus, with the region of residues 67 to 97 being of particular interest given sequence similarity with known PA-binding regions identified in mammalian systems [47]. We therefore analyzed the similarity of ABI2, HAB1, and PP2CA to this region in ABI1. Of these, ABI2 has the greatest similarity to ABI1 in the first 104 amino acids (36% identity and 47% similarity based on pair-wise alignment), including in the domain from residues 67 to 97 (33% identity and 83% similarity) making ABI2 the most likely additional candidate for inhibition by PA.

The model predicts an interplay between the stable motif and Ca 2+ c oscillations
As described earlier, regularly or irregularly repeated Ca 2+ c spikes yield the consistent activation of the closure-inducing stable motif, which then leads to closure, in all model versions where Ca 2+ c inhibits (directly or indirectly) at least one PP2C other than ABI1. In all model versions including the original model [4], ABA leads to sustained Ca 2+ influx through the membrane through the ABA-• AtRAC1 -• Actin reorganization ! CaIM path [32,48] and also contributes to Ca 2+ c induced Ca 2+ release from internal stores [49]; thus in the presence of ABA, Ca 2+ c exhibits sustained oscillations (see S2 Text), which then drive closure [20]. We next sought to determine whether it is possible that such closure-driving Ca 2+ c patterns can be generated in the absence of ABA.
Through simulations and causal logic analysis, we find that in the absence of ABA, a sufficient condition for Ca 2+ c elevation is an initial condition wherein at least one of cADPR, GHR1, AtRAC1 (small GTPase RAC1), PLC (in the model version where Ca 2+ c inhibits ABI2) or at least one of cADPR, GHR1, AtRAC1, PLC, PLDδ, DAG (in the model version where PA inhibits ABI2) is initiated in their state corresponding to stomatal closure (see S7 Text). Notably, AtRAC1 is inactivated during stomatal closure, while all the other nodes are activated (see S7 and S13 Tables). The conclusion that initial inactivity of AtRAC1 contributes to Ca 2+ c increase and to stomatal closure is supported by the experimental observation that a loss of function AtRAC1 mutant induced actin reorganization and stomatal closure [32]. Our analysis indicates that a necessary condition for repeated Ca 2+ c elevation is the sustained activation of at least one of the three positive feedback loops incident on Ca 2+ c . As shown in Fig 7A, these feedback loops are Ca 2+ For most of the edges participating in these feedback loops the activity of the source of the edge is sufficient to maintain the activity of the target (e.g. the activity of PLC is sufficient for the production of DAG). The PA! ROS edge is an exception; the production of ROS has multiple other conditions, which are simultaneously satisfied only when the stable motif associated with closure has stabilized (see Fig 7A and S7 Text). This edge, and thus the stable motif, participates in two of the three Ca 2+ c feedback loops. PLC participates in all three Ca 2+ c positive feedback loops, thus in the absence of ABA its simulated knockout leads to the stabilization of Ca 2+ c in the OFF state. This prediction is supported by the experimental observation that pharmacological inhibition of PLC abolished oscillations in Ca 2+ c induced by pathogenic elicitors [50,51]. We find that PLC and DAG are also important mediators of the effect of Ca 2+ c on the stable motif associated with closure: the closure-inducing effect of an externally imposed Ca 2+ c pulse or repeated Ca 2+ spikes are dependent on the Ca 2+ c ! PLC ! DAG ! PA pathway in both model versions. If this pathway is disrupted, the percentage of closure in the model reduces to its baseline value (obtained in the absence of any signal or imposed Ca 2+ c pattern). This model prediction is supported by the experimental observations that pharmacological inhibition of PLC impaired ABA-induced Ca 2+ c oscillations and abolished ABA-induced stomatal closure [20]. Our analysis indicates that there exists a mutual reinforcement between the Ca 2+ c patterns and the stable motif associated with closure in the absence of ABA, as summarized in Fig 7B. Sustained Ca 2+ c oscillations are sufficient to stabilize the stable motif (as described in S4 and S6 Texts). In turn, the stable motif mediates two positive feedback loops centered on Ca 2+ c , leading to sustained Ca 2+ c oscillations (as described above and in S7 Text). Stomatal closure necessitates the activation of the stable motif and at least a transient Ca 2+ c elevation, whose effect is then sustained through the self-regulatory stable motifs of the nodes MPK9/12 (Mitogen activated Protein Kinases 9 & 12) and CPK3/21 (CALCIUM DEPENDENT PROTEIN KINASES 3 & 21), Vacuolar acidification and Microtubule depolymerization. Thus, according c oscillations lead to stabilization of the stable motif associated with closure (see S4 Text). In turn, stabilization of that motif due to one of its internal driver sets mediates the red and blue feedback loops in panel A and leads to sustained Ca 2+ c oscillations. Initial activity of PLDδ or DAG can activate PA, which, when combined with Ca 2+ c elevation, can lead to the stabilization of the motif. This is represented by a composite node (grey circle) for which Ca 2+ c elevation for at least one time step and initial activity of PLDδ or DAG are necessary; this composite node is sufficient for stabilization of the stable motif. Initial activity of cADPR, GHR1, or PLC, or inactivity of AtRAC1 can induce Ca 2+ c elevation. Due to the stochasticity incorporated by various possible update orders, Ca 2+ c elevation has a low but non-zero probability (around 0.2) of activating the Ca 2+ c ! PLC ! InsP3/6 ! CIS loop, leading to Ca 2+ c oscillations. This is represented by the dashed edge from Ca 2+ elevation to Ca 2+ c oscillation. Stomatal closure necessitates the activation of the stable motif and at least a transient Ca 2+ c elevation (whose effect is sustained through the self-regulatory stable motifs MPK9/12, CPK3/21, Vacuolar acidification, and Microtubule depolymerization). This necessary requirement is depicted by a composite node (black circle) for which Ca 2+ c elevation and stable motifs are each necessary, while the composite node itself is sufficient for closure. https://doi.org/10.1371/journal.pcbi.1007429.g007 Calcium-related PP2C inhibition in stomatal signaling to the logic backbone of the network corresponding to stomatal closure in the absence of ABA, sustained Ca 2+ c oscillations are sufficient, but not necessary, to induce closure. This is consistent with the observation of closure induced by induced Ca 2+ c oscillations [11], and with both the Ca 2+ c dependent and Ca 2+ c independent aspects of ABA response [52][53][54].

Discussion
Here we evaluated our previously posited hypothesis that Ca 2+ c inhibits multiple proteins in the PP2C family [4]. We first derived an accessible reduced model that maintained the properties of our original model. We then found that the reduced model augmented with the assumption that Ca 2+ c inhibits one or more additional PP2C proteins agrees with all the observations that were previously discrepancies. Our combined computational and experimental analyses have led us to a refined hypothesis, in which Ca 2+ c elevation results in production of PA, which binds and inhibits ABI1 and a second PP2C, plausibly ABI2, thus engendering strong Ca 2+ c -induced stomatal closure in the absence of ABA. Comparing the performance of the two model versions in settings of constitutive activation of a node, we find that in each setting the model wherein PA is inhibiting ABI2 yields a larger or equal percentage of closure. Indeed, PA's inhibition of ABI2 may be more effective than Ca 2+ c inhibition of ABI2 since Ca 2+ c elevation cannot be sustained for an extended period because of the negative feedback loop with Ca 2+ ATPase. Because PA-binding motifs are not highly conserved [47], the possibility that one or both of HAB1 and PP2CA are inhibited by PA also cannot be ruled out. The existence of alternative or additional putative mediators of an inhibitory effect from Ca 2+ c to a PP2C also cannot be ruled out (see S3 Text). Denoting the putative mediator by X and considering, as an example, that the target PP2C is ABI2, there are two possible logic paths: Ca 2+ c activates X and X inhibits ABI2, or Ca 2+ c inhibits X, which is a required activator of ABI2. In the future, these new predictions can be further tested by a variety of wet bench approaches.
The hypothesized inhibition of additional PP2Cs by Ca 2+ c introduces a discrepancy with observations in case of knockout of the node RCARs, which denotes the ABA receptors PYR/ PYL/RCAR. In the presence of ABA, the combined knockout of multiple PYR/PYL/RCAR receptors is experimentally known to impose lack of ABA sensitivity, while the model simulations indicate hyposensitivity, reaching delayed closure (S13 Table). This discrepancy would decrease if one of the five ABA-regulated nodes with no known connections to RCARs, i.e. SPHK1/2, PI3P5K (PHOSPHATIDYLINOSITOL 3-PHOSPHATE 5-KINASE), PEPC (PHOSPHOENOLPYRUVATE CARBOXYLASE), Malate or AtRAC1 (which therefore are currently independent of RCARs by model logic) would be in fact RCARs-dependent. For example, if ABA's inhibition of AtRAC1 were mediated by RCARs, knockout of RCARs would yield reduced ABA sensitivity in the model (see caption of S14 Table). These results therefore point toward prioritization of wet bench experiments designed to test the ABA-dependence of these nodes. Alternatively, it is possible that hyposensitivity rather than insensitivity to PYR/ PYL/RCAR receptor knockout reflects biological reality. For example, a quadruple PYR/PYL/ RCAR knockout retained partial sensitivity to ABA in inhibition of stomatal opening [55], suggesting possible involvement of other ABA receptors. Whether these receptors are other members of the PYR/PYL/RCAR family or a different type of receptor remains an open question [56,57].
While the majority of the experimental analysis to date focused on stomatal closure in the presence of ABA, which was consequently the focus of our previous model [4] and logic backbone analysis [16], here we analyzed the conditions that allow closure in the absence of ABA. Our finding of inhibition of the PP2C proteins by Ca 2+ c fills a knowledge gap: augmenting the previous model with this finding restores agreement with experimental observations of closure in response to externally supplied Ca 2+ , PA, NO, S1P, cADPR, or InsP3, constitutive activity of PLDδ, or externally induced Ca 2+ c oscillations. The augmented models can serve as a foundation for future analysis of any signal or perturbation that yields stomatal closure.
For example, one interesting outcome of the augmented models is that they yield an increased baseline percentage of closure (e.g., 21.5% in Fig 3B), indicating that in a fraction of simulations the stable motif associated with closure is activated despite the absence of ABA and of constitutive activation of any nodes. In our models, seventeen nodes are initialized in a randomly selected state, simply because there is a lack of experimental data regarding their state in open stomata (see Methods). We identified and summarized in Fig 7B the multiple feedback mechanisms through which trajectories that start from initial activation of certain nodes (or inactivation in case of AtRAC1) may ultimately lead to closure. Future investigation of the resting state of these nodes will clarify the biological probability of these trajectories.
The augmented models and their analyses in response to various Ca 2+ c patterns provide insight into the mechanism through which sustained Ca 2+ c oscillations (whether regular or irregular) can drive the system to closure. Namely, the models suggest that repeated Ca 2+ c increase can be translated into a stable activation pattern through multiple self-or mutually reinforcing positive feedback loops (see Figs 2B and 7). Our computational results are consistent with experimental results indicating that defined Ca 2+ c oscillations are important for maintaining the closed state of stomata [11]. Our analysis implies that, if via other signaling pathways e.g. due to high CO 2 [58,59], pathogenic elicitors [50,51], or methyl jasmonate [60] the guard cell shows suitable Ca 2+ c patterns, it will likely stabilize in the closed stomate state. In depth elucidation of the temporal and quantitative aspects of Ca 2+ c dynamics will require a tight integration of experimental and computational investigation.

Network reduction
Source nodes are defined as nodes that regulate (have outgoing edges to) other nodes but themselves have no documented regulation, i.e. have no incoming edges. The original model contains 22 source nodes that are assumed to maintain a constant state (e.g. constant presence of a certain molecule). As proven in [18], elimination of a source node by substituting the assumed constant state into the regulatory function of its target(s) preserves the attractor repertoire of Boolean networks. Our previous study of the logic structure of a Boolean network [16] implies that this reduction does not affect the trajectories of the unperturbed system. We used logic analysis to identify the cases where the effects of perturbing a source node (e.g. the effect of depleting a normally abundant molecule) are equally captured by perturbing its target node(s). We only reduce a source node if the reduction preserves agreement with experimental evidence on perturbing the source node. There are two cases of source nodes that we reduce in this work; in each case we encode the reduced source node into the label of the resulting node. The first is the case where the source node has exactly one successor (i.e. exactly one outgoing edge). Here, the notation of the resultant node encodes whether the relationship from source node to the successor node in the original model is activating (positive) or inhibiting (negative). Examples: B[A] means that A is the source node and B is its target in the original model and there is a positive edge A ! B; C[A,~B] denotes that A and B are source nodes and C is their shared target with edges A ! C and B-• C. The list of the 17 reduced source nodes in this category along with the logic implication of each on its target is shown in S2 Table. Second is the case where the source node has multiple outgoing edges. These source nodes are eliminated only if in the original model the observed effect of knockout (KO) or constitutive activation (CA) of the source node is logically equivalent to the effect of the combined perturbation of its successor nodes. We use lowercase letters to mark these cases. For example, if A is a source node that has outgoing edges to B and C, and A is assumed to be ON in the original model, then A is removed from the network and the label B is changed to B{a} and the label C is changed to C{a}. These cases are detailed in S3 Table. Following [18] and [16], we use two ways of reducing internal nodes of the network while preserving the underlying logic of the system and hence conserving the dynamics. One of these methods operates by merging two nodes related by an equivalence relationship, that is, an edge A ! B such that A is the only regulator of B and B is the only successor of A. For example, f ADPRc = 8-nitro-cGMP, which means that the state of ADPRc follows the state of 8-nitro-cGMP with a time delay. These two nodes can be collapsed into one node, which we denote as 8-nitro-cGMP ! ADPRc. Perturbation (i.e. knockout or constitutive activation) of the resulting node equivalently represents the respective perturbation of each of the merged nodes. The 8 node pairs in this category are detailed in S4 Table. The other reduction method for internal nodes focuses on the case where a node has only one regulator, but the regulator has multiple outgoing edges. An example is an edge A ! B where A is the only regulator of B, but A has successors other than B.  Table. In addition to these reduction methods, we also merged the source nodes CPK6 and CPK23 into one node named CPK6/23 because they have the same successor nodes.

Network components
The strongly connected component (SCC) of a network represents a subnetwork in which every pair of nodes A and B has both a path from A to B and a path from B to A. The in-component of the SCC is the set of nodes that can reach the SCC through paths, and the out-component is the set of nodes that can be reached from the SCC through paths. For example, ABA is in the in-component of the SCC shown in Fig 2 while Closure is in its out-component.

Simulations
The simulations were executed using the Python library BooleanNet [61]. The SBML and Boo-leanNet files of all the model versions, along with the python code used to run the simulations, are available on the github repository /parulm/Stomata_PP2Cs. The initial state of the nodes in the reduced model is derived from the initial state of the nodes in [4], and reflects the existing knowledge regarding open stomata. Thirty-two of the 49 nodes are initialized in a particular state (ON or OFF); these nodes, their states, and the justification for the initial state are listed in S1 Text. The remaining 17 nodes, namely AnionEM, AtRAC1, cADPR, DAG, Depolarization, GHR1, KEV, KOUT, PEPC, PLC, PLDδ, QUAC1, SLAC1, SLAH3, TCTP, V-ATPase, V-PPase, are randomly initialized as ON or OFF. With these states as the starting point, the states of the nodes of the network are updated asynchronously, in a randomly selected order for 50 consecutive time steps. The order of node update is chosen randomly at each time step, giving equal probability to all update orders [62]. Update of a node means that the state of the node (ON or OFF) is updated using the Boolean regulatory function (update rule), evaluated using the current states of the regulator nodes. The Boolean update rules of the reduced network are given in S1 Text. Each additional inhibition is implemented by adding a term to the update rule of the inhibited node. Specifically, inhibition by Ca 2+ c is implemented by adding "and not Ca 2+ c " to the update rule of the target node, and inhibition by PA is implemented by adding "and not PA" to the update rule. We performed 4500 simulations, each starting with a randomly selected state of the above mentioned 17 nodes. Although this number of simulations covers a small fraction of the total possible states visited by the system, we verified that this amount of sampling is sufficient to reflect the decisions made by the system. We determined the fraction of simulations in which the node Closure is in the ON state at each time step. For easier understanding we represent this fraction as a percentage and refer to it as the percentage of closure.
To characterize the nine nodes whose state oscillates in the long-term dynamics of the system, we used an ensemble of 4000 simulations of 100 time steps. We determined the length of their ON and OFF periods (i.e. the number of consecutive time steps for which the node is ON or OFF, respectively) in the last 40 time steps of the simulation, when the system is certain to have reached the relevant attractor. We then calculated the average duration of the ON and OFF period, respectively. For example, in the particular Ca 2+ c time series [1, 1, 0, 1, 1, 0, 1, 0, 1, 1] the average ON period is 1.75 time steps and the average OFF period is 1 time step.

Simulation of node constitutive activation and knockout
The constitutive activation (CA) of a node is simulated by setting its state to ON and keeping it fixed to ON throughout the entire simulation-this can be thought of as equivalent to providing a biomolecule in non-limiting quantity or, for an enzyme, providing the enzyme in non-limiting quantity and in the active state, during an experiment. The knockout (KO) of a node is simulated by setting its state to OFF and maintaining the OFF state throughout the experiment. A separate set of 4500 simulations was run for each of the different fixed states (constitutive activation or knockout) of various nodes. These in silico mutations are categorized into different response categories by their steady state percentage of closure and by their cumulative percentage of closure (CPC) values. The CPC is the summed fraction of simulations in which Closure = ON over the 50 timesteps-this is identical to the area under the curve in the simulation results plot on, for example, Fig 3. To classify the various responses to node constitutive activation or knockout, we used our previously defined categories [4]. To evaluate the variance in the simulation results of the unperturbed system, we performed 10 sets of simulations (where each set consisted of 4500 simulations over 50 time-steps) and determined the mean and standard deviation of the percentage of closure and CPC. These values inform the range of the close to wild type (or close to baseline) response category. In applying this criterion, we also considered the natural gaps in the results and the logical implications in the model.
In the simulations in the presence of ABA, hypersensitive refers to the case when the simulation reaches 100% closure faster (in fewer time steps) compared to the simulation of the unperturbed (wild type) system. Close to wild type is the response category when the CPC lies within two standard deviations of the average wild type CPC. Hyposensitivity refers to reaching 100% closure more slowly than the wild type system. Reduced sensitivity to ABA is the case when the percentage of closure stabilizes at a value much less than 100%, while insensitivity means 0% closure in all simulations over all time-steps. In the simulations in the absence of ABA, we refer to the scenario where no node is kept constitutively active as the baseline. As 17 nodes in the simulation start from a random state, the stable motif corresponding to Closure = ON may have a chance to stabilize, and so can lead to a small non-zero percentage of closure in the baseline case. Significantly increased response is the case when the percentage of closure reaches 100%. Slightly increased response is the case when the final percentage of closure is higher than the baseline but less than 100%; these values are generally in the 40-60% range. Close to baseline is the case when the percentage of closure and CPC are within two standard deviations of the respective baseline values. Decreased response is noted when the percentage of closure is less than the baseline. Significantly decreased response is the case when there is no closure, i.e., the CPC is 0.

Evaluation of the consistency between the simulated and experimentallyobserved effect of constitutive activation or knockout of a node
Hypersensitivity to ABA is considered consistent with experimental observation of an increased degree of stomatal closure (i.e., aperture decrease) compared to wild type at the same concentration of ABA. Close to wild type response to ABA in a simulation is considered consistent with stomatal closure that is not statistically significantly different from the response of the wild type at the same concentration of ABA. Because the variability of experimental conditions does not allow a ranking of reductions in closure response, hyposensitivity, reduced sensitivity and insensitivity to ABA are all considered consistent with experimental observation of decreased stomatal closure compared to wild type.
As the few experimental studies performed in the absence of ABA focus on stomatal closure, they do not allow a clear separation into more than two response categories. Accordingly, we group the response categories into two broad classes: closure (which includes slightly and significantly increased response compared to baseline) and absence of closure (which includes close to baseline response, decreased response and significantly decreased response compared to baseline, i.e., opening). Slightly or significantly increased response compared to baseline are considered consistent with experimental observation of stomatal closure (decreased aperture compared to the wild type) in the absence of ABA. Close to baseline response, decreased and significantly decreased response are considered consistent with a stomatal aperture measured in the absence of ABA that is not statistically significantly different from, or is greater than that of the wild type.

Stable motif, driver node, and attractor analysis
Stable motifs can be identified by finding strongly connected components (SCCs) in an expanded network that expresses the regulatory function of each node [19]. We used the software library available on https://github.com/jgtz/StableMotifs to determine the stable motifs of all the model versions and the attractors determined by these stable motifs. We verified the obtained attractors by simulations. We used the logic theory developed in [16] and the associated algorithms, available on https://github.com/parulm/suff_necc, to determine the internal and external drivers of stable motifs.
For the preparation of cell lysates, mesophyll protoplasts (2 ml, 4 x 105) were settled in 10 ml WI (0.5 M mannitol, 20 mM KCl, 4 mM MES, pH 5.7) solution in a 10 cm Petri dish for 30 min at room temperature. Protoplasts were collected by centrifugation at 100 X g for 2 min, and washed with 2 ml MMG (0.4 M mannitol, 15 mM MgCl 2 , 4 mM MES, pH 5.7) solution three times to remove extracellular calcium. To prepare the cell lysate for the PP2C assay, 100 μl of lysis buffer (20 mM Tris-HCl, pH 7.5, 20 mM KCl, 1 mM EDTA, 10 mM DTT, 0.5% Triton X-100, 50% glycerol, Roche protease inhibitor cocktail) was added to the protoplast pellet and vortexed for 5 sec. Lysis was allowed to proceed for 15 min on ice [41].
For each PP2C reaction, 2.5 μl of protoplast lysate was added in 22.5 μl reaction mix (50 mM Tris-HCl pH 7.5, 10 mM MgCl 2 , 5 mM DTT, 0.1 μg 32 P-casein with various CaCl 2 concentrations, or 1 μM Okadaic acid or 100 nM Calyculin A (to inhibit PP1 and PP2A), followed by incubation in a shaker at 30˚C, 600 rpm for 50 min. After the reaction, 2 ml cold 20% TCA was added, the reaction was kept on ice for 10 min, and then spun for 1 min in a Denville 260D microfuge at 12 Krpm and all supernatant, containing free 32 P (released from 32 P-casein by PP2C), was removed. The casein pellet was resuspended in 1 ml scintillation cocktail (UNI-VERSOL TM -ESLIQUID, MP Biomedicals) and each sample was counted for 3 min in a Beckman LS6500 multi-purpose scintillation Counter [41].
Supporting information S1 Table. The full names corresponding to the abbreviated node names in the ABA induced closure network, reproduced from [4] for convenience. (DOCX) S2 Table. List of cases where a source node that has a single successor (target) node is reduced and merged with its successor node. (DOCX) S3 Table. List of cases where the reduced source node has multiple target nodes.