- © 2012 American Society of Plant Biologists. All rights reserved.

## Abstract

Stomatal guard cells play a key role in gas exchange for photosynthesis while minimizing transpirational water loss from plants by opening and closing the stomatal pore. Foliar gas exchange has long been incorporated into mathematical models, several of which are robust enough to recapitulate transpirational characteristics at the whole-plant and community levels. Few models of stomata have been developed from the bottom up, however, and none are sufficiently generalized to be widely applicable in predicting stomatal behavior at a cellular level. We describe here the construction of computational models for the guard cell, building on the wealth of biophysical and kinetic knowledge available for guard cell transport, signaling, and homeostasis. The OnGuard software was constructed with the HoTSig library to incorporate explicitly all of the fundamental properties for transporters at the plasma membrane and tonoplast, the salient features of osmolite metabolism, and the major controls of cytosolic-free Ca^{2+} concentration and pH. The library engenders a structured approach to tier and interrelate computational elements, and the OnGuard software allows ready access to parameters and equations ‘on the fly’ while enabling the network of components within each model to interact computationally. We show that an OnGuard model readily achieves stability in a set of physiologically sensible baseline or Reference States; we also show the robustness of these Reference States in adjusting to changes in environmental parameters and the activities of major groups of transporters both at the tonoplast and plasma membrane. The following article addresses the predictive power of the OnGuard model to generate unexpected and counterintuitive outputs.

Stomatal guard cells surround pores in the epidermis of plant leaves and regulate the pore aperture. They open the pore in response to low CO_{2} and light to facilitate CO_{2} access for photosynthesis, and they close the pore in the dark, under drought stress, and in the presence of the water-stress hormone abscisic acid to minimize water loss through transpiration. Stomata have a profound impact on the water and carbon cycles of the world (Gedney et al., 2006; Betts et al., 2007). Their dynamics have been incorporated into models for transpiration and water use efficiency (Farquhar and Wong, 1984; Ball, 1987; Williams et al., 1996; Eamus and Shanahan, 2002; West et al., 2005), successfully reproducing the gas exchange, CO_{2}, and transpirational characteristics of experiments at the plant and community levels. To date, these models have taken a top-down approach. They subsume stomatal movements within a few empirical parameters of linear hydraulic pathways and conductances without reference to the molecular mechanics of the guard cell. No generalized guard cell model has yet to be developed from the bottom up, drawing on the wealth of knowledge available for guard cell transport, signaling, and homeostasis. It is clear that such a model is now needed. The depth and breadth of information available for stomatal guard cells has made them the premier cell system in plants for studies of membrane transport, signaling, and homeostasis; yet, in face of the complexity of the guard cell system, there remains a very large gap in quantitative understanding about how guard cell transport works together to modulate solute flux and regulate stomatal aperture.

A very large body of experimental evidence supports the collective role of ionic fluxes across both plasma membrane and tonoplast, and the metabolism of Suc and malate (Mal) in shaping the changes in osmotic load and turgor pressure that drive stomatal opening and closing (Willmer and Fricker, 1996; Blatt, 2000; Schroeder et al., 2001; Hetherington and Woodward, 2003; Sokolovski and Blatt, 2007). All of the predominant transporters—the major K^{+}, Cl^{−}, and Mal-permeable channels and the H^{+}-ATPases at the plasma membrane and tonoplast, as well as Ca^{2+}-permeable channels at both membranes (Allen and Sanders, 1997; MacRobbie, 1997; Sanders et al., 2002; White and Broadley, 2003; Dreyer et al., 2004; Peiter et al., 2005; Sokolovski and Blatt, 2007)—have each been isolated and studied in sufficient depth to provide detailed and accurate flux equations with parameters fully constrained by experimental data. Much detail is available relevant to the activities of soluble enzymes, the metabolic pathways, and their kinetics that determine the synthesis, interconversion, and breakdown of Suc and Mal within guard cells (Willmer and Fricker, 1996).

The information this body of data represents and their nature poses a major challenge for its incorporation and integration within any quantitative kinetic model. Not only is there, for each transport process, a unique set of kinetic and regulatory descriptors, but for the majority of transporters the process itself is inherently recursive, acting on one or more of these descriptors. For example, gating of the outward-rectifying K^{+} channels is sensitive to membrane voltage as well as K^{+} concentration; depolarizing the membrane promotes the current (Blatt, 1988b, Blatt and Gradmann, 1997; Blatt and Gradmann, 1997), yet the current generated on activating these channels normally draws K^{+} out of the cytosol, thereby driving the membrane voltage negative and suppressing channel gating. Indeed, each transport process carrying charge across a membrane affects—and is affected by—the voltage across that membrane, if only as a consequence of mass action and the movement of charged ions it carries. The problem is all the more complex because, for charge transporters operating across a common membrane, voltage is both substrate and product and is shared between all of these charge-carrying transport processes (Weiss, 1996; Blatt, 2004a). Thus, the challenge becomes one of systematically integrating each and every one of these processes in a way that accommodates the recursive nature of transport and within an overarching strategy that is sufficiently flexible to allow parameter modifications, even substitutions for the equations representing each process, on the fly during experimentally guided model refinements.

To this end, we have expanded on an integrated approach pioneered in models of mammalian epithelia and erythrocytes (Lew et al., 1979, 2003; Lew and Bookchin, 1986; Mauritz et al., 2009), incorporating in an iterative computational strategy the ensemble of transport and buffering equations and their associated variables assigned by the operator. We describe the concepts behind our development of the programming library for Homeostasis, Transport and Signalling (HoTSig) and construction of the OnGuard software for dynamic modeling of the guard cell. An important feature of the HoTSig library is its open structure, standardized with the major sets of equations and descriptors for each of the various types of transporters. The OnGuard software incorporates a graphical user interface (GUI) that integrates with the Microsoft Windows platforms and a Reference State Wizard for defining the starting point for in silico experimentation. Here we show that OnGuard-generated models readily achieve stability in a set of physiologically relevant Reference States associated with the open and closed states of stomata, and we explore the subsets of transporters and their parameters to which the models are especially sensitive. We show that these models are robust in adjusting to changes in environmental parameters. The predictive power of the OnGuard modeling approach is demonstrated with selected kinetic simulations in the companion article (Chen et al., 2012).

## MODELING

### The Modeling Strategy

