Global Sensitivity Analysis of OnGuard Models Identifies Key Hubs for Transport Interaction in Stomatal Dynamics1[CC-BY]

Incorporating global sensitivity analysis to mechanistic models of stomatal guard cells highlights previously unexpected interactions in transport driving stomatal aperture. The physical requirement for charge to balance across biological membranes means that the transmembrane transport of each ionic species is interrelated, and manipulating solute flux through any one transporter will affect other transporters at the same membrane, often with unforeseen consequences. The OnGuard systems modeling platform has helped to resolve the mechanics of stomatal movements, uncovering previously unexpected behaviors of stomata. To date, however, the manual approach to exploring model parameter space has captured little formal information about the emergent connections between parameters that define the most interesting properties of the system as a whole. Here, we introduce global sensitivity analysis to identify interacting parameters affecting a number of outputs commonly accessed in experiments in Arabidopsis (Arabidopsis thaliana). The analysis highlights synergies between transporters affecting the balance between Ca2+ sequestration and Ca2+ release pathways, notably those associated with internal Ca2+ stores and their turnover. Other, unexpected synergies appear, including with the plasma membrane anion channels and H+-ATPase and with the tonoplast TPK K+ channel. These emergent synergies, and the core hubs of interaction that they define, identify subsets of transporters associated with free cytosolic Ca2+ concentration that represent key targets to enhance plant performance in the future. They also highlight the importance of interactions between the voltage regulation of the plasma membrane and tonoplast in coordinating transport between the different cellular compartments.

