Systems dynamic modeling of the stomatal guard cell predicts emergent behaviors in transport, signaling, and volume control.

The dynamics of stomatal movements and their consequences for photosynthesis and transpirational water loss have long been incorporated into mathematical models, but none have been developed from the bottom up that are widely applicable in predicting stomatal behavior at a cellular level. We previously established a systems dynamic model incorporating explicitly the wealth of biophysical and kinetic knowledge available for guard cell transport, signaling, and homeostasis. Here we describe the behavior of the model in response to experimentally documented changes in primary pump activities and malate (Mal) synthesis imposed over a diurnal cycle. We show that the model successfully recapitulates the cyclic variations in H+, K+, Cl−, and Mal concentrations in the cytosol and vacuole known for guard cells. It also yields a number of unexpected and counterintuitive outputs. Among these, we report a diurnal elevation in cytosolic-free Ca2+ concentration and an exchange of vacuolar Cl− with Mal, both of which find substantiation in the literature but had previously been suggested to require additional and complex levels of regulation. These findings highlight the true predictive power of the OnGuard model in providing a framework for systems analysis of stomatal guard cells, and they demonstrate the utility of the OnGuard software and HoTSig library in exploring fundamental problems in cellular physiology and homeostasis.

The guard cells, which surround stomatal pores in the epidermis of plant leaves, regulate the pore aperture to balance the often conflicting demands for CO 2 in photosynthesis with the need to conserve water by the plant. Stomatal transpiration accounts for much of the nearly 70% of global water usage associated with agriculture and has a profound impact on the water and carbon cycles of the world (Gedney et al., 2006;UNESCO, 2009). Recent studies have associated increases in continental water runoff with the rise in available CO 2 and decreases in stomatal transpiration (Gedney et al., 2006) and have suggested that stomatal behavior skews the impact of greenhouse gasses on fresh water resources (Betts et al., 2007). The past half century has generated a vast wealth of knowledge for guard cell transport, signaling, and homeostasis, resolving the properties of all of the major transporters and many of the signaling pathways that control them (Blatt, 2000a;Schroeder et al., 2001;Blatt et al., 2007;Wang and Song, 2008;McAinsh and Pittman, 2009). Even so, resolving many aspects of stomatal dynamics remains a challenge. These studies have yet to yield any detail about how the entire network of transporters works as a unit to modulate solute flux and regulate stomatal aperture. Quantitative systems analysis offers one approach to this problem that is now much needed. Efforts to model stomatal function to date generally have been driven by a top-down approach: The mechanics of stomatal movements are subsumed within a few empirical parameters of linear hydraulic pathways and conductances (Farquhar and Wong, 1984;Ball, 1987;Williams et al., 1996;Eamus and Shanahan, 2002;West et al., 2005). These models have proven useful at the plant and community levels; but they have not incorporated the essential detail to support an understanding of the molecular and cellular mechanics that drive stomatal movements.
In the previous article (Hills et al., 2012) we introduced a computational approach to developing a dynamic model of the stomatal guard cell based on the HoTSig library and OnGuard software. We resolved an OnGuard model that takes account of all of the fundamental properties for transporters at the plasma membrane and tonoplast, the salient features of osmolite metabolism, and key homeostatic and dynamic signaling characteristics that have been described in the literature. The model successfully integrated a number of the steady-state characteristics of guard cells, recapitulating the patterns in guard cell response to the extracellular variables of KCl and CaCl 2 concentrations and to extracellular pH. Here we explore the capacity of the model to reproduce diurnal oscillations in guard cell membrane transport and malate (Mal) metabolism, and its consequences for the dynamics of guard cell volume, turgor pressure, and stomatal aperture. We demonstrate the true predictive power of the OnGuard model in generating a number of unexpected and counterintuitive outputs. Among these, the model yields counterintuitive changes in cytosolic-free [Ca 2+ ] ([Ca 2+ ] i ) and a daily exchange of Cl 2 with Mal that are well documented in the literature, but have been suggested to require additional and complex levels of regulation. These behaviors are accounted for entirely by the known kinetic features of the transporters encoded in the model. Thus, the results demonstrate the predictive power of the OnGuard model as a framework from which to test the basic tenets of the stomatal behavior and to explore the interactions of transport and metabolism in the guard cell system.

The Diurnal Model
By nature guard cells define a closed cellular system within the surrounding apoplastic volume of the leaf tissue (Wille and Lucas, 1984). Transport of K + and Cl 2 , both major contributors to the osmotic content of guard cells, and of H + and Ca 2+ must take place across the plasma membrane and coordinate with the metabolism of the predominant, osmotically active organic compounds, notably of Suc and Mal. A number of excellent reviews summarize the very large body of experimental data that pertain to transport, metabolism, and their regulation in guard cells Willmer and Fricker, 1996;MacRobbie, 1997;Thiel and Wolf, 1997;Blatt, 2000b;Hetherington, 2001;Schroeder et al., 2001;Outlaw, 2003;Dreyer et al., 2004;Blatt et al., 2007;Pandey et al., 2007). Details of the kinetic parameters incorporated for each of the transporters and for Mal and Suc metabolism are provided in the preceding article (Hills et al., 2012). We draw attention here to a few key conceptual points that bear on dynamic behavior of the guard cell.
First and foremost, electrophysiological and radiotracer flux studies (Willmer and Fricker, 1996;MacRobbie, 1997;Blatt, 2000b;Schroeder et al., 2001;Dreyer et al., 2004) have shown that the dominant membrane transporters are active under most physiological conditions, although their operation may be kinetically limited by substrate (ion) availability and/or membrane voltage. This observation contrasts with a common misconception that transporters either activate or shut down fully in response to various stimuli and, thus, may be modeled kinetically as simple on-off or Boolean processes (Li et al., 2006). For example, currents through both inward-(I K,in ) and outward-rectifying (I K,out ) K + channels are clearly evident under voltage clamp both in the absence and presence of the water-stress hormone abscisic acid (ABA) that triggers stomatal closure (Blatt and Armstrong, 1993;Romano et al., 2000;Garcia-Mata et al., 2003). [A list of abbreviations can be found in supplemental appendix S5 of this article's companion (Hills et al. 2012).] ABA reduces the activity of I K,in and promotes that of I K,out through qualitative changes in gating, although membrane depolarization is an overriding factor in determining the K + flux (Thiel et al., 1992;Romano et al., 2000). Thus, one conclusion to be drawn from these observations is that guard cells are unlikely to be found in a state in which the net flux of an osmotically active solute is zero. Instead, the observations suggest that guard cells transit between situations of osmotic solute uptake and loss, to achieve by time-averaged approximation a dynamic range of (quasi-)steady states in solute contents and, consequently, of stomatal apertures.
This conclusion resonates with a second set of observations that guard cell membrane voltage can oscillate with periods from a few tens of seconds to many minutes (Thiel et al., 1992;Gradmann et al., 1993;Grabov and Blatt, 1998) between two quasi-stable states, a depolarized state associated with K + and Cl 2 efflux, and a hyperpolarized state that can be shown to drive K + uptake Clint and Blatt, 1989;Blatt and Thiel, 1993;Blatt, 2000b). Such oscillations have been observed to occur spontaneously; they can be induced by changes in external solute (e.g. KCl) concentration as well as hormone treatments, including ABA and auxin (Thiel et al., 1992;Blatt and Armstrong ,1993;Blatt and Thiel, 1994;Sokolovski and Blatt, 2007); and they have been associated with substantial changes in [Ca 2+ ] i (Grabov and Blatt, 1998;Blatt, 2000b). These characteristics lead to an expectation of considerable temporal dynamism in membrane transport and its regulation in guard cells.
Finally, a third point that helped inform our approach to modeling arises from the general finding that the ion fluxes facilitating stomatal movements reflect a small fraction only of the potential activities of several transporters, notably of the predominant voltage-activated K + channels and Ca 2+ channel at the plasma membrane (compare with Thiel et al., 1992Thiel et al., , 1993Blatt et al., 1990bBlatt et al., , 1992Hamilton et al., 2000), the CLC H + -Cl 2 antiporter, the SV (TPC) and FV cation channels, and the VCL and VMAL (ALMT) anion channels at the tonoplast (compare with Pottosin and Schönknecht, 2007;DeAngeli et al., 2009). In general, the activities of these transporters are kinetically limited by their inherent gating properties within the range of voltages typically found across the plasma membrane and tonoplast. Thus, the channel activities recorded under voltage clamp at the voltage extremes do not represent the typical current amplitudes in vivo, even if they are useful as measures of the biophysical and regulatory properties associated with the currents.