Cellular homeostasis is especially well suited to an integrative, bottom-up mathematical modeling approach. The physicochemical relations that constrain the behavior of homeostatic variables—conservation of mass, charge, the electroneutrality of intracellular and extracellular solutions, the osmotic transport of water linked to all solute concentrations and their gradients, and of membrane potential linked to all ion gradients and permeabilities—are simple quantitative relations, all easy to model. For the guard cell, the most relevant homeostatic variables comprise cell volume, cell osmolality, water potential and turgor, membrane potential, cytosolic and vacuolar K^{+}, Cl^{−} and to lesser extents total and free Ca^{2+} contents and concentrations, pH, unidirectional and net fluxes of all ionic or neutral solutes through each transporter, water permeability, and flux. Also important from the standpoint of transport and metabolic regulation are the intracellular proton and Ca^{2+} buffering systems, the cell content of impermeant solutes, their osmotic coefficients, charge, and its dependence on pH. For most solutes, the lipid bilayer provides an effectively impermeable barrier. Thus transport aside, the relation between the free and total concentrations of each solute depends on factors such as buffering, macromolecular binding, and the degree of ionization. For all of the major solutes, quantitative data including buffering constants are available or can be estimated for endogenous buffer systems as well as for all of the experimentally applied buffers and chelators (e.g. HEPES, EGTA) in common use (Ferreira and Lew, 1976; Tsien and Tsien, 1990; Föhr et al., 1993; Grabov and Blatt, 1997; Tiffert and Lew, 1997).

The biophysical relations of the different types of transport that occur across eukaryotic membranes—mediated by pumps, channels, and carriers—are all well understood and for the predominant transporters have been studied in sufficient depth to provide detailed and accurate flux equations (Weiss, 1996). The predominant ATP-driven ion pumps, H^{+}-coupled transporters, and passive ion channels all have been characterized with regard to stoichiometry and mechanism, either in guard cells or in other plant cell types (Sanders, 1990; Willmer and Fricker, 1996; Blatt, 2004b), and their operation is readily described by sets of kinetic equations. Much of this kinetic detail applies directly to the transporters in situ at the plasma membrane. Our knowledge of transport across the tonoplast is less well developed, largely because of difficulties to gain access in vivo. In principle, the relative ignorance of transport at the tonoplast leaves quantitative modeling open to indetermination with unknown parameters. Nonetheless, a large body of solid experimental results can be translated as constraining equations on the system, including data on vacuolar ion contents and fluxes (MacRobbie, 1995, 2000, 2002; Willmer and Fricker, 1996; Gobert et al., 2007) as well as the biophysical constraint of osmotic equilibrium between cytosol and vacuole. These constraints link the elemental transport and chemical processes to the macroscopic variables of the system such as the total guard cell volume (*VT*), turgor pressure (*P*_{T}), and stomatal aperture (*A*_{S}). (A complete list of abbreviations is provided in Supplemental Appendix S5.) Thus, by careful design of the boundary equations it is possible to minimize the range of values for the set of unknown parameters that complies with the macroscopic behavior of the system. This approach guided the design of the modeling reported here.

Key information is available also for the regulation of transport as well as that of Mal and Suc metabolism (Willmer and Fricker, 1996). Again, gaps exist in our understanding of these controls but, from the context of a modeler seeking to understand how the system responds to perturbation, the only relevant biology is encapsulated by how one model variable is connected to another. It is frequently the case that this phenomenology is directly accessible to experiment, whereas the underlying mechanistic details are not, or are accessible only qualitatively. For example, it is well known that an elevated [Ca^{2+}]_{i} inactivates I_{K}_{,in} of the guard cell, but we can only surmise that this may occur by its activation of a CBL-dependent protein kinase (Xu et al., 2006; Cheong et al., 2007) and we do not have sufficient quantitative information to model kinase activation or K^{+} channel phosphorylation. Nevertheless, the quantitative relation between [Ca^{2+}]_{i} and K^{+} channel activity is known (Grabov and Blatt, 1999) and, by applying a mathematical description of this relation, we can safely place the mechanistic details in a black box that subsumes the intermediate kinetics. There are many other examples of the successful use of black-box phenomenology, including the Hodgkin-Huxley equations, which described the fundamental physiological processes of the Na^{+} and K^{+} channels responsible for action potentials in axons long before the underlying molecular mechanisms were elucidated (Hille, 2001). Such black boxes effectively parameterize phenomenological modules (Endy and Brent, 2001) and may be opened if, and when any elements subsumed within a module become a target of the modeling. This phenomenological approach also reduces complexity and computational burden.

### The HoTSig Library and OnGuard Software

By contrast with analytical approaches, numerical computation allows flexibility in modifying these representations during experimentally guided model refinements. It is therefore important to resist temptations to seek solutions that may appear to provide simplified or explicit analytical equations for voltage and the sum of membrane ion fluxes, but that sterilize the main objectives of the modeling exercise. We expanded on an iterative computational strategy applied successfully in the past (Lew et al., 1979, 2003; Lew and Bookchin, 1986; Mauritz et al., 2009), developing a library and software that incorporates the ensemble of transport and buffering equations and their associated variables assigned by the operator. These equations were included together with a set of macroscopic descriptors of guard-cell-specific features in construction of the OnGuard software. In operation, the OnGuard software calculated and logged the dynamic adjustments of ion flux and compartmental composition and of membrane voltage from a set of operator-defined starting conditions; it used the sets of nonlinear differential flux equations for the transporters; and it obeyed the fundamental physical constraints of mass and charge conservation. We built into the software the facility to adjust the integration interval, taking account of the sums of each of the ionic fluxes such that the duration of iteration steps and the frequency of data output could be set to adjust with the rate of change in the system (Lew and Bookchin, 1986; Lew et al., 1991), although in practice even short integration periods of 1 ms are handled with little loss of speed when run on a quad-processor-based computer typical of those now widely available commercially. We provide a brief explanation here. Further details of the HoTSig library and the OnGuard software are provided in Supplemental Appendices S1 to S4.

#### Compartmental Analysis

The utility of a model relies on its facility for comparison of its output with experimentally observed behaviors. The vast majority of studies at the cellular level continue to be carried out on guard cells isolated from the leaf in epidermals peels, or as protoplasts and isolated vacuoles, and maintained in controlled external media. As a starting point, we therefore treated the apoplast surrounding guard cells as an infinite reservoir, uninfluenced by material entering or emanating from the guard cell, leaving its composition to be defined by the operator as required for each particular simulation. We also assumed a single endomembrane compartment, hereafter referred to as the vacuole. This simplification avoided the need to define additional sets of poorly constrained transporters for other endomembrane compartments, although these may be added later once sufficient experimental detail becomes available. Extracellular medium, guard cell cytosol, and vacuole therefore defined a semiclosed system of compartments in series.