The physical requirement for charge to balance across biological membranes means that the transmembrane transport of each ionic species is interrelated, and manipulating solute flux through any one transporter will affect other transporters at the same membrane, often with unforeseen consequences. The OnGuard systems modeling platform has helped to resolve the mechanics of stomatal movements, uncovering previously unexpected behaviors of stomata. To date, however, the manual approach to exploring model parameter space has captured little formal information about the emergent connections between parameters that define the most interesting properties of the system as a whole. Here, we introduce global sensitivity analysis to identify interacting parameters affecting a number of outputs commonly accessed in experiments in Arabidopsis (Arabidopsis thaliana). The analysis highlights synergies between transporters affecting the balance between Ca 2+ sequestration and Ca 2+ release pathways, notably those associated with internal Ca 2+ stores and their turnover. Other, unexpected synergies appear, including with the plasma membrane anion channels and H + -ATPase and with the tonoplast TPK K + channel. These emergent synergies, and the core hubs of interaction that they define, identify subsets of transporters associated with free cytosolic Ca 2+ concentration that represent key targets to enhance plant performance in the future. They also highlight the importance of interactions between the voltage regulation of the plasma membrane and tonoplast in coordinating transport between the different cellular compartments.
Stomata form the major pathway for CO 2 entry across the leaf epidermis for photosynthesis and for water loss by transpiration from the leaf tissues. Pairs of guard cells surround the stomatal pore and regulate the aperture. These cells balance the demand for CO 2 with that for water conservation. Guard cells expand and contract to open and close the pore. They take up and lose solutes, notably K + and Cl 2 , which, together with the synthesis and metabolism of organic solutes, especially malate, provide the osmotic driving force for these changes in aperture (Kim et al., 2010;Roelfsema and Hedrich, 2010;Lawson and Blatt, 2014). Thus, membrane transport comprises the principal set of effectors of a regulatory network that ensures the homeostatic control of the guard cell for stomatal aperture.
Environmental signals, notably light, CO 2 , water availability, and the hormone abscisic acid, affect this network, modulating transport and solute accumulation. Research at the cellular and molecular levels has focused on these inputs and their contributions to stomatal movements. A large body of work has highlighted Ca 2+ -independent and Ca 2+ -dependent signaling, the latter including elevations in free cytosolic Ca 2+ concentration ([Ca 2+ ] i ), protein kinase and phosphatase activities, which inactivate inward-rectifying K + channels and activate Cl 2 (anion) channels, as well as the changes in cytosolic pH that promote the outwardrectifying K + channels and solute loss (Keller et al., 1989;Blatt et al., 1990;Thiel et al., 1992;Lemtiri-Chlieh and MacRobbie, 1994;Blatt, 1998, 1999;for review, see Blatt, 2000;Hetherington and Brownlee, 2004;Kim et al., 2010;Lawson and Blatt, 2014;Jezek and Blatt, 2017).
Despite this extensive body of knowledge about the individual transporters and their regulation, relating the transport capacity of guard cells to stomatal movements in quantitative mechanistic terms poses a number of difficulties. Most important, the physical requirement in transport for charge to balance across each membrane means that the transport of any one ionic species is necessarily joined to that of all others across the same membrane (Jezek and Blatt, 2017). As a result, the capacity for flux through most transporters is rarely limiting, especially through the individual ion channels that facilitate K + , Cl 2 , and malate flux. Instead, transport is generally limited by the balance in charge flux across the plasma and tonoplast membranes. As a corollary, manipulating solute flux through any one transporter inevitably affects this balance in transport and, thereby, directly affects other transporters at the same membrane, often with unforeseen consequences.
Systems modeling offers one approach for overcoming these difficulties to quantify the mechanics of stomatal movements and to predict and interpret the effects of experimental manipulations. It enables the detailed knowledge available for the individual transporters to be reconstructed within the physiological framework of the cell, fully constrained by fundamental physical laws and the known kinetic relationships, ligand binding, and related regulatory properties for each transporter. The development of the OnGuard platform for modeling guard cells Hills et al., 2012; freely available at www.psrg.org.uk) has proven most successful to date, demonstrating true predictive power in uncovering previously unexpected and emergent behaviors of stomata (Wang et al., 2012(Wang et al., , 2014aBlatt et al., 2014;Minguet-Parramona et al., 2016).
Despite these successes, resolving models with the OnGuard platform to date has required manual exploration with over 200 parameters that define the characteristics of the dominant ion transporters and organic solute metabolism. Roughly 85% of these parameters are known experimentally to within a factor of 3, and half of these are known accurately to within a margin of 5% to 10% of the parameter value. Nonetheless, refining OnGuard models has remained an arduous task, often demanding many weeks to explore the parameter space across each selection of experimental variables. Furthermore, this process captures little formal information about the connections between subsets of parameters that determine transport, defining its most interesting characteristics that emerge from interactions within the system as a whole. We have now implemented a semiautomated approach to estimate the importance of transporters and the interactions between them in OnGuard, making use of global sensitivity analysis methods to determine the parameter subsets that, together, account for the largest proportions of variance in several model outputs that are commonly reported in experimental studies. Our findings highlight the emergent synergy inherent in transport and the several hubs of interaction between transport functions. This information is likely to help in selecting potential targets for enhancing plant performance in the future. It also demonstrates the importance of the trans-network of interactions between the plasma membrane and tonoplast that coordinates transport between the different cellular compartments. Figure 1 shows the standard Arabidopsis model outputs for stomatal aperture and [Ca 2+ ] i over the diurnal cycle with stepwise light transitions. Also shown is an exploded view of the [Ca 2+ ] i time course at the end of the daylight period and following the transition to dark when the stomata close Wang et al., 2012). The oscillations in [Ca 2+ ] i are coupled to oscillations in membrane voltage of several minutes duration, with transients in [Ca 2+ ] i from values near 200 nM to maxima in excess of 1 mM. This entire sequence of oscillations is stably repeated with each diurnal cycle for all Figure 1. Macroscopic outputs of the OnGuard Arabidopsis model. Outputs resolved over a diurnal cycle with stepped transitions between dark and light (the dark period is indicated by the black bar above) and with 10 mM KCl, 1 mM CaCl 2 , and pH 6.5 outside. The full set of model parameters and initializing variables can be found in Supplemental Appendix S1 and may be downloaded with the OnGuard software at www.psrg.org.uk. A, Model output of stomatal aperture for the diurnal period 22 to 17 h relative to the start of the diurnal cycle. Secondary parameters derived from this output and used in the Sobol analysis are indicated for the initial increase (rise) and decrease (fall) in stomatal aperture, the minimum (min) and maximum (max) apertures, and the dynamic range of the aperture. B, Model output of [Ca 2+ ] i . Secondary parameters derived from this output and used in the Sobol analysis are indicated for the number of [Ca 2+ ] i oscillation cycles and the oscillation period. In the latter case, the median of the oscillation periods was used.  Figure 1, with the Sobol index scale (right). The full list of transporter identities and functions can be found in Jezek and Blatt (2017) and Hills et al. (2012). Note that the highest indices across the vertical and horizontal associated in each case with parameters for the endomembrane Ca 2+ release channel VCa in and the Ca 2+ sequestration pump VCa-ATPase. The high indices here indicate that parameters for virtually every other transporter interact synergistically with those of the VCa in and VCa-ATPase. In affecting the maximum aperture (A), this interdependence is highest for the K d values for Ca 2+ of the VCa in and VCa-ATPase. Substantial interdependence also is seen between the K d for Ca 2+ of the VCa in and the TPK K + channel of the tonoplast and parameters of the plasma membrane Ca 2+ channel (Ca in ), the Cl 2 channel (SLAC), and the Ca 2+ export pump (Ca-ATPase). Note, too, that a similar set of interdependencies is seen in the aperture dynamic range (C), albeit primarily with the VCa-ATPase, and that the interdependencies affecting the minimum aperture (B) are generally of lower index values and include hubs associated with the plasma membrane Ca 2+ -ATPase and the tonoplast Ca 2+ /H + antiporter CAX. known model parameter sets Minguet-Parramona et al., 2016). The major peaks of voltage and [Ca 2+ ] i occur in antiparallel fashion and show a mean frequency of 1.9 mHz (=8.9 min) for the standard Arabidopsis parameter set (Minguet-Parramona et al., 2016), thus closely matching previously published data for [Ca 2+ ] i and voltage oscillations in these and other guard cells (Blatt and Armstrong, 1993;McAinsh et al., 1995;Grabov and Blatt, 1998;Allen et al., 2001). These oscillations also are coupled to the periodic, almost stepwise decreases in stomatal aperture evident in Figure 1. For purposes of analysis, we derived the secondary outputs, as indicated in Figure 1, that are commonly reported in experimental studies.