The Reference Cycle
Formulating a dynamic model of cellular homeostasis normally begins with the definition of a Reference State representing the resting physiological condition of the system under study. However, when the model describes periodic cellular responses, it is convenient to circumscribe the associated parameters within a Reference Cycle. Previously (Hills et al., 2012), we introduced paired Reference States corresponding to the guard cell of the closed and open stomata and reflecting the pattern of nocturnal closure and diurnal aperture found in most plant species (Willmer and Fricker, 1996). We envisaged the Reference Cycle as a logical extension of these paired states, with the nocturnal closed state equating with our earlier Closed Reference State (Hills et al., 2012) in which the guard cells retained a baseline of osmotic load and a minimum of ion flux across the tonoplast and plasma membrane. To define the diurnal cycle, the turnover rates of all primary pumps (H + -ATPases, Ca 2+ -ATPases, H + -PPase) across tonoplast and plasma membrane were assigned values in the nocturnal (closed) state of 5%, 10%, or 20% of their maximum output consistent with experimental estimates of the known light-stimulated activities (see Assmann et al., 1985;Gotow et al., 1985;Goh et al., 1995Goh et al., , 1996Kinoshita et al., 2001;Shimazaki et al., 2007). The rates of Suc and Mal synthesis were set to zero in the dark.
Transition to the light was introduced by incorporating a hyperbolic dependence on light for these components. For simplicity we made use of a common L 1/2 value of 50 mmol m 22 s 21 for the fluence dependence of all of the light-dependent processes and we simulated the light cycle as a 6-h ramp from 0 to 2,200 mmol m 22 s 21 , a further 6-h ramp down to 0 mmol m 22 s 21 , and a 12-h dark period. All other model parameters were kept constant in any one simulation. As described in the following paragraphs, one consequence of using a simple ramp cycle was to delay the final steady state achieved in daytime outputs until late in the daylight period. A steady state in these outputs was achieved much earlier in the daylight period when a Gaussian-shaped light regime was used (not shown, available to the user in OnGuard). We chose nonetheless to use the ramp cycle here because it yields simple hyperbolic functions in the activation/deactivation of the pumps and metabolism in relation to elapsed time and, as a result, it is much easier to appreciate the interactions between transporters. The parameter ensemble for the osmotic composition of the guard cell, buffering constants for H + and Ca 2+ , the relevant transporters, Suc and Mal metabolism and their kinetic descriptors are those described in the preceding article (Hills et al., 2012), and are based on published data for Vicia. This final model was resolved through fine adjustments to the ensemble, introduced between simulations. The process was constrained throughout by experimental data, including electrophysiological data on single channels and pumps in isolated guard cells, data from measurements of stomatal apertures and calculated turgor pressures in populations of stomata, and from metabolic studies on Mal and Suc metabolism (Dittrich and Raschke, 1977;Talbott and Zeiger, 1993;Willmer and Fricker, 1996;Hills et al., 2012). Tabulations of sample adjustments and effects on simulation outputs are summarized in figures 4 and 5 of the preceding article (Hills et al., 2012).

Diurnal Changes in Macroscopic Outputs
In vivo (Raschke, 1979;MacRobbie, 1987;Willmer and Fricker, 1996), the diurnal cycle is associated with a progressive accumulation of K + with Cl 2 , and subsequently with Mal, leading to an increase in turgor pressure and stomatal aperture in the first half of the day; only later do lengthening of the spectral wavelength and other factors lead to a decline in H + -ATPase and other transport activities. Mal synthesis and accumulation delays the decline in osmotic load, often effectively exchanging with Cl 2 during the daylight period. Daylight events are thus accompanied by large solute and water fluxes in and out of guard cells, mostly between vacuole and apoplast passing through the cytosolic compartment of the guard cells, and by an accelerated metabolic activity to generate and degrade additional osmolites, all with a highly regulated time course. Figures 1 to 7 summarize the principle outputs of simulation with the OnGuard model and its capacity to reconstruct these, and related characteristics of the guard cell diurnal cycle.
As anticipated, modulating the activities of the primary ion pumps, Suc and malic acid synthesis resulted in a cycle of diurnal stomatal opening and nocturnal closure. As shown in Figure 1, stomatal apertures varied over a physiological range between roughly 4 and 13 mm and were paralleled by physiologically reasonable changes in guard cell volume, turgor, and vacuolar volume percent (compare with Willmer and Fricker, 1996), the latter defined as the percentage of the cell volume occupied by the vacuole. Stomatal closure at the end of the day was followed by a small and gradual rise in aperture and the associated macroscopic outputs, in effect anticipating the start of the next day much as has been observed in vivo (Gorton et al., 1993;Meidner and Willmer, 1993). In the simulation conditions-10 mM KCl and 1 mM CaCl 2 , pH 6.5, typical of many studies with epidermal peels-the onset of daylight was associated with hyperpolarization of the plasma membrane to voltages near 2130 mV and the dark period was accompanied by depolarization of the plasma membrane to voltages near the equilibrium voltage for K + , consistent with the diurnal cycle in energetic outputs of the ATP-driven pumps (see Raschke, 1979;Spanswick, 1981;Blatt, 1987b;Blatt and Clint, 1989;Clint and Blatt, 1989;McClure et al., 1989;Kinoshita et al., 1995 and below). The tonoplast showed much smaller diurnal variations and an inverse pattern, with voltages of 220 mV in daylight and near 250 mV at night, indicating a dominance of secondary conductances especially during the daylight hours (Goldsmith and Goldsmith, 1978;Gobert et al., 2007). Transition to the dark was marked by a period of voltage excursions, or action potentials, at the plasma membrane, and was followed at first by a period during which the membrane rose to voltages near 240 mV. We return to these observations later.
Osmotic Contents and Flux of K + Figure 2 shows that stomatal opening was accompanied by increases in K + concentration in the cytosol and in the vacuole, these increases following monophasic kinetics, and decreased again in the first hours of dark during stomatal closure, consistent with experimental measurements (see Hills et al., 2012). The K + concentrations in both compartments gradually rose again later at night in parallel with stomatal aperture (Fig. 1). In the cytosol, the K + concentration varied between approximately 110 mM in the closed state and 210 mM in the open state; in the vacuole, K + concentrations ranged between approximately 20 and 120 mM ( Fig. 2A). Analysis of the total K + flux (Fig. 2B) showed that the larger proportion of K + influx across the plasma membrane was shunted across the tonoplast to the vacuole during the day and this pattern reversed in the first hours of dark, as expected to accommodate the large changes in vacuolar volume (Fig. 1). (Note that here, and in all subsequent analyses, we have defined a positive flux as movement of the ionic species [not the charge] out of the cytosol, either across the plasma membrane or the tonoplast.) At the plasma membrane ( Fig. 2C), K + influx was dominated by I K,in in the first half of the day, this flux relaxing to approximately 20% of its maximum, roughly equivalent to that through the H + -K + symport in the second half of the day. Closure was marked by the predominance of K + efflux Figure 1. Macroscopic outputs from the OnGuard model resolved over the diurnal Reference Cycle (see text) with the standard environmental parameters of 10 mM KCl, 1 mM CaCl 2 , and pH 6.5. Data represent a 5-d window of the stable Reference Cycle (12 h light:12 h dark, indicated by bars above). OnGuard model parameters are those described in the preceding article (Hills et al., 2012). Shown are plasma membrane and tonoplast voltage (A), stomatal aperture, turgor pressure, and total guard cell volume (B), and the percentage of cell volume occupied by the vacuole (C). Figure 2. K + contents and analysis of K + fluxes at the plasma membrane and tonoplast resolved over the diurnal Reference Cycle as described in Figure 1 (12 h light:12 h dark, indicated by bars above). Shown are cytosolic and vacuolar [K + ] (A), the net K + flux across the plasma membrane and tonoplast (B), the K + flux through the K +permeable transporters at the plasma membrane, comprising the two K + channels and the H + -K + symporter (C), and the K + flux through the K + -permeable transporters at the tonoplast, comprising the TPK and FV channels (D). K + flux through the TPC channel accounted for less than 1% of either of the other channel fluxes, and has therefore been omitted for purposes of clarity. Note that positive flux is defined as movement of the ionic species (not the charge) out of the cytosol, either across the plasma membrane or the tonoplast.
through I K,out , which relaxed to a near-zero value before the second half of the night. At the tonoplast (Fig.  2D), K + flux in both directions was mediated largely by the TPK channel with approximately 10% passing through the FV channel, consistent with experimental evidence for the importance of TPK1 for K + homeostasis (Gobert et al., 2007). Less than 1% of the K + flux passed through the TPC channel (not shown), consistent with evidence that in Arabidopsis (Arabidopsis thaliana) the tpc1 mutant has little effect on tissue K + contents (Peiter et al., 2005). These fluxes largely mirrored the diurnal pattern of those at the plasma membrane, with the exception that both influx and efflux were observed, primarily through the TPK channel, later during stomatal closure as the channel provided the major flux pathway for charge balance across the tonoplast (see Figs. 3, 4, and 6).