#### Preservation of Electroneutrality

The capacitative currents that charge the membrane are both transient and many orders of magnitude smaller that those mediated by the ion fluxes relevant to cellular homeostasis (Findlay and Hope, 1976; Jack et al., 1983), and can thus be neglected for purposes of modeling here. Electroneutrality is preserved by the instant value of the membrane potential across each of the membranes in the system, *V*_{pm} and *V*_{ton}. These values ensure that the sums of the individual ionic currents through each of the electrogenic transporters on the tonoplast and plasma membrane are zero at all times. Thus, by setting the sum of all currents *I*_{i} at time *t* to zero, that is Σ*I*_{i(t)} = 0 at each membrane, and solving these implicit equations for *V*_{pm} and *V*_{ton}, respectively, we derived their values in each iteration, thereby ensuring that electroneutrality was preserved throughout all simulations. For each membrane, Σ*I*_{i(t)} is necessarily a complex function containing all of the different kinetic representations of the various pump, channel, and carrier-mediated transporters.

#### Formulating the Constraining Equations

Stomatal dynamics are determined by the regulated gain and loss of osmotically active solutes, *Q*_{i}, within the guard cells and the associated changes in osmotic and turgor pressure (Raschke, 1979; MacRobbie and Lettau, 1980; Willmer and Fricker, 1996). The computational end product of the ensemble of all individual transport and metabolic processes is the instantaneous value of the sum of all osmotically active solutes within the guard cell [Σ*Q*_{i}(*T*)], contained both within vacuole [Σ*Q*_{i}(*v*)] and cytoplasm [Σ*Q*_{i}(*c*)]:Two main constraints apply here, based on reliable data (Hill and Findlay, 1981; Willmer and Fricker, 1996): (1) the tonoplast cannot sustain hydrostatic pressure differences, and (2) the water permeability of both tonoplast and plasma membrane are high. Therefore, the osmotic pressure difference between vacuole and cytoplasm can be assumed to be zero and the turgor pressure of the guard cell related solely to the osmotic potential difference across the plasma membrane. If, for simplicity, we rename *QT* = Σ*Q*_{i}(*T*), *Qv* = Σ*Q*_{i}(*v*), and *Qc* = Σ*Q*_{i}(c), and define *VT*, *Vv*, and *Vc* as the total volume of the guard cell, the vacuole, and the cytosol compartments, respectively, then the new labels and the two constraints define the equalities:The osmotic pressure (Π) at equilibrium with turgor pressure *P*_{T} at each instant of time is given by the Van’t Hoff equation:where ΣC_{apo} is the sum of the concentrations of all osmotically active solutes in the apoplast. Solving for *VT* yields:Here, *VT* is linked with *P*_{T} for each value of *QT. QT* in turn represents the end result of a large computational sequence in each iteration of the model and reflects the osmotic load of the guard cell at that time. The two unknowns in Equation 8, *VT* and *P*_{T}, will be complex functions of cell wall plasticity and its elastic modulus (Cosgrove, 1987). However experimental data links these variables empirically through the properties of the guard cell wall and the geometry of the stomata (Raschke, 1979; MacRobbie and Lettau, 1980; Willmer and Fricker, 1996; Franks et al., 2001) and satisfy the requirement for a constraining equation that describes the link between *VT* and *P*_{T}.

The relations between *P*_{T} and stomatal aperture, *A*_{S}, and between *VT* and *A*_{S} have been measured for a number of species, including *Vicia*, *Commelina*, tobacco (*Nicotiana tabacum*), and Arabidopsis (*Arabidopsis thaliana*; MacRobbie and Lettau, 1980; Franks et al., 2001; Shope et al., 2003; Gay, 2004; Meckel et al., 2007; Sokolovski et al., 2008). Both relations can be approximated by linear expressions with varying slopes and baseline stomata closure values such that:andwhere *m*, *n*, *r*, and *s* are empirically determined constants.

From these equations, we can derive the desired relation between *VT* and *P*_{T} using *A*_{S} as a parametric variable. Substituting from Equation 9 thus gives:and from Equation 10:Equating Equations 11 and 12 gives:where *p = m*/*r* and *q = n* − *m·s*/*r*. Replacing *P*_{T} from Equation 13 in Equation 9 builds in the constraints encoded by the experimentally determined Equations 9 and 10, and enables a solution for *VT* and *QT* in each cycle of computation. Replacing *P*_{T} from Equation 8 in Equation 9 yields:and rearranging gives:a quadratic equation in *VT* with the single solution:With the values of *QT* and *VT*, the values of *Vc*, *Vv*, *P*_{T}, and *A*_{S} at any given time can be derived as:In this strategy, all the particular geometrical features and wall properties of the guard cells are phenomenologically encoded within the linear coefficients in Equations 9 to 13, allowing the operator to test species-related behavior based on the availability of experimental data that define these coefficients. This strategy delivers the value of all the macroscopic variables of the system (*VT*, *P*_{T}, *A*_{S}) at each instant of time, allowing comparisons between predicted and experimentally observed results. The critical micromacro link in this strategy is the value of the osmotic load at each instant of time, *QT*(t), and its computation is the focus of the sections below.

#### Computing the Osmotic Load

The osmotic load, *QT*, of the guard cell responsible for stomatal dynamics arises from four main processes comprising (1) membrane transport, (2) buffering reactions, (3) metabolic synthesis and degradation reactions, and (4) signaling and regulatory reactions. There are many subgroups within each of these main categories. Supplemental Tables S1 to S6 list the macroscopic constraints and the associated processes encoded in this initial OnGuard formulation. Details of the software implementation and additional explanations behind the choices in formulation will be found in Supplemental Appendices S1 to S4 and Figs. S1-S3. Supplemental Tables S1 and S2 summarize the cellular and compartmental characteristics typical of guard cells, with specific reference to those of *Vicia*. Supplemental Tables S3 to S6 outline the major transmembrane ion transporters at the plasma membrane and tonoplast, along with their fundamental biophysical and regulatory characteristics, their basic kinetic descriptors, and relevant literature. In each case, these transporters divide between primary, energy-driven ion pumps—the plasma membrane and vacuolar H^{+}- and Ca^{2+}-ATPases, and the vacuolar H^{+}-PPase—H^{+}-driven solute pumps such as H^{+}-driven anion symporters (Meharg and Blatt, 1995), and ion channels including the slow vacuolar channel of the tonoplast identified by the *TPC1* gene of Arabidopsis (Peiter et al., 2005). Supplemental Appendix S1 includes descriptions of the buffering reactions for each compartment relating to H^{+}, Ca^{2+}, and Mal, and the encapsulated metabolic reactions for Mal and Suc. The sequential steps followed in the computation of *QT(t)* for each iteration during simulations are illustrated in the flow diagram of Figure 1.