RESULTS AND DISCUSSION
Global variance-based analysis with Sobol indices (Sobol, 2001) addresses problems inherent to work with model parameters that may range across many orders of magnitude and gives a detailed, high-level view of parameter interactions. The method places emphasis on the variance in parameter values and thereby obviates issues of parameter weighting. In other words, if a transporter becomes a target of interest to improve plant behavior, the variation needed will be highlighted independent of the parameter value itself. The Sobol analyses of interacting parameters for minimum and maximum stomatal aperture and the associated dynamic range are shown in Figure 2. Here, and in the following figures, the Sobol indices are shown as heat maps to highlight those transporter parameters with the strongest interdependencies: in other words, those parameters that affect the output synergistically (see "Materials and Methods"). A complete list of transporter definitions will be found in Hills et al. (2012) and Jezek and Blatt (2017), and the acronyms are summarized in Table I. We discounted trivial interactions for any one transporter, such as that between the Hill coefficient and K d for Ca 2+ of the endomembrane Ca 2+ release channel (VCa in ), to focus on interdependencies between transporters.
A review of the maximum aperture interactions ( Fig. 2A) showed a substantial dependence on the K d for Ca 2+ of VCa in and its interactions with several other transporters, most strongly with the K d for Ca 2+ of the endomembrane Ca 2+ sequestration pump VCa-ATPase but also with parameters for the plasma membrane Ca 2+ channel Ca in , the Ca-ATPase, the SLAC anion channel, and the Hill coefficient for Ca 2+ control of the tonoplast TPK K + channel. These interactions yielded Sobol index values from 0.28 to 0.42, supporting the conclusion of an interdependence between these parameters and the importance of endomembrane Ca 2+ circulation in determining the maximum aperture. They also highlight the substantial interaction between transport at the plasma membrane and tonoplast that was identified previously and associated with Ca 2+ Wang et al., 2012;Minguet-Parramona et al., 2016).
The minimum aperture (Fig. 2B) also showed an interdependence between the K d values for Ca 2+ of VCa in and the VCa-ATPase, albeit weaker. However, nearly equivalent interdependencies also were evident for the minimum aperture between the K d for Ca 2+ VCa-ATPase with parameters for the plasma membrane Ca-ATPase and, to a lesser extent, with several other plasma membrane and tonoplast transporters, including the TPK K + channel and the CAX Ca 2+ /H + antiporter. Finally, a comparison with these heat maps showed that the aperture dynamic range (Fig. 2C) was most sensitive to the interactions between the K d values Table I. Transporter acronyms used in Figures 2 to 4 Further details of the transporters, their characteristic parameters, and their inclusion in the OnGuard model can be found in Jezek and Blatt (2017) and Hills et al. (2012). Parameters examined for each transporter listed in the figures are for the transporter number (Transporter), its primary K d (K d ), and the corresponding Hill coefficient (Hill). for Ca 2+ of VCa in and the VCa-ATPase, but unlike the maximum aperture, the dominant interactions otherwise were between the K d for Ca 2+ of the VCa in with parameters of the plasma membrane Ca 2+ channel, the SLAC anion channel, and the Ca-ATPase. Figure 3 shows the heat maps for parameter interactions affecting the initial rates of stomatal opening and closing. The most notable interactions in stomatal opening were between parameters defining the plasma membrane Ca-ATPase and the ALMT anion channel and in stomatal closing between endomembrane VCa in and the plasma membrane H-ATPase. The Sobol indices in these instances ranged between approximately 0.08 and 0.15, appreciably lower than those associated with minimum and maximum apertures. Nonetheless, these interactions highlight the importance of membrane voltage, and the transporters that are most strongly affected by its variation, in driving the rates of solute accumulation and loss (Thiel et al., 1992;Minguet-Parramona et al., 2016).
Previous work highlighted the importance in stomatal closure of oscillations in membrane voltage and  (McAinsh et al., 1990;Fricker et al., 1991;Gradmann et al., 1993;Grabov and Blatt, 1999), and they predicted the dependence of oscillation frequency on transport parameters such as those defining the Ca 2+ sensitivity of the SLAC anion channel (Chen et al., 2010;Minguet-Parramona et al., 2016).
Sobol analyses for [Ca 2+ ] i oscillation frequency and number are shown in Figure 4. As might be anticipated from this previous study of [Ca 2+ ] i oscillation structure, the analysis for oscillation frequency (Fig. 4A) showed substantial interdependencies between parameters for  (2017) and Hills et al. (2012). Note the difference in Sobol index scales. The highest indices across the vertical and horizontal are associated with the rate of closure (B) and interacting parameters for the endomembrane VCa in channel and the primary plasma membrane H + -ATPase. However, an appreciable hub also is associated with the tonoplast TPK channel. For the initial rate of stomatal opening (A), hubs of interdependence are associated with parameters for the plasma membrane Ca 2+ -ATPase and the ALMT anion channel, which plays a key role in the control of membrane voltage and stomatal closing (Meyer et al., 2010;Wang et al., 2014a;Minguet-Parramona et al., 2016;Jezek and Blatt, 2017).