Suc and Malic Acid Metabolism and the Mal Osmoticum
Mal is known to be a major contributor to the osmotic contents of the guard cell, notably in the vacuole, and may accumulate to concentrations in excess of 150 mM during the day in some circumstances (Van Kirk and Raschke, 1978a;Outlaw, 1990;Talbott and Zeiger, 1993;Willmer and Fricker, 1996). For modeling purposes, we neglected Suc transport, instead using Mal as a proxy for Suc as an osmoticum (Hills et al., 2012), consistent with the predominance of Mal in many circumstances (Willmer and Fricker, 1996). Thus, with Suc synthesis activated in the OnGuard model in the daylight, malic acid production was engaged and then relaxed to a near-steady-state rate of production near 12 amol s 21 , sufficient to accumulate over 40 mM h 21 in the cytosol (Fig. 3A); in the first 2 to 3 h of dark, a fraction of this Mal was metabolized to Suc and broken down as it was recovered from vacuolar stores (see below, and Van Kirk and Raschke, 1978a;Outlaw, 1990;Talbott and Zeiger, 1993). During the day, the bulk of Mal production was diverted by transport of Mal 22 across the tonoplast, leading to a rise in total vacuolar Mal concentration from approximately 20 mM in the closed state to over 100 mM in the open state; in parallel, Mal in the cytosol ranged from 0.7 to values near 25 mM, respectively (Fig. 3, B and C), much as has been estimated from experimental data (see Hills et al., 2012). Mal accumulation during the day was accompanied by a small increase in export across the plasma membrane, culminating in a substantial efflux from the vacuole and cytosol during closure in the first hours of dark ( Fig. 3, C-E). Again, these characteristics match well the experimentally documented evidence for Mal flux (Allaway, 1973;Raschke, 1977, 1978b;Gotow et al., 1985;Tarczynski and Outlaw, 1990). At the plasma membrane, the OnGuard model predicted Mal efflux to be mediated largely by the SLAC-type anion channel during daylight and through the R-(ALMT-)type anion channel during stomatal closure in the first hours of dark (Fig. 3C), much as has been deduced from experimental data (Hedrich et al., 1990;Schmidt and Schroeder, 1994;Dietrich and Hedrich, 1998;Negi et al., 2008;Vahisalu et al., 2008). Mal flux in both directions across the tonoplast was dominated by the VMAL (ALMT) anion channel, as previously postulated (Meyer et al., 2011), the flux thus driven by the balance of electrochemical driving forces across the tonoplast and plasma membrane, and the rapid decline in cytosolic Mal early during stomatal closure (see Fig. 3B).

Osmotic Contributions from Cl 2
Both cytosolic and vacuolar Cl 2 concentrations exhibited more complex, biphasic responses to the diurnal cycle in the OnGuard model ( Fig. 4A), consistent with experimental observation (Raschke and Schnabl, 1978;Zeiger, 1993, 1996). In the cytosol, the Cl 2 concentration rose during the first hours of daylight to a maximum near 15 mM, declining thereafter to below 10 mM at the end of the day. This pattern was repeated in reverse during the night: the cytosolic Cl 2 concentration falling in the first 2 h of dark, thereafter rising from a minimum of 2 to 6 mM at the end of the dark period. Vacuolar Cl 2 concentration, by contrast, decreased monotonically from 45 to less than 20 mM during daylight, declining further to near 15 mM in the first 2 to 3 h of dark, before rising again in parallel with the cytosolic Cl 2 concentration. The dark rise in both vacuolar and cytosolic Cl 2 concentrations was suppressed when the external Cl 2 concentration was reduced from the standard 12 mM (10 mM KCl + 1 mM CaCl 2 ) and the decline in cytosolic and vacuolar Cl 2 concentrations during the second half of the daylight period was absent at external Cl 2 concentrations of 22 and 32 mM (not shown).
Analysis of the Cl 2 flux (Fig. 4B) showed a net influx across the plasma membrane that gradually declined during the dark period, and a persistent but smaller efflux throughout much of the day. The first 2 to 3 h of dark were marked by excursions in Cl 2 flux between efflux and influx across the plasma membrane. Across the tonoplast, the inverse pattern was observed: A small Cl 2 influx occurred throughout the daylight, the first 2 to 3 h of dark were marked by large excursions in Cl 2 influx, and during the remainder of the dark Cl 2 was transported outward (into the vacuole). Thus, stomatal opening was accompanied by a net flux of Cl 2 from the vacuole to the apoplast; closure at the start of the dark period was marked by much larger fluxes of Cl 2 from the vacuole to the cytosol and export across the plasma membrane; and this pattern reversed after the first 2 to 3 h of dark. The rise in cytosolic Cl 2 concentration during the first hours of the day arose from the rapid Cl 2 influx across the tonoplast and a slower rise in the rate of Cl 2 export across the plasma membrane.
At the plasma membrane (Fig. 4C), Cl 2 influx was mediated by the H + -Cl 2 symport throughout the daylight period and to a lesser extent during the night. The decline and reversal in net flux arose as the consequence of Cl 2 efflux through the SLAC channel during the day and the first 2 to 3 h of the dark period. Cl 2 export was augmented by efflux through the R-(ALMT-)type anion channel in the first 2 to 3 h of dark. By contrast, Cl 2 flux at the tonoplast was dominated throughout the diurnal cycle by the VCL channel, the CLC H + -Cl 2 antiport making a lesser contribution to Cl 2 accumulation in the vacuole throughout the 24-h cycle (Fig. 4D). As a consequence, reversal of the net flux of Cl 2 across the tonoplast was driven by changes in the electrochemical driving force for Cl 2 and a passive shift in VCL flux, which followed the decline in Cl 2 concentration in the cytosol and changes in vacuolar membrane voltage (Figs. 1 and 4A). One conclusion to be drawn from this analysis is that net fluxes of Cl 2 and Mal 22 are not linked. It is certainly the case that cytosolic Mal affects both the major anion channels at the plasma membrane, but the effect is a moderate suppression in channel activity (Wang and Blatt, 2011), which could be expected to enhance Cl 2 retention rather than enhancing its efflux. Thus, the OnGuard model successfully reproduced the apparent daytime exchange between Cl 2 and Mal (Van Kirk and Raschke and Schnabl, 1978), but without a requirement for direct control of Mal on Cl 2 flux as previously postulated (Hedrich and Marten, 1993). . Chloride contents and analysis of Cl 2 fluxes at the plasma membrane and tonoplast resolved over the diurnal Reference Cycle as described in Figure 1 (12 h light:12 h dark, indicated by bars above). Shown are total cytosolic and vacuolar [Cl 2 ] (A), the net flux of Cl 2 across the plasma membrane and tonoplast (B), the flux of Cl 2 through the Cl 2 -permeable transporters at the plasma membrane, comprising the SLAC and R-(ALMT-)type anion channels and H + -Cl 2 symporter (C), and the flux of Cl 2 through the Cl 2 -permeable transporters at the tonoplast, comprising the VCL channel and CLC H + -Cl 2 antiporter (D). Note the difference in scales between C and D. Note that positive flux is defined as movement of the ionic species (not the charge) out of the cytosol, either across the plasma membrane or the tonoplast.