#### The User Interface

The operational core of the OnGuard software comprises a GUI with a set of real-time displays of the standard, steady-state current-voltage (IV) curves at each membrane, tabular flux data, and a chart recorder. The software incorporates typical, multiple-document interface construction, with each document saved to disk representing a complete OnGuard model. The GUI is built on the Microsoft Windows platform and Microsoft Foundations Classes and is written in the C++ programming language and provides extensive, contextual help functionality. During simulations, IV displays provide visual representations of the kinetic characteristics for each of the various charge-carrying transporters in a scalable format similar to that of the Henry Electrophysiology Suite (Hills and Volkov, 2004). Current-voltage curves are updated in real-time along with the free-running membrane voltages at intervals specified by the operator. These IV representations offer immediate feedback on the evolving kinetic and thermodynamic features for each of the transporters and they enable quantitative assessment of the underlying causal relationships determining the temporal behavior of each ion flux. Illustrations of the IV displays are shown in Figure 2A and included with the results in the following article (Chen et al., 2012). Also shown, the tabular data summarizes the net fluxes of each ionic species at both membranes as well as the contents of each ionic species in the three compartments, vacuole, cytosol, and apoplast (Fig. 2B). Again, these data are available in real time and, optionally, can be logged separately in comma-delimited format that can be read by spreadsheet programs such as Microsoft Excel and SigmaPlot. Finally, OnGuard provides a running graphical display—an on-screen chart recorder—of the most important, user-selectable data (Fig. 2C).

The parameters of any model are accessed by means of a series of property pages, drop-down lists, and pop-up dialogue boxes (Fig. 3A). These pages include general buffering and environmental parameters such as turgor pressure/aperture relations and ion concentrations, the parameters specifying the various transporters at each membrane, parameters defining the characteristics for Suc and Mal metabolism, the light:dark cycle, and the time. Dialogue boxes for each transporter include operator-selectable controls that define the inherent biophysical properties of the transporter. For example, these controls enable access to the voltage dependence of the inward-rectifying K^{+} channel and its gating charge, its ion selectivity and conductance, as well as options to define regulatory parameters such as the kinetics for its inactivation by [Ca^{2+}]_{i} and its activation by extracellular H^{+} (Schroeder et al., 1987; Blatt et al., 1990b; Blatt, 1992; Thiel et al., 1993; Grabov and Blatt, 1997, 1999); similarly, for the various pumps and ion-driven transporters these controls provide access to the (pseudo-) rate constants for the carrier model, and the facility to define allosteric (regulatory) ligand and light sensitivities. For transporters that do not carry charge and for transport activities that are less well defined with respect to voltage or related kinetic properties—notably H^{+}-driven ion-exchange activities—simple, concentration-driven transport definitions give access to transport stoichiometries, maximum transport rates, and apparent *K*_{1/2} values as well as any allosteric ligand or other regulatory sensitivities.

## RESULTS

### Model Construction with OnGuard

The usual approach to formulating dynamic models of cellular homeostasis begins with the definition of an initial state—the Reference State—representing a physiological resting condition of the system under study (compare with Lew et al., 1979; Lew and Bookchin, 1986). Once a Reference State is established, the operator introduces one or more perturbations that represent new physiological, pathological, or experimental conditions to be explored and follows the response of all system variables as they evolve over time.

#### The Reference State Wizard

To establish a Reference State, we devised a software wizard that allows the operator to specify the underlying biophysical status of the system—standard ion concentrations in each compartment and voltages across the plasma membrane and tonoplast—and then query the model for net charge, driver ion (for plants, H^{+}), and solute fluxes. The wizard allows the operator to review and edit parameters such as membrane voltage and compartment solute composition, and to determine parameters for metabolic equilibrium in Suc and Mal synthesis and catabolism; most important, the wizard gives access to pages that compare the fluxes of ionic species across each membrane and permit their balancing by adjustment to the population(s) of transporters and, as necessary, to their underlying kinetic descriptors (see Fig. 3B). Thus, the operator is able to interrogate and adjust the subsets of transporters affecting the flux of each ion to satisfy the requirements of a Reference State.

Obviously such manipulations imply a knowledge of the (likely) unit density and/or limiting current amplitudes and their kinetic parameters for each transporter against which the operator can base judgment of biological validity. For most transporters of guard cells this information is available, often to within a factor of three and in many cases with an accuracy of a few percent of the mean for the parameter (see Supplemental Tables S3–S6). Even when information for some transporters was uncertain, we found their characteristics sufficiently constrained by the properties of other processes to ensure minimal indetermination. For example, with the H^{+}-ATPase as the sole pathway for H^{+} export across the plasma membrane, achieving charge and net H^{+} flux balance in the Reference State required a close match with H^{+} consumption through Mal catabolism and with the dominant H^{+} return pathways, and required similar balancing of the associated ion fluxes. Defining the H^{+}-Cl^{−} and H^{+}-K^{+} symporters included in our model, but for which there is limited information in guard cells (see below), required coordinating their relative contributions to this H^{+} balance. The H^{+}-coupling ratios of these transporters (H^{+}-Cl^{−} symport returns two H^{+} per charge while H^{+}-K^{+} symport returns 0.5 H^{+} per charge; Sanders et al., 1985; Blatt and Slayman, 1987) implies a 2:1 ratio in the activities of the two currents at the free-running membrane voltage to achieve charge balance in 1:1 ratio with H^{+} export via the H^{+}-ATPase. Additionally, H^{+}-coupled Cl^{−} influx required an equal efflux of Cl^{−} through the sum of the anion channels, and H^{+}-coupled K^{+} influx required an equal efflux of K^{+} through the outward-rectifying K^{+} channels. Finally, each of these components added a new current and therefore required balancing with opposing currents of the H^{+}-ATPase and K^{+} channels such that the total membrane current is zero at the free-running voltage. The consequence was to constrain each of the currents to within a narrow range of values relative to one another and within the constraints of the known densities and current amplitudes for the H^{+}-ATPase, K^{+}, and Cl^{−} channels characteristic of the guard cells (see Supplemental Tables S3 and S4).