VCa in and virtually every other transporter examined.
The strongest interdependencies were observed with parameters for the VCa-ATPase and yielded Sobol indices near 0.2 and 0.32; however, all of the other transporter parameters yielded appreciable Sobol indices, with values around 0.16 to 0.18. Minguet-Parramona et al. (2016) also noted previously that the oscillation number was subject to the dynamic range of stomatal apertures for a given set of conditions: in other words, the number of cycles needed to transit from the open to the closed state. Thus, it was not surprising that we found weak interdependencies only, largely in the range of 0.02 to 0.04 and including the interaction between VCa in and the VCa-ATPase. Unexpectedly, however, this analysis also predicted at the hub of these interactions both the VCa-ATPase and the TPK channel. Furthermore, the strongest interdependence was predicted between parameters defining TPK and VCa in . This interaction is particularly noteworthy. Past studies showed that the Arabidopsis tpk1 mutant exhibits a much reduced rate of stomatal closure, even though the K + content of the mutant is seemingly unaffected (Gobert et al., 2007). Given the emergent interdependence of VCa in with the TPK channel, we predict that the effects of the tpk1 mutant are mediated through their interaction with endomembrane Ca 2+ release and their impact on [Ca 2+ ] i during stomatal closure.
Complex biological processes are defined by nonlinear relationships and are never equal to the sum of the properties of their components. New properties emerge from the interactions between components that cannot be predicted a priori, even with knowledge of the underlying processes behind each component in isolation. These so-called emergent properties (Novikoff, 1945) are amply represented by the behavior of stomata, such as in the SLAC1 Cl 2 channel, nominally associated with stomatal closure, and its unexpected effects on the K + channel activities in stomatal opening identified through the OnGuard modeling platform. In this case, Wang et al. (2012) predicted a connection in the Arabidopsis slac1 mutant between the Cl 2 and K + channels at the plasma membrane that was subsequently validated experimentally. OnGuard modeling showed how slac1 slowed K + uptake and stomatal opening, even though the Cl 2 channel contributes directly only to solute loss and stomatal closure. Indeed, the OnGuard platform Hills et al., 2012; freely available at www.psrg.org.uk) has proven most successful to date, demonstrating true predictive power in uncovering a number of previously unexpected and emergent behaviors of guard cells Wang et al., 2014a;Minguet-Parramona et al., 2016). These studies illustrate how quantitative modeling is essential as an approach to physiology that otherwise confounds intuitive understanding and leads to misinterpretations. To date, however, reviews of the interactions inherent to this complex system of transport processes have been limited. Wang et al. (2014a) carried out a local sensitivity analysis, extracting the key macroscopic outputs, their dependencies on a selection of parameters for individual transporters, varying these parameters one at a time. Such first-order analysis, although enlightening, is not well equipped to capture global information about the interdependencies between transporters. This information is essential for any critical evaluation of how the relevant transporters are coordinated and which characteristics are most important for such coordination.
Our screen here for parameter interdependencies underscores the central importance of a number of transport interaction hubs. We selected for analysis the subset of transporters with known and direct connections to [Ca 2+ ] i . These transporters comprise some 70% of the combined transport activities at the plasma membrane and tonoplast. The majority do not transport Ca 2+ or affect [Ca 2+ ] i per se, although they respond to [Ca 2+ ] i (Jezek and Blatt, 2017). Not surprisingly, the analysis highlights interactions between transporters affecting the balance between Ca 2+ sequestration and Ca 2+ release pathways, notably those associated with internal Ca 2+ stores and their turnover. The K d values for Ca 2+ of the endomembrane Ca 2+ channel VCa in and the VCa-ATPase are predicted to be central for interactions affecting minimum and maximum aperture (Fig. 2), the initial rate of stomatal closure (Fig. 3), and the median period and number of [Ca 2+ ] i oscillations (Fig. 4). Each of these overlapping hubs categorizes substantial synergies between parameters across virtually every transporter at both the plasma membrane and tonoplast and, thereby, underscores the intrinsic functional network engendered by transport across each membrane and also between membranes.
Other hubs of synergy surface for the ALMT anion channel associated with the initial rate of stomatal opening and for the H-ATPase associated with the initial rate of stomatal closure (Fig. 3). These were unexpected, as they are not associated directly with [Ca 2+ ] i regulation. They appear counterintuitive at first, because the ALMT channel is normally active in stomatal closure and the H-ATPase is normally active in stomatal opening. However, they make sense in light of the importance of membrane voltage in determining net solute flux (Thiel et al., 1992;Chen et al., 2012;Jezek and Blatt, 2017). The ALMT channel activates, positivegoing, with a midpoint voltage near 260 mV (Dietrich and Hedrich, 1998;Meyer et al., 2010;Hills et al., 2012) and, therefore, imposes a strong, voltage-dependent brake on membrane hyperpolarization for initial solute uptake. Conversely, H-ATPase activity is key to hyperpolarizing the plasma membrane, and its activity must be suppressed for rapid stomatal closure (Merlot et al., 2007;Minguet-Parramona et al., 2016).
What is more important, therefore, is that this systematic approach offers global insights into the manner in which transport interacts to facilitate ion efflux for stomatal movements. Most useful, then, is to recognize that these emergent interactions arise through the physical requirements for the conservation of mass and charge balance in transport. Just as the oscillations in voltage and [Ca 2+ ] i described previously (Minguet-Parramona et al., 2016) simply reflect a spectrum of frequencies that emerge from the balance of intrinsic transport activities of the guard cell, so too the scope for variation within each parameter defining the various solute transporters arises from the constraints of charge coupling that are intrinsic to all transmembrane activities. In short, each hub and interdependent parameter is the by-product of the combination of characteristics that determine each transport process.