Cytosolic pH, the H + -ATPase, and H + -Coupled Transport
Direct measurements in guard cells (Irving et al., 1992;Blatt and Armstrong, 1993;Thiel et al., 1993;Armstrong et al., 1995;Grabov and Blatt, 1997;Zhang et al., 2001) have indicated resting pH i values near 7.5 to 7.7, with one notable exception , and vacuolar pH generally has been estimated to situate between pH 4.5 and 6 (see Hills et al., 2012). Figure 5A shows that the OnGuard model faithfully reproduced the characteristics anticipated for the pH, both in the cytosol (pH c ) and vacuole (pH v ). Values for pH c remained close to 7.6 to 7.7 throughout much of the day, dropped to 7.5 during the night, and rose again above 7.6 before the start of the subsequent daylight period. This daily cycle was punctuated by an acid-going transition in the first 1 to 2 h of daylight and a similar period at the end of the day when the model showed a 0.1-unit pH i rise to a value near 7.8. Details of the diurnal changes in pH i in the guard cell remain to be explored in vivo. However, a 0.1 to 0.3 unit rise in pH i is known to occur during the first 10 to 30 min of stomatal closure evoked by ABA (Irving et al., 1992;Blatt and Armstrong, 1993;Zhang et al., 2001), so it is significant that a similar rise in pH c coincided with the initial and very rapid fall in stomatal aperture (compare Figs. 1 and 5).
Direct measurements of pH i indicate substantial static buffering, as incorporated in the OnGuard model (Hills et al., 2012), and a capacity for dynamic pH i control (Grabov and Blatt, 1997). However, static buffering will accomodate finite changes in H + concentration and cannot account for long-term H + loads such as imposed by metabolism. At 1 to 3 pmol h 21 (guard cell) 21 , maximum rates of malic acid synthesis alone in Vicia guard cells are sufficient to generate a H + load near 0.5 M h 21 (Gotow et al., 1985;Outlaw, 1990;Tarczynski and Outlaw, 1990). Thus, H + elimination from the guard cell is essential for pH i homeostasis. As expected, the OnGuard model predicted the vast bulk of the H + production associated with daytime Mal synthesis (Fig. 3A) to be exported via the plasma membrane H + -ATPase, with roughly 20% transported to the vacuole (Fig. 5B). The model returned a net H + efflux that reached a maximum close to 30 amol s 21 in the first 2 h of the day and settled to a rate near 20 amol s 21 in the second half of the daylight period before declining to near zero, consistent with the diurnal pattern of H + flux in vivo (Willmer and Fricker, 1996). A substantially higher H + flux was estimated through the H + -ATPase than was generated as net H + export, much as can be concluded from its activity in vivo (Blatt, 1987a(Blatt, , 1988Blatt and Clint, 1989;Fricker and Willmer, 1990;Lohse and Hedrich, 1992), such that 30% to 40% of the H + export via the H + -ATPase was balanced by H + entry coupled with Cl 2 uptake and, to a lesser extent, with K + uptake (Fig. 5C). One prediction to come from the OnGuard analysis, therefore, is that pH i homeostasis is likely to be affected by the absence of inorganic ions available for symport with H + , and this conclusion accords with analogous considerations in several plant cells (Smith, 1973;Sanders et al., 1989;Guern et al., 1991) and in fungi (Blatt and Slayman, 1987).

Vacuolar pH, H + , and Mal Transport
Malic acid comprises the major pH buffer in the vacuole and its accumulation is associated with acidification of the vacuolar contents (Van Kirk and Raschke, 1978a;Outlaw, 1990;Talbott and Zeiger, 1993;Willmer and Fricker, 1996). In the OnGuard model these diurnal changes resulted in pH v oscillating between 4.6 units during the day and around 5.0 to 5.4 units during the second half of the night (Fig.  5A). All evidence indicates that the organic acid is Figure 5. Cytosolic and vacuolar pH, and analysis of H + fluxes across the plasma membrane and tonoplast resolved over the diurnal Reference Cycle as described in Figure 1 (12 h light:12 h dark, indicated by bars above). Shown are cytosolic and vacuolar pH, pH c , and pH v , respectively (A), the net H + flux across the plasma membrane and tonoplast (B), the H + flux through the H + -permeable transporters at the plasma membrane, comprising the H + -ATPase, and the H + -K + and H + -Cl 2 symporters (C), and the H + flux through the H + -permeable transporters at the tonoplast, comprising the VH + -ATPase, VH + -PPase, the CLC H + -Cl 2 antiporter, and the CAX H + -Ca 2+ antiporter (D). Note that positive flux is defined as movement of the ionic species (not the charge) out of the cytosol, either across the plasma membrane or the tonoplast.
transported as the fully deprotonated (Mal 22 ) formwith the VMAL channel as the primary pathway for tonoplast Mal 22 flux (see Fig. 3D; for review, see Hills et al., 2012)-thereby implying that a component of charge balance is achieved with parallel H + transport via the tonoplast VH + -ATPase and H + -PPase. In the OnGuard model (Fig. 5D), H + transport to the vacuole was roughly divided between the VH + -ATPase and H + -PPase, as estimated from experimental data in several plant species and cell types (Rea and Poole, 1993;Martinoia et al., 2007). The model predicted H + reentry to the cytosol to be dominated by its exchange with Ca 2+ mediated by the CAX H + -Ca 2+ antiport, primarily during stomatal closure in the first 2 to 3 h of dark; H + exchange with Cl 2 was indicated to contribute a smaller fraction to H + return from the vacuole, at least within the context of the extracellular parameters of the simulation (above). In total, however, the H + flux represents a small fraction of the total H + load estimated for the guard cell: A comparison of fluxes and ions accumulated in the simulation indicates that net H + flux through the combined pathways of the VH + -ATPase and H + -PPase contributed approximately 30% to vacuolar charge balance with the remainder assumed by K + flux (see Fig. 2).
One underlying assumption in these estimations is of the diurnal variations in energy-dependent H + transport. To validate this assumption, we turned to Arabidopsis, carrying out quantitative PCR for the VHA-a2, VHA-a3, A, and C subunits of the VH + -ATPase on mRNA isolated at 4-h time intervals from leaves; additionally, we isolated microsomal membranes from the same leaf samples to assay for Concanamycin-A-sensitive ATPase activity (Palmgren, 1990;Brüx et al., 2008). The results of these experiments (Supplemental Figs. S1 and S2 and Supplemental Materials and Methods S1), demonstrated that both the transcription and enzymatic activity of the VH + -ATPase pass through a diurnal cycle, with a maxima during midday and a minimum close to 10% of this value during the nocturnal period. Furthermore, analysis of the model outputs showed that transport of the major osmotica and changes in pH v followed the pattern consistent with experimental observation (above). Thus, our prediction of a diurnal cycle in VH + -ATPase activity as built into the model are borne out, and the results gave us confidence to assume a comparable cycle of activity for the H + -PPase. Figure 6A shows that stomatal opening was accompanied by an increase in total [Ca 2+ ] in the vacuole, rising from a mean value near 15 mM in the dark to 25 mM at the end of the daylight period; total [Ca 2+ ] in the cytosol remained between 0.2 and 0.4 mM throughout much of the daylight period. [Ca 2+ ] i rose from a resting value near 180 nM in the dark to a quasi-steady state close to 300 nM in the second half of the daylight period (Fig. 6B). In the final hour of daylight and the first 2 to 3 h of dark, the OnGuard model generated a series of voltage and [Ca 2+ ] i excursions that culminated with stomatal aperture relaxing to a minimum closed value and recovery of dark levels in total [Ca 2+ ] and [Ca 2+ ] i . These characteristics are broadly consistent with a number of previous observations: They mirror the action-potential-like oscillations in membrane voltage and [Ca 2+ ] i elevations that have been associated with stomatal closure (Gilroy et al., 1991;Irving et al., 1992;McAinsh et al., 1992;Thiel et al., 1992;Blatt and Armstrong, 1993;Gradmann et al., 1993; Grabov and , the net flux of Ca 2+ across the plasma membrane and tonoplast (C), the Ca 2+ flux through the Ca 2+ -permeable transporters at the plasma membrane, comprising the hyperpolarization-activated Ca 2+ channel and the Ca 2+ -ATPase (D), and the flux of Ca 2+ through the Ca 2+ -permeable transporters at the tonoplast, comprising the Ca 2+ -ATPase, the CAX H + -Ca 2+ antiporter, and the TonVCa Ca 2+ channel (E). Note the difference in scales between D and E. Note that positive flux is defined as movement of the ionic species (not the charge) out of the cytosol, either across the plasma membrane or the tonoplast. Blatt, 1998;Staxen et al., 1999), and they recapitulate diurnal variations in [Ca 2+ ] i with resting values elevated in the daytime relative to the night (Dodd et al., 2005(Dodd et al., , 2006.

Diurnal Variations in Cytosolic-Free [Ca 2+ ]
Detailed analysis of the OnGuard output (Fig. 6C) showed a net Ca 2+ influx across the plasma membrane and export to the vacuole throughout the daylight period, although the bulk of this flux occurred in the first 8 h of the day. The flux direction was reversed during the first 2 to 3 h of dark, albeit with excursions in tonoplast flux to positive (outward, directed to the vacuole) values, before both membrane fluxes relaxed to near-zero values. At the plasma membrane (Fig. 6D), this Ca 2+ flux was dominated by the inward-rectifying Ca 2+ channel during daylight hours and by the Ca 2+ -ATPase in the first hours of dark. Much of Ca 2+ exported to the vacuole (Fig. 6E) by the vacuolar Ca 2+ -ATPase was returned through the TonVCa Ca 2+ channel during daylight hours. Activity of the CAX H + -Ca 2+ antiport was evident in the first 2 to 3 h of dark, and all three transporters underwent a series of excursions, albeit with Ca 2+ return via the Ca 2+ channel dominating, before the component Ca 2+ fluxes relaxed to near-zero values in the second half of the dark period. Comparison of the flux amplitudes predicted by the OnGuard model for the plasma membrane and tonoplast shows that vacuolar Ca 2+ transport dominated by at least one order of magnitude throughout the diurnal cycle. This prediction is generally in accord with the major roles for endomembrane sequestration and release in Ca 2+ homeostasis (Sanders et al., 2002;Blatt et al., 2007;McAinsh and Pittman, 2009) and with direct evidence in guard cells for its importance in potentiating and shaping [Ca 2+ ] i signals (Gilroy et al., 1991;McAinsh et al., 1991;Grabov and Blatt, 1997Garcia-Mata et al., 2003).
Of the model outputs, the diurnal variation predicted for resting [Ca 2+ ] i is provocatively counterintuitive. The OnGuard model incorporated a 5-fold increase in the activities of the Ca 2+ -ATPases during daylight hours at both the plasma membrane and the tonoplast. Yet resting [Ca 2+ ] i rose substantially in the light, despite the enhanced Ca 2+ export from the cytosol. How can we explain such behavior? The OnGuard model outputs yield a simple, if unexpected answer. Figure 7A illustrates the underlying characteristics for Ca 2+ transport at the plasma membrane, cross-referenced to membrane voltage and Ca 2+ flux (Fig. 7, C and D), with data captured from screenshots at time points throughout the diurnal cycle (curves for other transporters are not shown for clarity, see also Hills et al., 2012). The IV curves for the Ca 2+ channel and Ca 2+ -ATPase shown here encapsulate the kinetic descriptors for each of these transporters as a function of the membrane voltage, and they correspond to the component IV curves that would be recorded under voltage clamp. Comparing the currents at the resting (free-running) membrane voltages, indicated by the voltage axis intercept of the total membrane IV curve, it is evident that the balance of Ca 2+ transport led to a net Ca 2+ influx into the cytosol early in the day, a gradual decay in this flux toward the end of the day, and a greater (and oscillating) Ca 2+ efflux early in the night before the net Ca 2+ flux relaxed back to a value near zero. Thus, at the plasma membrane, stimulating the pumps and the resulting membrane hyperpolarization during the day affects directly the kinetic limits for Ca 2+ export by the Ca 2+ -ATPase, suppressing its ability to export Ca 2+ and at the same time promoting Ca 2+ entry through the hyperpolarization-activated Ca 2+ channels.
Analysis of current at the tonoplast (Fig. 7B) shows a complementary pattern of Ca 2+ flux. In this case, enhancing Ca 2+ -ATPase activity and the small membrane depolarization during the day can be seen to have promoted Ca 2+ efflux to the vacuole sufficient to keep pace with Ca 2+ entry from the apoplast. Additionally, as [Ca 2+ ] i rose during the day it enhanced Ca 2+ release, which largely balanced efflux via the tonoplast Ca 2+ -ATPase. Finally, the decline in pump activities at the end of the day left the predominant flux for Ca 2+ directed inward across the tonoplast over the same time period that depolarization of the plasma membrane favored Ca 2+ export through the Ca 2+ -ATPase (Fig. 7, C and D). Thus, stimulating the ATP-and pyrophosphate-dependent pumps in the OnGuard model, including the Ca 2+ -ATPases at the plasma membrane and tonoplast, translates to a net influx of Ca 2+ from the apoplast, through the cytosol, to the vacuole during the day; this net flux reverses as the stomata closed during the early hours of the night. In short, the counterintuitive changes in [Ca 2+ ] i during the day resulted from membrane hyperpolarization and consequent alterations in kinetic restrictions on the Ca 2+ -ATPase and Ca 2+ channel at the plasma membrane; these changes were paralleled by a progressive decline in the capacity for Ca 2+ sequestration. It is worth noting, too, that the predicted variations in Ca 2+ flux and [Ca 2+ ] i were achieved without encoding any extrinsic feedback such as has been suggested to affect the diurnal capacity for endomembrane Ca 2+ release (Dodd et al., 2007). Of course, the simulation does not preclude this or other feedback regulation, but it demonstrates the ability of the model to encapsulate otherwise unexpected biological phenomena within the emergent properties of a relatively simple system. As a corollary, it also suggests more subtle roles for control mechanisms such as posttranslational regulation affecting the capacity for Ca 2+ release.

Predicting Endomembrane Ca 2+ Flux and [Ca 2+ ] i Oscillations
One requirement for the OnGuard model was our inclusion of a self-limiting Ca 2+ flux pathway for Ca 2+ release from the vacuole. Such pathways are essential for evoked Ca 2+ release and, like the inositol trisphosphate (IP 3 )-and ryanodine-receptor channels of animals (Bezprozvanny et al., 1991;Hille, 2001), must incorporate an element of self inhibition to prevent prolonged [Ca 2+ ] i elevation and exhaustion of the Ca 2+ store. Ca 2+ channels that are activated by Ca 2+ and ligands such as IP 3 , cADP ribose, and nitric oxide are known or have been demonstrated indirectly to reside in the endomembranes and vacuoles of plant cells, including those of Vicia guard cells (Alexandre et al., 1990;Leckie et al., 1998;Garcia-Mata et al., 2003). With the exception of the plasma membrane Ca 2+ channel of guard cells (Hamilton et al., 2000), self inhibition has not been demonstrated in plants per se; nevertheless, it is implicit in evoked [Ca 2+ ] i oscillations such as those observed in Vicia guard cells (McAinsh et al., 1995;Grabov and Blatt, 1998;Staxen et al., 1999;Garcia-Mata et al., 2003;Sokolovski et al., 2005). The SV (TPC1) channel defines a Ca 2+permeable pathway in the guard cell tonoplast with gating characteristics well constrained experimentally (Schulzlessdorf and Hedrich, 1995;Peiter et al., 2005;dacz-Narloch et al., 2011), but this channel is unsuited to the task of Ca 2+ release using any of the range of gating and permeation parameters that are reported in the literature (see Hills et al., 2012). Without self limitation, trials of the SV channel in the OnGuard model either failed to release Ca 2+ and to drive stomatal closure, or yielded prolonged elevations, lasting tens of minutes to hours, with [Ca 2+ ] i pinned to values in excess of 50 to 100 mM until the vacuolar Ca 2+ content dropped below 3% of its mean daytime value. Such behavior is not consistent with the physiology of the guard cell.
To overcome this difficulty we included in the model (Hills et al., 2012) a Ca 2+ channel (TonVCa) with [Ca 2+ ] i activation analogous to the IP 3 -receptor channels of the Beta vulgaris vacuole (Alexandre et al., 1990). We added to this channel self inhibition with time-dependent inactivation and reactivation above and below fixed [Ca 2+ ] i thresholds, respectively, similar to the slow inactivation and latency of the animal IP 3 -receptor channels (Bezprozvanny et al., 1991). Finally, gating of this hypothetical Ca 2+ channel was assigned a [Ca 2+ ] v sensitivity to its voltage dependence to avoid depleting the endomembrane Ca 2+ store. These characteristics remain undocumented, and therefore should be viewed as placeholders predicting the characteristics of one or more endomembrane Ca 2+ release pathways. They may represent features of a tonoplast channel, but need not be restricted to this membrane; in vivo intracellular Ca 2+ release depends on other compartments including the endoplasmic reticulum (Clapham, 1995;Bootman et al., 1997;Hamilton et al., 2000;Navazio et al., 2001;Garcia-Mata et al., 2003;Sokolovski et al., 2008). In short, the TonVCa characteristics are a predictive summation of the dynamic outputs essential for the [Ca 2+ ] i physiology of the guard cell, rather than the literal descriptors for a Ca 2+ channel yet to be identified.
These caveats aside, the TonVCa channel successfully constrained [Ca 2+ ] i excursions to maxima near 1 to 2 mM, and it resulted in simulated [Ca 2+ ] i oscillations with periods of 5 to 10 s to 10 to 20 min, characteristics qualitatively and quantitatively similar to experimental measurements in vivo (Gradmann et al., Figure 7. [Ca 2+ ] i and analysis of the energetics for Ca 2+ flux across the plasma membrane and tonoplast resolved over the diurnal Reference Cycle as described in Figure 1 (12 h light:12 h dark, indicated by bars above). Shown in A and F are the current-voltage curves for the dominant Ca 2+ transporters at the plasma membrane and tonoplast, respectively. Total membrane current is indicated in each case by the dotted line; the free-running membrane voltage is defined by the point at which this line crosses the voltage axis. Also summarized are the plasma membrane and tonoplast voltages (B), the [Ca 2+ ] i (C), and the Ca 2+ channel and Ca 2+ -ATPase fluxes at the plasma membrane (D) and tonoplast (E). Current-voltage curves are cross-referenced to time points in B to E by letter. Time point b is characterized by plasma membrane hyperpolarization, which favors Ca 2+ influx through the Ca 2+ channels at the plasma membrane (A); time point c is characterized by plasma membrane depolarization, which favors Ca 2+ efflux through the plasma membrane Ca 2+ -ATPase. Note the difference in scales between D and E, which underlines the predominance of Ca 2+ circulation between the cytosol and the endomembrane store of the vacuole. Note that positive flux is defined as movement of the ionic species (not the charge) out of the cytosol, either across the plasma membrane or the tonoplast. 1993; Blatt, 2000b;Blatt et al., 2007;McAinsh and Pittman, 2009). Figure 8 shows the logged output from the train of [Ca 2+ ] i and voltage excursions simulated by the OnGuard model at the end of the daylight period (see also Fig. 7) and associated IV curves taken from screenshots at the time points indicated. Also plotted are the associated fluxes of Ca 2+ through the Ca 2+ channels and Ca 2+ -ATPases at the two membranes and the change in stomatal aperture. Analysis of the flux and voltage over time indicated within each cycle a clear sequence of events that may be summarized as follows.
(1) Initially, the negative membrane voltage drives Ca 2+ influx through Ca 2+ channels at the plasma membrane, raising [Ca 2+ ] i (time point a).
(2) This rise in [Ca 2+ ] i promotes TonVCa channel activity and Ca 2+ influx across the tonoplast. (3) A threshold is reached near the K Ca for TonVCa channel activation, beyond which the positive feedback of [Ca 2+ ] i leads to self activation of the TonVCa channel (timepoint b). (4) Elevated [Ca 2+ ] i stimulates SLAC and R-(ALMT-)type channels and suppresses the H + -ATPase at the plasma membrane, consistent with the K Ca for each transporter (see Hills et al., 2012); initially the effect is to generate an N shape and two quasi-stable voltages (two inflections of positive slope with the voltage axis) in the total plasma membrane IV curve before the membrane collapses to the more positive of these two voltages (timepoint c). (5) Elevated [Ca 2+ ] i enhances the activities of the Ca 2+ -ATPases at both membranes and, with depolarization, the kinetic restriction on the plasma membrane Ca 2+ -ATPase is relieved. (6) With [Ca 2+ ] i above its inactivation threshold, the TonVCa channel passes slowly into inactive latency, thereby shifting the flux balance toward Ca 2+ export and a gradual recovery in [Ca 2+ ] i (time point d). (7) As [Ca 2+ ] i falls, the activities of the SLAC and R-(ALMT-)type anion channels decline and the H + -ATPase recovers until the plasma membrane hyperpolarizes again. (8) Finally, the TonVCa channel escapes latency as [Ca 2+ ] i falls below its minimum threshold.
The parallels to neuromuscular action potentials are striking (Jack et al., 1983;Hille, 2001), notably the underlying feedback between membrane voltage, Ca 2+ entry and Ca 2+ -induced Ca 2+ release, [Ca 2+ ] i , and membrane voltage recovery. In this case, however, each depolarizing cycle at the plasma membrane and tonoplast promoted K + and anion (Cl 2 and Mal 22 ) flux from the vacuole, through the cytosol, and out across the plasma membrane, resulting in a loss in turgor pressure and volume. As a consequence, the OnGuard model yielded an oscillatory decline in stomatal aperture, each step-like fall in aperture arising from the loss in osmotic solutes and water during the periods of depolarization (Fig. 8A). This prediction accords with experimental observations of stomatal closure following a Ca 2+ stimulus (Wang et al., 2006), and it may also relate to oscillations in transpiration long known to occur in leaves of several plant species (Lang et al., 1969;Cowan, 1972;Farquhar and Cowan, 1974;Cardon et al., 1994).
Other features of the oscillations generated by OnGuard model are equally noteworthy. We highlight three of these here. First, the model predicts that the initial rise in [Ca 2+ ] i is linked to hyperpolarization of the plasma membrane, much as was postulated a decade ago from experimental data for the Ca 2+ channels and [Ca 2+ ] i recordings (Grabov and Blatt, 1998;Hamilton et al., 2000Hamilton et al., , 2001Blatt, 2000b). These earlier studies showed that Ca 2+ channel activity is strongly dependent on membrane voltage and promoted by membrane hyperpolarization, thus favoring [Ca 2+ ] i elevation; they also indicated slower, Ca 2+ -sequestering and export activities that must underpin the recovery phase of the [Ca 2+ ] i and voltage. Two conclusions to be drawn from the experimental data, and from the simulations above, are that the rise in [Ca 2+ ] i is strongly dependent on membrane voltage and that, once elevated, membrane voltage is dependent on [Ca 2+ ] i recovery. In other words, [Ca 2+ ] i and plasma membrane voltage form two major arms of a cyclic control loop that regulates plasma membrane transport. The simulations also predict that the recovery phase, but not the initial transient depolarization and [Ca 2+ ] i elevation phase, is dependent on Ca 2+ export at the plasma membrane and especially at the tonoplast. Thus, suppressing or eliminating these activities-either of the Ca 2+ -ATPases or of CAX Ca 2+ -H + antiport activity via block of vacuolar H + transport-is predicted to slow the recovery. Some experimental evidence supports this prediction: Notably, mutations affecting the vacuolar H + -ATPase and tonoplast energization greatly prolong [Ca 2+ ] i elevations (Allen et al., 2000). Finally, a third conclusion to arise from the simulations is the fundamental importance of H + -ATPase inhibition at the plasma membrane by elevated [Ca 2+ ] i . Trials with the H + -ATPase uncoupled from [Ca 2+ ] i -that is, either independent of [Ca 2+ ] i or lacking an appreciable inhibition with [Ca 2+ ] i elevations above 300 to 400 nMfavored prolonged nonphysiological increases in steady-state [Ca 2+ ] i . Analysis of the corresponding IV curves showed the behavior to arise from insufficient inward current to drive the voltage positive and bias the membrane for Ca 2+ efflux. In general, the H + -ATPase held the membrane voltage sufficiently negative that the Ca 2+ -ATPase was kinetically restricted, unable to keep pace with the Ca 2+ influx passing through the Ca 2+ channels at the plasma membrane. Thus the model predicted a progressive accumulation of Ca 2+ in the vacuole and a pronounced elevation of [Ca 2+ ] i .

DISCUSSION
Stomatal dynamics have long been incorporated into models, notably to predict gas-exchange characteristics at the level of the plant community (Ball, 1987;Williams et al., 1996;Eamus and Shanahan, 2002). These models subsume the cellular mechanics of the guard cell within a few empirical parameters of linear hydraulic pathways and conductances. However, no model of stomata to date has taken full advantage of the wealth of knowledge available at the cellular level for guard cell transport, signaling, and homeostasis. The complexity of the guard cell transport in itself precludes any quantitative description of how guard cells regulate stomatal aperture, let alone a clear understanding of the emergent properties of the guard cell system as a whole. To address this gap in understanding, we have taken a computational approach to dynamic modeling of the guard cell. The OnGuard software and model described in the preceding article (Hills et al., 2012) and elaborated above incorporate all of the fundamental properties for transporters at the plasma membrane and tonoplast, and the salient features of osmolite metabolism. Here we demonstrate its capacity for true predictive power in generating a number of counterintuitive outputs in the emergent behavior of the system, using the diurnal cycle as a testbed. Many of these outputs find support in experimental data already extant in the literature, and we provide additional validation in support of our starting parameters for primary transport activity from an analysis of gene and protein expression for the tonoplast H + -ATPase. Of particular note, we show that the OnGuard model faithfully reproduces the diurnal variations in K + , Cl 2 , and Mal flux and content, in cytosolic and vacuolar pH, and in [Ca 2+ ] i that have been reported in the literature. Furthermore, it , the Ca 2+ influx through the Ca 2+ channels at the plasma membrane and tonoplast (C), the Ca 2+ efflux through the plasma membrane and tonplast Ca 2+ ATPases (D), and the [Ca 2+ ] i (E). Current-voltage curves for these Ca 2+ transporters at the plasma membrane and tonoplast are shown in E and F, respectively. Total membrane current is indicated in each case by the dotted line; the free-running membrane voltage is defined by the point at which this line crosses the voltage axis. Current-voltage curves are cross-referenced to time points in B to E by letter. As in Figure 7, time point a is characterized by plasma membrane hyperpolarization, which favors Ca 2+ influx through the Ca 2+ channels at the plasma membrane; time point c is characterized by plasma membrane depolarization, which favors Ca 2+ efflux through the plasma membrane Ca 2+ -ATPase. Note the Ca 2+ influx at the plasma membrane is replaced by Ca 2+ influx across the tonoplast shortly before plasma membrane depolarization (C); thereafter, the characteristics of the [Ca 2+ ] i rise are determined by Ca 2+ flux from the tonoplast and, as the tonoplast Ca 2+ channels inactivate, by [Ca 2+ ] i recovery-mediated Ca 2+ -ATPases at the two membranes (D). Note that positive flux is defined as movement of the ionic species (not the charge) out of the cytosol, either across the plasma membrane or the tonoplast.
generates action-potential-like variations in membrane voltage and in [Ca 2+ ] i similar to spontaneous oscillations observed experimentally and postulated to play a role in control of stomatal aperture, and it predicts these oscillations to facilitate changes in the osmotic solute content of the guard cell. Significantly, the OnGuard model arrives at each of these predictions without ad hoc assumptions of signaling pathways to connect the various transport and metabolic activities. Other model outputs will now fuel substantive research projects in their own right and are beyond the scope of this study. In short, the OnGuard software establishes a framework for the systems biological analysis of stomatal guard cells, and it sets out a flexible modeling environment that should find application in explorations of similar physiological and related problems in the future.

Guard Cell Systems Modeling
Of the few models developed, bottom up from knowledge of guard cells, the only recent effort (Li et al., 2006) surprisingly made deliberate use of Boolean network analysis, thus ignoring the large body of detailed kinetic information available. This choice was not an obvious one. The difficulty with a Boolean approach is its inability to encompass the kinetic dimension essential to dynamic modeling of physiological processes. The power of such network analysis is as a tool to study systems for which there is little quantitative information and a large number of components and possible connections between them. Its most common application is to identify critical components-for example, to answer the question whether signal transmission might be blocked by eliminating one component or its connection to another-to map the major nodes and causal pathways of a network. However, the simplicity of Boolean models, in which components and their links can only be on or off, precludes any useful application in studies that aim to explore emergent behaviors arising from dynamic interactions within a well-defined network. Of course, it is possible to apply pseudotemporal characteristics to Boolean models, but the outputs are inevitably disconnected from any meaningful comparison with physiology. That the Boolean guard cell model (Li et al., 2006) yielded an output that had been demonstrated experimentally more than a decade earlier (Blatt and Armstrong, 1993), but went without mention, does little to support the claim to its predictive power.
Some past success was achieved in recapitulating guard cell voltage oscillations based on a computational approach and provided an impetus for the work we present here. Gradmann et al. (1993) incorporated the predominant K + and Cl 2 channels, H + -ATPase, and H + -Cl 2 symport activities to demonstrate the minimum of transporters necessary to simulate the oscillations in voltage observed to arise spontaneously in guard cells. This modeling effort was limited in several respects. Among others, it did not include components necessary for Ca 2+ transport, or parameters for the sensitivities of the various transporters to [Ca 2+ ] i and pH. The latter omissions are especially important, because pH and [Ca 2+ ] i are known to exert control on a number of these transport processes. Furthermore, it failed to satisfy the basic requirements for a stable Reference State or the equivalent-lacking a sufficient complement of pathways to balance the modeled ion fluxes, a minimum that might be dictated by our Reference State analysis-and thus incorporated substantial predictive indetermination. Nonetheless, the Gradmann model was successful in demonstrating voltage oscillations similar to those observed in vivo, and it led to the proposal of a time-averaging mechanism controlling K + and anion fluxes across the plasma membrane, elements of which are evident in the OnGuard model (see Figs. 7 and 8). Analogous results were achieved in an application to ion transport in complex tissues of the plant root (Shabala et al., 2006), although similar limitations and lack of indetermination undermine its predictive utility.

Predictive Power of the OnGuard Model
The few instances in which mathematical modeling has been applied to cellular homeostasis with sufficient rigor have been remarkably successful both in reproducing known cellular physiology and in predicting unexpected behaviors. For example, dynamic models of cellular homeostasis (Lew et al., 1979) correctly predicted a transient cell shrinkage and protracted fall in epithelial short-circuit current following ouabain inhibition of the Na + /K + -ATPase. Both predictions were counterintuitive but experimentally confirmed (Macknight et al., 1975a(Macknight et al., , 1975b) and, at the time, appeared to contradict the general validity of Ussing's now widely accepted, two-barrier description of mammalian epithelia (Ussing, 1982). The OnGuard model demonstrates a similar, if not a greater, degree of success in integrating knowledge of stomatal physiology and generating novel insights. Validation of the model using a set of basic parameters within a simple, diurnal regime gave a remarkably close approximation to the guard cell behavior in vivo and its variation with the diurnal cycle; physiologically sensible outputs were obtained for guard cell turgor pressure, volume, partial volumes of the cytosol and vacuole, stomatal aperture, and values for total osmotic solute content as well as compartmental ion and metabolite concentrations, cytosolic and vacuolar pH, plasma membrane, and tonoplast voltage (Figs. 1-6).
In addition, the OnGuard model generated a surprising wealth of predictions, many of which we have highlighted above. Among these, the model yielded an apparent exchange of Cl 2 with Mal, and counterintuitive changes in [Ca 2+ ] i and pH over the diurnal cycle, each of which finds direct support in independent experimental data. It is especially significant that these outputs were achieved without additional feedback controls, such as those suggested to affect the diurnal capacity for endomembrane Ca 2+ release (Dodd et al., 2007). Thus, these simulations demonstrate the capacity of the OnGuard model to encapsulate otherwise unexpected biological phenomenology within the emergent properties of the known network of intrinsic kinetic parameters of the guard cell. Other predictions that find support in experimental literature (Blatt, 2000b;Schroeder et al., 2001;Blatt et al., 2007;McAinsh and Pittman, 2009) arose from simulated trains of short-term oscillations in [Ca 2+ ] i that took place over periods of seconds to minutes. These demonstrated (1) a coupling to plasma membrane voltage during the early phase of [Ca 2+ ] i elevation and, conversely, the dependence of membrane voltage on [Ca 2+ ] i recovery, (2) the central role for anion channels in driving the transition between [Ca 2+ ] i elevation and its recovery, and (3) the fundamental importance of H + -ATPase inhibition at the plasma membrane by elevated [Ca 2+ ] i to engage the recovery phase of the oscillatory cycle (see . Significantly, all of these predictions arise from the basic Reference Cycle that we implemented as a starting point for exploratory simulations.
We have yet to extend our studies to the effects, for example, of environmental and hormonal challenges or of single transporter mutations and genetic ablations, and therefore anticipate many more fundamental insights to come of this work in the future.

Limitations of the OnGuard Model
These successes notwithstanding, the OnGuard model is subject to a number of restrictions that will need revisiting as our modeling efforts advance. Of these, the inclusion of additional intracellular compartments with different dynamic Ca 2+ and H + buffering (sequestration) capacities and K d values is likely to prove important in fine tuning signaling events and, for example, could allow for the midday depression in stomatal aperture previously documented in several species (Dodge et al., 1992;Hirasawa and Hsiao, 1999). We have refrained from adding compartments other than the vacuole at present, because there remains insufficient experimental data to adequately constrain the parameters without predictive indetermination. Similarly we have yet to introduce mechanisms that accommodate protein (de-) phosphorylation cascades or changes in redox status (Wang and Song, 2008;Kim et al., 2010). These, and similar omissions do not mean that such processes are unimportant, but rather that specific kinetic detail has yet to become available that would justify their explicit inclusion. Subsuming their effects within the functioning of the individual target processes offers an alternative, effectively black boxing the underlying mechanism(s) within the phenomenological behavior of the targets (Endy and Brent, 2001). This approach is useful, too, in reducing complexity and computational burden without a loss in predictive power. There are many examples of the successful use of black-box phenomenology, including the Hodgkin-Huxley equations used to represent the operation of the Na + and K + channels responsible for neuronal action potentials; the Hodgkin-Huxley equations explained the fundamental physiological processes of channel gating long before the underlying molecular mechanisms were elucidated (Hille, 2001). We have devised mathematical black boxes that effectively subsumed ensembles of passive conductances across the plasma membrane of Neurospora (Blatt and Slayman, 1987), Chara (Blatt et al., 1990a), and Vicia guard cells (Blatt, 1987a), and the same approach has been used in the past (Ferreira and Lew, 1976;Raftos et al., 1999), and in this OnGuard model (Hills et al., 2012), to accommodate cytosolic proton and divalent cation buffering reactions.
In conclusion, we have built on the concept of a Reference Cycle, using the diurnal cycle of activities known for the primary, ATP-dependent pumps, Suc, and Mal metabolism to generate a number of quantitative, and often counterintuitive, predictions emergent in the behavior of the stomatal guard cell. Significantly, the OnGuard model yielded each prediction from a single, common framework of model parameters and without ad hoc assumptions of signaling pathways to connect the various transport and metabolic activities. Thus, we are able to demonstrate the capacity of the OnGuard software to generate a model with true predictive power for the stomatal guard cell.

OnGuard Modeling
The OnGuard software (www.psrg.org.uk) was driven through a diurnal 12:12 h light:dark cycle as described in the text. Light intensity was elevated in a linear ramp from zero at the start of the cycle to 2,200 mmol m 22 s 21 at 6 h into the cycle, and then was reduced in the same manner to zero at 12 h into the cycle. All light-sensitive processes-the transport ATPases and PPase, and Suc synthesis-were assigned fluence dependencies with a K 1/2 of 50 mmol m 22 s 21 with minimum (dark) and maximum (light) activities as described previously (Hills et al., 2012). All other model parameters were fixed. Parameter values will be found in supplemental appendix S6 of the companion article (Hills et al., 2012) and are available with the OnGuard software (www.psrg.org.uk).

Supplemental Data
The following materials are available in the online version of this article.
Supplemental Figure S1. Diurnal variation in VH + -ATPase transcript levels from Arabidopsis leaf tissue.
Supplemental Figure S2. Diurnal variation in VH + -ATPase activity from Arabidopsis leaf tissue.
Supplemental Materials and Methods S1.