In practice, a systematic approach dictated that fluxes be balanced first at the tonoplast and then at the plasma membrane. Flux adjustment at each membrane began with nondriver ions (K^{+}, Cl^{−}, Ca^{2+}, Mal^{2−}, Suc) and culminated with fine adjustments to the H^{+}-ATPases (and for the tonoplast, the H^{+}-PPase) densities. Transporter numbers were assessed against the available experimental data and the process was repeated if the approximation failed. We found it sufficient to bring all of the fluxes of each of the transported species in balance to within ±10^{−15} mol s^{−1}; thereafter trial simulations achieved a stable Reference State, generally within 8 to 10 h of simulation time, and maintained this condition even when allowed to run free over periods equivalent to many months.

#### The Choice of Transporters and Their Parameters

Most of the relevant quantitative kinetic information available comes from studies of the guard cells of *Vicia* (Supplemental Tables S3–S6). Therefore, we used this cell system as our starting point for model construction. Nonetheless, a complete set of parameters for all of the pumps, carriers, and channels is not available for any one species. A comparison of transport parameters listed in the tables shows a substantial degree of similarity between species, and even between cell types. For example, the fundamental gating characteristics for the plasma membrane K^{+} channels of *Vicia*, tobacco, and Arabidopsis guard cells (see Supplemental Table S4) include overlapping values for gating charge (δ) and half-activation voltage (*V*_{1/2}), as well as a similar sensitivities to extracellular [K^{+}], cytosolic pH, and [Ca^{2+}]_{i} (Dreyer and Blatt, 2009). Furthermore, Supplemental Table S3 highlights the close quantitative relationships between the plasma membrane H^{+}-ATPases of *Vicia* guard cells (Blatt, 1987a, 1988a), the giant alga *Chara* (Blatt et al., 1990a), and the fungus *Neurospora* (Gradmann et al., 1978; Sanders and Slayman, 1982; Blatt and Slayman, 1987). For these reasons we felt confident to borrow characteristics from other species—even other cell types—when information was lacking for *Vicia*, as necessary scaling the transporter for *Vicia* guard cells. For example, kinetic detail of H^{+}-coupled K^{+} symport is missing for guard cells. Nonetheless, an estimate of the probable symport current can be derived from measurements of K^{+} uptake against its electrochemical gradient in *Vicia* guard cells (see Supplemental Table S3; Blatt and Clint, 1989; Clint and Blatt, 1989) and a current similar to that documented for the H^{+}-K^{+} symport of Arabidopsis roots (Maathuis and Sanders, 1994). We could have modeled this transporter as a simple current source—that is, independent of membrane voltage—but we reasoned that the kinetic constraint of voltage was likely to be important for the dynamics of K^{+} and H^{+} balance, much as has been demonstrated in *Neurospora* (Blatt and Slayman, 1987). H^{+}-coupled K^{+} transport was originally described in *Neurospora* (Rodriguez-Navarro et al., 1986; Blatt and Slayman, 1987; Blatt et al., 1987) and, although the current is typically 10-fold greater in K^{+}-starved *Neurospora* than in the plant cells, the qualitative kinetic dependencies on external [K^{+}], pH, and voltage are similar. A relative weighting of reaction-kinetic constants is available for the carrier cycle of the H^{+}-coupled K^{+} symport in *Neurospora* (Blatt et al., 1987); hence, we used this weighting and the current amplitudes from Arabidopsis and estimated for *Vicia* as a guide in assigning values to the reaction-kinetic constants for the H^{+}-K^{+} symport in the guard cells. Much less kinetic information is available for plasma membrane Ca^{2+}-ATPases. Again, it was unrealistic to assume its voltage independence: Best estimates place the equilibrium voltage for the plasma membrane Ca^{2+} pumps near −200 mV with 1 mm [Ca^{2+}] outside and, hence, within the physiological voltage range. Instead, we borrowed the carrier cycle established for the H^{+}-ATPase, scaling it to estimates of the typical Ca^{2+} flux and measured values for the *K*_{m} with respect to [Ca^{2+}]_{i}, thereby ascribing a significant dependence on membrane voltage over much of the physiological range. Finally, to avoid some uncertainties in defining parameters we concatenated transport activities associated with Cl^{−} and NO_{3}^{−}, subsuming NO_{3}^{−} with Cl^{−} as a single, anionic species. Both anions contribute to the osmotic content of the guard cells, their relative presence depending to a large extent on availability, both are permeable through several of the anion channels present at the two membranes, and coupled transport for Cl^{−} and NO_{3}^{−} show similar thermodynamic and kinetic properties, to the extent that they are known (Sanders et al., 1989; Meharg and Blatt, 1995). By this maneuver we could draw on the available kinetic detail for H^{+}-coupled NO_{3}^{−} transport at both membranes (Supplemental Tables S3 and S5) and avoid redundancy with less-well-defined characteristics, notably for Cl^{−} transport at the plasma membrane.

A case-by-case summary of the reasoning behind the selection of each transporter will be found in Supplemental Appendix S2. The complete set of OnGuard parameters used for the modeling described in this article will be found in Supplemental Appendix S6 and is available for download with a demonstration of the OnGuard software (see www.psrg.org.uk). Our inclusive approach taken in developing the model described below expands the number of model parameters and, correspondingly, the model complexity. It is often a choice in modeling to use a skeleton construction to determine the minimum number of components, and associated parameters, needed to simulate a particular set of biological phenomena. In this case, however, the parameter space for the ensemble of transporters is exceptionally well constrained experimentally, generally to within a factor of three and frequently with an accuracy of a few percent (see Supplemental Tables S3–S6). This information, and the fundamental requirements for charge and ionic balance, placed narrow limits on the characteristics of the remaining transporters, even when their parameterization was less certain. As we demonstrate here, and in the companion article (Chen et al., 2012), a deterministic solution is readily found that recapitulates a very wide range of physiological behaviors. Nonetheless, the OnGuard software enables the user to generate and test skeleton models with equal ease.

#### Closed and Open Reference States

To test the parameter space of OnGuard models and their robustness, we used the results of initial simulations as a guide in defining a pair of states of the guard cell associated with the closed and open stomata, namely the Closed and Open Reference States. The Closed Reference State was envisaged as typical of stomata in the dark, in which the guard cells retained a baseline of osmotic load and a minimum of ion flux across the tonoplast and plasma membrane. For purposes of modeling, the currents of all primary pumps (H^{+}-ATPases, Ca^{2+}-ATPases, H^{+}-PPase) at the tonoplast and plasma membrane were assigned values of 5%, 10%, or 20% of their outputs in the Open Reference State, consistent with experimental estimates of the known light-stimulated activities of plasma membrane ATPases (Assmann et al., 1985; Shimazaki et al., 1986; Serrano et al., 1988; Goh et al., 1995, 1996; Kinoshita et al., 2001; Chen et al., 2012). Again, for the Closed Reference State we sought parameter sets such that the net fluxes of all solutes, Φ_{i}, across both tonoplast and plasma membrane were zero, that is ΣΦ_{i} = 0 for each membrane. We then used the same model parameters for the Open Reference State, allowing only the step up in primary pump currents, Suc and Mal synthesis. All other model variables were kept constant between these paired Reference States, and fine adjustments to the properties of individual transporters were introduced between simulations. Thus, Reference State pairing offered a convenient test of the capacity of the parameter ensemble to support stomatal movement across a range of common environmental variables.