MATERIALS AND METHODS
Local sensitivity analysis is a popular tool for identifying the influence that individual parameters have on the output(s) of a model. The analysis involves repeatedly running a simulation, with a single parameter either increased or decreased, and monitoring outputs. The approach can identify parameters that have a strong influence on model behavior. Coupling between model components (in OnGuard, arising from the fundamental model construction that includes the requirement for charge balance across each membrane) means that changing parameters in isolation does not yield directly any information about the interactions between their parameters imposed by charge balance or other indirect factors such as varied ionic conditions.
To overcome these limitations, we applied Sobol sensitivity analysis (Sobol, 2001), a powerful tool for performing global sensitivity analysis. A detailed description of the procedure is beyond the scope of this article, but in essence, the analysis decomposes the variance in a particular model output into contributions from individual parameters and, crucially for this application, contributions from groups of parameters. Consider a particular model output Y, for example, maximum stomatal aperture, which is a function of a set of parameters, X 1 , X 2 ,.,X N , which define several different transporters. The first-order Sobol index S i for parameter X i is given as and can be interpreted as the proportion of the total variance in Y(V(Y)) that can be attributed to changes in X i . In this expression, V Z denotes the variance (with respect to Z) and E Z (f(Z)) is the expected value of f(Z) with respect to Z. X -i is the set of all parameters excluding X i . The larger this value, the more important X i is considered to be with respect to the output Y. The strength of Sobol analysis is that it can be extended beyond first-order contributions (e.g. S i ) to second and higher orders. Of particular interest here are second-order contributions. In particular, S ij gives the joint contribution of parameters i and j to the variance where S ij gives the proportion of the total variance that can be attributed directly to interactions between these two parameters: in other words, the proportion attributed to i and j minus their individual contributions. Most important, a high value of S ij indicates that large changes in the output can be effected by varying both i and j together: in other words, there is synergy in modifying both parameters that is greater than the sum of their respective effects on the output. Simulations using OnGuard were performed on the model described previously for the Arabidopsis (Arabidopsis thaliana) guard cell (Wang et al., 2012;Blatt et al., 2014;Minguet-Parramona et al., 2016) over the standard diurnal period, but with two changes. First, we incorporated step changes in light intensity instead of ramped intensity variations. Second, we introduced descriptors for light quality, dividing the incident light between red and blue spectral components, and we assigned to blue light control of the various ATPases that previously had been conferred to the total incident light. Using step changes in light intensity simplified the analysis. Incorporating the spectral definitions brings the OnGuard platform into alignment with current knowledge of the actions of light, especially of blue light, on the H + -ATPases (Takemiya et al., 2013;Wang et al., 2014b;Yamauchi et al., 2016) and will allow future refinements in model design, but it had no material effect on model output in our formulation.
For Sobol sensitivity analysis, we stepped the light intensity from 0 to 1,000 mmol m 22 s 21 at the start of the daylight period using a mixture of 900 mmol m 22 s 21 red light and 100 mmol m 22 s 21 blue light, and we stepped the light back to 0 mmol m 22 s 21 at the end of the daylight period. All other environmental conditions were kept constant as described previously (Wang et al., 2012). A subset of 15 of the 22 transporters at the plasma membrane and tonoplast, all contributing to Ca 2+ transport or subject to regulation by [Ca 2+ ] i in the model, was chosen to assess their impact on emergent properties in stomatal control. These transporters and their parameters are identified in an accompanying review (Jezek and Blatt, 2017) and the original descriptions of the OnGuard platform Hills et al., 2012). The transporters were reflected in 45 parameters comprising, for each transporter, the transporter number, the dissociation constant, and the Hill coefficient with reference to the Ca 2+ ligand. The parameters were varied over a range of 620% of their default values. Efficient Monte-Carlo sampling using a Sobol sequence (Sobol, 2001;Tissot and Prieur, 2015) was used to generate parameter combinations that cover evenly the mathematical space associated with the individual parameters. In keeping with standard Sobol analysis, desired parameter ranges were scaled to the unit interval (0-1) for sampling and then rescaled to run through the model. The result was a total of 188,498 parameter combinations to test independently within OnGuard. Simulations were conducted in parallel on an in-house computer cluster. In each instance, we extracted system outputs for seven separate variables: minimum and maximum stomatal aperture, initial rates of stomatal opening and closing, the amplitude of [Ca 2+ ] i oscillations, their median periodicity, and the total number associated with stomatal closure (Minguet-Parramona et al., 2016).
Global sensitivity analysis was performed against the variance in these outputs as well as other first-order outputs derivable from them, such as the dynamic range in stomatal aperture. The analysis was performed using the R function sobolroalhs from the sensitivity package. This function implements the estimation of the Sobol sensitivity indices introduced by Tissot and Prieur (2015) using two Orthogonal Array-Based Latin Hypercubes, which allows the estimation of all closed second-order indices containing the sum of the secondorder effect between two inputs and the individual effects of each input. Such analysis yields indexed values that vary in direct relation to the importance of the sensitivity of the model to the corresponding parameters. In general, Sobol index values should be positive, with a sum equal to 1, and are considered significant when the lower edge of the confidence interval is above 0 and the value itself is above 0.05 (5% risk of a zero parameter value). The code for calculating the Sobol sensitivity indices with OnGuard is available from www. psrg.org.uk.