As an initial test of the parameter ensemble, we explored model robustness by varying systematically the surface densities of individual transporters within bounds of ±1.3-, 2-, and 5-fold of their starting values as dictated by experimental knowledge and uncertainty. Figure 4 summarizes the results of simulations within the ±1.3-fold boundary reporting both macroscopic (aperture, osmotic content, total Ca^{2+}) and microscopic (pH_{i}, pH_{v}, [Ca^{2+}]_{i}) outputs. Displacements of each output were normalized to the corresponding value obtained with the starting parameters for purposes of comparison. In general, we found the model to be relatively insensitive to several of the predominant osmotic solute transporters at the plasma membrane. For example, a 1.3-fold increase or decrease in the population of R- (ALMT-) type anion channels (Keller et al., 1989; Meyer et al., 2010) resulted in <5% variation in any of the outputs, either in the closed or open states. Only outside the physiological limits for the channel population (see Supplemental Table S4) were variations in output returned that, in principle, might be detectable through biological experimentation (not shown). Similar sensitivities were recovered with variations in densities for the outward-rectifying K^{+} channel, and H^{+}-Mal symport at the plasma membrane, and for analogs of the TPK1, TPC1, FV, VCL, and VMAL channels at the tonoplast. These results demonstrated a considerable robustness despite the corresponding ranges in densities associated with these transporters. They also suggest a substantial functional tolerance for divergent transporter densities in the guard cell.

The OnGuard model proved substantially more sensitive to variations in the densities of transporters directly affecting [Ca^{2+}]_{i} and cytosolic pH (see Fig. 4). As central signaling and energetic elements, both [Ca^{2+}]_{i} and cytosolic pH affect an extended network of transporters at both membranes (Blatt, 2000; Schroeder et al., 2001) and, thus, were anticipated to impose greater restrictions on performance. Substantially greater relative variations in macro- and/or microscopic outputs were recovered, even with 1.3-fold changes in densities of the H^{+}-ATPase, Ca^{2+}-ATPase, and Ca^{2+} channel at the plasma membrane, and for the Ca^{2+}-ATPase and Ca^{2+} channel at the tonoplast. Notable sensitivities, especially in cytosolic and/or vacuolar [H^{+}] were evident also to changes in the densities of the H^{+}-Cl^{−} and H^{+}-K^{+} symporters at the plasma membrane and the CLC, NHX, and CAX antiporters at the tonoplast. Aperture and osmotic content sensitivities to H^{+}-Cl^{−} symport density is consistent with its central role in anion uptake and, for reasons noted above, variations in either the H^{+}-Cl^{−} symport, H^{+}-K^{+} symport, or H^{+}-ATPase densities had the greatest impact in the Closed Reference State, especially on pH_{i} balance. The roles for the endosomal cation and anion antiporters in regulating pH are broadly consistent with observation documented in the literature (Pardo et al., 2006; Padmanaban et al., 2007; Braun et al., 2010; Smith and Lippiat, 2010; Weinert et al., 2010).

For purposes of comparison, we also varied the biophysical and kinetic parameters for transporters with appreciable effects on both the macro- and microscopic parameters in Figure 4. As examples, Figure 5 summarizes the outputs from simulations with the plasma membrane Ca^{2+} channel and the tonoplast Ca^{2+}-ATPase, both of which contribute to the Ca^{2+} circuit of the model. For the Ca^{2+} channel, we varied by ±2-fold each of the key gating parameters affecting its voltage sensitivity—*V*_{1/2} and δ that define the voltage-yielding half-maximal conductance and the sensitivity of the channel gate to changes in voltage, respectively—and the apparent *K*_{Ca} and associated Hill coefficient, *h*_{Ca}, for channel inactivation by [Ca^{2+}]_{i}. For the Ca^{2+}-ATPase, the same strategy was applied to the apparent *K*_{Ca} and *h*_{Ca} determining [Ca^{2+}]_{i}-dependent activation. We also varied by ±2-fold the relative magnitudes of the reaction-kinetic constants *k*_{12}^{o} and *k*_{21}^{o} for the Ca^{2+}-ATPase, while maintaining their ratio constant to avoid any thermodynamic bias. The constants *k*_{12}^{o} and *k*_{21}^{o} define the voltage dependence of the transport cycle, and their magnitudes relative to the reaction-kinetic constants for the rest of the cycle determines the range of voltages over which transport flux is voltage sensitive (Hansen et al., 1981). The fundamental parameters for the Ca^{2+} channel are constrained through direct experimental analysis to values well within this ±2-fold range (Hamilton et al., 2000, 2001; Pei et al., 2000; Köhler and Blatt, 2002; Supplemental Table S4), whereas this is not the case for the Ca^{2+}-ATPase. So the comparison also offered a useful context for the parameterization of the pump. As expected, altering any of the parameters, either for Ca^{2+} channel or the Ca^{2+}-ATPase, affected all macroscopic outputs as well as [Ca^{2+}]_{i}, notably in the Open Reference State. Less obvious was their influence on cytosolic and vacuolar pH, most evident in each case in the Closed Reference State; this connection between Ca^{2+} and pH arises in part through their overlaps in [Ca^{2+}]_{i}-mediated regulation of transport across the tonoplast, and we return to the relationship between [Ca^{2+}]_{i}- and pH-dependent transport in the following article (Chen et al., 2012). These details aside, the similarities between pump and channel in scale and distribution of the outputs indicates a similar degree in sensitivity to variations in parameters for the Ca^{2+}-ATPase.

#### Testing Reference State Behavior with Environmental Challenge

Ultimately, the validity of any homeostatic model lies in its ability to recapitulate physiological behavior. There is much quantitative experimental detail relating to the physiological and macroscopic responses of guard cells to their ionic environment (see Supplemental Tables S1 and S2; Willmer and Fricker, 1996). Therefore, we varied external solute concentrations for both the closed and open reference states to compare the effects on a range of outputs. Simulations were run with KCl concentrations ranging from 1 to 30 mm, with CaCl_{2} concentrations from 0.3 to 3 mm, and pH values from 6.0 to 7.0. The outputs demonstrated that the model successfully recapitulates stomatal behavior under each of these conditions. For example, increasing KCl concentration resulted in a progressive rise in guard cell volume and turgor in the open reference state, with stomatal aperture rising from around 8 μm in 1 mm KCl to almost 15 μm in 20 to 30 mm KCl (Fig. 6A). In the closed reference state, guard cell volume and turgor showed a slight decrease with increasing KCl concentration and stomatal aperture declined below 4 μm at the higher KCl concentrations, consistent with the increased osmotic load outside. Complementary and physiologically sensible outputs were evident for every other experimentally measurable variable, including membrane voltage (Fig. 6B), cytosolic and vacuolar K^{+}, Cl^{−}, and Mal concentrations (Fig. 7, A and B), pH_{i}, pH_{v}, and [Ca^{2+}]_{i} (Fig. 8, A and B).

These outputs also demonstrate a number of less obvious, but equally important trends, again recapitulating well-documented experimental data (Willmer and Fricker, 1996; see also Raschke and Schnabl, 1978; Van Kirk and Raschke, 1978; Blatt, 1987b, 2000; Lohse and Hedrich, 1992; Talbott and Zeiger, 1996; Dodd et al., 2005; Wang and Blatt, 2011). Among these, plasma membrane voltage showed a steeper dependence on KCl concentration in the Closed than in the Open Reference State (Fig. 6B); voltages in both states, and aperture in the Open Reference State, declined with CaCl_{2} concentration (Fig. 6B) while [Ca^{2+}]_{i} rose (Fig. 8B); increasing KCl concentration was accompanied by biphasic, and opposing changes in Cl^{−} and Mal in the vacuole in the Closed Reference State; and decreasing pH outside promoted cytosolic K^{+} concentration, vacuolar K^{+}, and Mal concentrations (Fig. 7, A and B), and stomatal aperture (Fig. 6A) in the open reference state, but had little influence on membrane voltage, pH_{i}, or pH_{v} (Figs. 6B and 8A). Finally, we note that, in every case, [Ca^{2+}]_{i} was elevated in the open compared with the Closed Reference State (Fig. 8B). The bases for these observations are fully explained by emergent interactions between the several transporters, the dominance of the plasma membrane H^{+}-ATPase in the Open Reference State, and the influence of membrane voltage on transport, especially across the plasma membrane. We explore in depth these interactions, and their consequences, in the following article (Chen et al., 2012).

## DISCUSSION

A major challenge in constructing any quantitative kinetic model of transport and homeostasis is to integrate the often substantial body of data in a systematic representation of the cell. Although basic physicochemical relationships are all simple quantitative relations easily incorporated in a quantitative description of transport, the recursive nature of charge transport presents a challenging dimension to any modeling effort that generally defies simple analytical solution. We have developed an open structured approach to software construction that should prove widely applicable to modeling cellular transport and homeostasis. HoTSig provides a standardized library that incorporates all of the major sets of equations and descriptors for transport across biological membranes, and a minimum set of equations defining the metabolism of the major organic solutes. The OnGuard software, constructed from this library, incorporates a GUI for real-time feedback on the individual transport and homeostatic processes to model the physiological behavior of stomatal guard cells; it incorporates a set of empirically defined equations to relate the output of solute content to cell volume, turgor, and stomatal aperture. Additionally, we introduce the Reference State Wizard and outline its use in defining a starting point for in silico experimentation with sensible outputs for all known variables. Finally, we show that an OnGuard model, with a realistic ensemble of transporters, is able to recapitulate the known characteristics of guard cells and stomata in the face of common experimental manipulations.

### A Quantitative Modeling Approach

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. Even so, the subject is often perceived as difficult or inaccessible. Most modeling efforts have been implemented on a case-by-case basis without a standardized format and, as a consequence, mastering all of the individual metabolic, transport, and buffering mechanisms presents a challenge, especially as much of the quantitative data relating to membrane transport are to be found in older literature and in formats that are not readily accessible except to the specialist. Utilities such as the Virtual Cell (Loew and Schaff, 2001) provide for modeling intracellular events that encompass reaction-diffusion processes in arbitrary geometries, but offer limited flexibility in defining underlying behaviors, for example as dictated by specific transport equations. These limitations are addressed in our development of the HoTSig approach. The software contains an expandable library for transporter kinetics, chemical buffering, and capacity for macromolecular binding and metabolic reactions, all accessible to input and modification by the operator. This open structure should make HoTSig adaptable for wide variety of single-cell systems.

The HoTSig library and OnGuard software greatly expands on earlier efforts of our own (Lew et al., 1979; Lew and Bookchin, 1986; Lew, 1991; Gradmann et al., 1993) and others (Loew and Schaff, 2001; Hunter and Nielsen, 2005; Shabala et al., 2006) in providing real-time feedback for the operator as well as a detailed output log. We adapted a number of generalized utilities that greatly reduce computational load and time, notably the Newton-Raphson chord approximation to solutions of implicit equations (Lew and Bookchin, 1986), and we introduced a look-ahead utility to adjust time increments to the pending dynamics of the model (Supplemental Appendices S3 and S4). We have placed considerable emphasis on an intuitive GUI for the entry of initial conditions, simulated perturbations, and for tabular and graphical displays to make simulation outputs readily apparent to the operator. Unlike the earlier efforts, the OnGuard graphical interface provides detailed flux and kinetic information in tabular and chart-recorder formats as well as real-time electrophysiological displays with the component and ensemble current-voltage relations at each membrane. We anticipate users of the OnGuard software will start with the defaults specified by the parameter set available for download with the demonstration (see www.psrg.org.uk), thereafter customizing individual transporters and/or environmental conditions to address specific physiological questions, species variations, and the consequences of genetic manipulations. Interpreting the current-voltage outputs requires some knowledge of electrophysiology; nonetheless, the tabular and chart-recorder formats offer useful reference points for comparison in real time.

To make this process as rapid and intuitive as possible, each descriptor in a model is editable and parameters are accessible on the fly during modeling sessions. As editable modules, these descriptors serve as phenomenological black boxes to be opened, or reduced, whenever the internal workings become a desirable or necessary part of a modeling project. From the point of view of a modeler seeking to understand how the homeostatic variables of a given system respond to a physiological or experimental perturbation, the only relevant biology is encapsulated by how one model variable is connected to another. It is frequently the case that this phenomenology is directly accessible to experiment, whereas the underlying mechanistic details are not, or are accessible only qualitatively. By fitting the experimental data with a mathematical relation, they are safely placed within a black box that represents a parameterized module with adjustable levels of resolution (Endy and Brent, 2001) and may be opened if, and when the included regulatory pathway becomes a target of study for the purposes of modeling. In effect, this phenomenological approach serves as a place holder for unknowns and as a means to reduce complexity and computational burden.

We anticipate that the HoTSig library will find general applications in exploring cell systems, in addition to guard cells, for which sufficient biophysical and kinetic detail for transport is now available. For example, it should prove useful in exploring the physiology of the plant root epidermis and its interaction with the soil environment. Root epidermal cells, including root hairs, function in the uptake of essential mineral nutrients from the soil, but also provide entry pathways for toxic ions, including heavy metals and Na^{+}, that arise through intensive irrigation and pollution, and present a major challenge to modern agriculture. Aspects of Na^{+} toxicity as well as pH, Ca^{2+}, and K^{+} transport oscillations have found their way into computational assessments of ion transport in roots (Amtmann and Sanders, 1999; Shabala et al., 2006), but will now benefit from quantitative analysis incorporating appropriate reference states that help minimize indetermination and validate predictive power. Finally, the OnGuard software is equally applicable to guard cells of species other than *Vicia*. As we noted above (see “The Choice of Transporters and Their Parameters”), there exists substantial quantitative similarity between the guard cells of a number of plant species, including *Vicia*, *Nicotiana*, and Arabidopsis. In large measure, adapting the model we resolved for *Vicia* to other species requires only an accounting for differences in cell geometry and the relationships between cell surface area, volume, and turgor pressure.

### Modeling Guard Cell Dynamics

In practice, the OnGuard software returned an ensemble of model parameters that, in simulations of limiting stomatal behavior, were robust and recapitulated the predominant characteristics of guard cell physiology with a bare minimum of intrinsic assumptions. Once resolved with the Reference State Wizard, we found little difficulty in establishing a model using parameters within experimentally determined constraints. Often Monte Carlo methods are used to determine the best solution(s) to simulations that must deal with an extended, *n*-dimensional parameter space when a deterministic solution cannot be found (Moskowitz and Caflisch, 1996; Robert and Casella, 2011). However, for all but a few of the guard cell transporters, the parameter space is exceptionally well constrained experimentally, generally to within a factor of three and frequently with an accuracy of a few percent (see Supplemental Tables S3–S6). This information, and the fundamental requirements for charge and ionic balance, placed narrow limits on the characteristics of the remaining transporters, even when their parameterization was less certain.

We found this model recapitulated the guard cell in the closed (dark) and open (light) states of the stomata with the single assumption of a defined, light-dependent increase in the activities of all ion-translocating ATPases and pyrophosphatase at the plasma membrane and tonoplast, and in Suc and Mal synthesis (Willmer and Fricker, 1996; Shimazaki et al., 2007). These paired reference states showed substantial intrinsic stability in the face of systematic changes to a range of parameters, including transporter densities at both membranes and the gating and regulatory characteristics of channels and pumps, the densities of which had the greatest effect on model outputs (Figs. 4 and 5). Perturbations in these parameters resulted primarily in differences in the dynamic range of macroscopic outputs, on [Ca^{2+}]_{i}, pH_{i}, and pH_{v}. Nonetheless, these outputs compare favorably with the literature (see Supplemental Tables S1 and S2; Figs. 6–8; Raschke, 1979; MacRobbie and Lettau, 1980; Hill and Findlay, 1981; Zeiger, 1983; Davies and Jones, 1991; Willmer and Fricker, 1996; Blatt, 2000; Schroeder et al., 2001; Shimazaki et al., 2007). Furthermore, the model successfully recapitulated stomatal characteristics (see Figs. 6–8) over a wide range of extracellular ion concentrations and pH for which there is substantial experimental documentation. These results fulfill a set of minimum criteria in validating the model and they also offer unexpected perspectives on the consequences of these manipulations, especially in the relationship between [Ca^{2+}]_{i}, cytosolic, and vacuolar pH. In the following article (Chen et al., 2012) we review these, and additional interactions, underlining the power of the OnGuard model in providing (otherwise counterintuitive) explanations for published experimental data with a minimum of underlying assumptions.

## MATERIALS AND METHODS

Details of the model assembly and computational components will be found in the Supplemental Appendices S1 to S4.

### Supplemental Data

The following materials are available in the online version of this article.

**Supplemental Figure S1.**Generalized four-state reaction kinetic cycle.**Supplemental Figure S2.**Malic acid (de-) protonation reactions.**Supplemental Figure S3.**Global metabolic reactions for synthesis and degradation of sucrose and malic acid.**Supplemental Table S1.**Basic biophysical parameters of*Vicia*stomatal guard cells.**Supplemental Table S2.**Compartmental contents of*Vicia*stomatal guard cells.**Supplemental Table S3.**Predominant plasma membrane pumps and carriers.**Supplemental Table S4.**Predominant plasma membrane ion channels.**Supplemental Table S5.**Predominant tonoplast pumps and carriers.**Supplemental Table S6.**Predominant tonoplast channels.**Supplemental Appendix S1.**The HoTSig Platform and OnGuard software.**Supplemental Appendix S2.**A brief summary of transporter selections.**Supplemental Appendix S3.**Implementing Newton-Raphson approximations.**Supplemental Appendix S4.**Determing Δt.**Supplemental Appendix S5.**Abbreviations.**Supplemental Appendix S6.**RCF1 model parameters.

## Acknowledgments

We thank Simon Rogers,Yizhou Wang, and Christopher Grefen (University of Glasgow) for comments during the development of the software and model.

## Footnotes

The author responsible for distribution of materials integral to the findings presented in this article in accordance with the policy described in the Instructions for Authors (www.plantphysiol.org) is: Michael R. Blatt (michael.blatt{at}glasgow.ac.uk).

↵1 This work was supported by the UK Biotechnology and Biological Sciences Research Council (grant nos. BB/F001630/1, BB/F001673/1, and BB/H024867/1 to M.R.B.)

↵2 These authors contributed equally to the article.

↵3 Present address: School of Science and Health, University of Western Sydney, Hawkesbury Campus, Richmond, NSW 2753, Australia.

↵4 These authors contributed equally to the article.

↵[OA] Open Access articles can be viewed online without a subscription.

↵[W] The online version of this article contains Web-only data.

### Glossary

- Mal
- malate
- HoTSig
- Homeostasis, Transport and Signalling
- GUI
- graphical user interface

- Received March 16, 2012.
- Accepted May 20, 2012.
- Published May 25, 2012.