Compound stress response in stomatal closure: a mathematical model of ABA and ethylene interaction in guard cells

Background Stomata are tiny pores in plant leaves that regulate gas and water exchange between the plant and its environment. Abscisic acid and ethylene are two well-known elicitors of stomatal closure when acting independently. However, when stomata are presented with a combination of both signals, they fail to close. Results Toshed light on this unexplained behaviour, we have collected time course measurements of stomatal aperture and hydrogen peroxide production in Arabidopsis thaliana guard cells treated with abscisic acid, ethylene, and a combination of both. Our experiments show that stomatal closure is linked to sustained high levels of hydrogen peroxide in guard cells. When treated with a combined dose of abscisic acid and ethylene, guard cells exhibit increased antioxidant activity that reduces hydrogen peroxide levels and precludes closure. We construct a simplified model of stomatal closure derived from known biochemical pathways that captures the experimentally observed behaviour. Conclusions Our experiments and modelling results suggest a distinct role for two antioxidant mechanisms during stomatal closure: a slower, delayed response activated by a single stimulus (abscisic acid ‘or’ ethylene) and another more rapid ‘and’ mechanism that is only activated when both stimuli are present. Our model indicates that the presence of this rapid ‘and’ mechanism in the antioxidant response is key to explain the lack of closure under a combined stimulus.


Plant material
Arabidopsis thaliana seeds of Columbia ecotype (Col-0) were sown on Levington F2+S soil and grown under constant conditions in a growth chamber (Sanyo Gallenkamp, Loughborough, UK) with a 10 hour photoperiod a light intensity of 100 − 150µE/m 2 /s, 22 • C and 70% of relative humidity. After 7 days, seedlings were transplanted individually into new pots and were maintained in the same conditions previously detailed. Leaves of 4-6 weeks old plants, which had not yet formed flower bolts, were harvested for the aim of these experiments.

ROS fluorescence assays
Stomatal H 2 O 2 concentration was measured as described in Ref. [8]. Leaves from 4 weeks old plants (approximately 2 leaves per condition) were blended in deionised water and epidermal fragments were collected with a 40µm sterile cell strainer (Fisher Scientific). Epidermal fragments were incubated in Petri dishes containing MES/KCl buffer (5mM KCl, 50 µM CaCl 2 , 10mM MES buffered to pH 6.15 with KOH) inside the growth chamber for 3 hours. This allowed the stomata to open prior to the application of the treatments. Epidermal fragments were collected and equally distributed into Petri dishes loaded with ethanol (as control), 10µM ABA (2-cis, 4-trans abscisic acid 98%, synthetic, Aldrich), 10µM the ethylene precursor ACC (1-aminocyclopropane-carboxylic acid hydro-chloride, Sigma), and a combination of 10µM ABA and 10µM ACC. The treatments, which had durations of 5, 15, 30, and 60 minutes, were performed. Following treatment for the appropriate time, fragments were incubated with 50µM H 2 DCF − DA (2, 7 -dichlorodihydrofluorescein diacetate, Invitrogen) for 10 minutes for H 2 O 2 detection. After a washing step in MES/KCl buffer for 20 minutes, epidermal fragments were placed onto a slide and observed under a microscope. All steps were carried out under dark conditions, as the dye is light-sensitive. H 2 O 2 was visualised with a fluorescence microscope (Axioskop2 plus, Carl Zeiss Ltd., UK) with Zeiss filter set 3 (excitation light filter: 450-490nm, emission light filter: 515-565nm). Images were captured with Axiovision software v3.1 (Carl Zeiss Vision GmbH, UK). Images were processed and fluorescence intensities (as mean of the pixel intensities) were measured with ImageJ software [1].

Stomatal assays
Stomatal aperture bioassays were performed on 4 week old leaves, as described in Ref. [8]. Leaves from 4-6 week-old plants were then cut from the plants using scissors. Excised leaves were floated for 3 hours inside the growth chamber in Petri dishes with buffer (5mM KCl, 50µM CaCl 2 , 10mM MES, buffered to pH 6.15 with KOH). After the initial treatment in buffer the leaves were exposed to 1µM , 10µM ABA, 1µM , 10µM ACC, a combination of 10µM ABA with 10µM ACC, and ethanol as a control. The treatments were left in the growth chamber for 15, 30, 45 and 60 minutes. Two leaves were blended in water for 1-2 minutes and epidermal fragments collected on a 100µm nylon mesh (Spectra-Mesh, BDH-Merck, Nottingham, UK) and transferred to a microscope slide. Measurements of individual stomatal aperture were conducted (see Fig. 1 of the main text) using a Leica DME light microscope, connected to a Leica DFC290 camera imaging system (Leica, Milton Keynes, UK). Leica QWinV3 software (Leica QWIN software, Leica, Milton Keynes, UK) was used to measure the apertures. Data for each data point are the mean of 3 × 30 measurements.

Construction of the stomatal closure model
In this section we describe in detail the steps and assumptions taken to construct the ODE model of stomatal closure under single and combined ABA and ethylene stimuli.

Signal perception and ROS production
In the main text we describe the events in ABA and ethylene perception and production of ROS. These early events are represented in Fig. S1A by the subnetwork comprised of the nodes ABA, Ethylene, PYR/PYL, ABI1, OST1, AtrbohD/F, ROS, AOX, and ETR1, whose concentrations we denote by [ where a subscripted T indicates total concentration of an enzyme (a non-negative constant in R), a subscripted P indicates phosphorylation (or, more generally, being active), and variable names joined with a dash are complexes. For example, the constant [P Y R T ] is the total amount of the ABA receptor PYR/PYL in any form, the variable [P Y R] is the concentration of "available" PYR/PYL molecules, [P Y R-ABA] is the concentration of ABA-bound PYR/PYL molecules, and [P Y R-ABI1] is the concentration of ABI1-bound PYR/PYL. We describe the events known to occur between ABA and ethylene perception and ROS production with the following ODE model: Equations (S1)-(S3) describe the events of ABA perception and ABI1-phosphatase inhibition [12,18], and OST1 phosphorylation [15]. Equations (S4) and (S5) describe the activation of AtrbohD/F where OST1 phosphorylates AtrbohF, unbound ETR1 inactivates AtrbohF, and AtrbohD is activated by ABA [8,11]. Equation (S6) shows ethylene binding and subsequent inactivation of ETR1 [6,7]. Equation (S7) shows the production of ROS by AtrbohD/F and by other Figure S1: Simplification of signal perception and ROS production model. A: Network representation of the ROS production model (S1)-(S8). B: Simplification of the ROS production model given in equations (S9)-(S13). C: Further simplification of ROS production, given in equations (S15)-(S17). cellular processes (summarised in the constant rate k 13 ), scavenging by antioxidants (summarised in the variable AOX), removal or decay of ROS (k 16 [ROS]); equation (S8) shows that AOX has endogenous production (k 18 ), production in response to ROS (k 19 [ROS]), and decay or inactivation (k 20 [AOX]).
There are some things to note about the model above: the way in which the treatments are given (floated on a Petri dish, see Materials and Methods) allows us to assume that the concentration of ABA and ethylene does not change, hence [ABA] and [ACC] are constants. The amount of ATP for the phosphorylation reactions and of NADPH for the production of ROS are also considered abundant enough within guard cells and they are implicitly included in the parameters of the model. We assume that when the complex PYR-ABA-ABI1 dissociates, it does so completely (i.e., into PYR, ABA, and ABI1), to make the equations simpler; we also assume that unbound ETR1 negatively interacts with AtrbohF (or an upstream activator). For now, the variable [AOX] clusters together the group of antioxidants that are active during stomatal closure. The model described in (S1)-(S8) has 8 equations and 20 kinetic parameters plus 6 parameters representing the total amount of the enzymes involved, which makes the task of determining parameter values from the ROS measurements a rather difficult one. We reduce the number of equations and parameters as much as possible by making a series of assumptions, but avoiding oversimplification of the system.
As noted in Refs. [8,11], AtrbohD has a limited role in ABA-induced stomatal closure and none in ethylene-induced closure, therefore we assume that ROS is exclusively produced by AtrbohF. This assumption is supported by our ROS measurements which show similar initial increases in ROS upon an ABA or ethylene stimulus (main text, Fig. 3), for if there were two sources of ROS the pattern of initial increase would most likely differ. We use the quasi-steady-state assumption (QSSA) on the dynamics of the ABA and ethylene receptors, so equations (S1) and (S6) are assumed to reach equilibrium before the other variables have changed significantly. Under the QSSA, the expressions for the ligand-bound receptors become: , .
The expression for [ET R1-ACC] is a standard Michaelis-Menten term with a maximum rate given by the total amount of ETR1. The expression for [P Y R-ABA] is a Michaelis-Menten-type term with a dependency on the variable [P Y R-ABI1] (which holds PYR/PYL molecules). The system in (S1)-(S8) becomes (S13) Figure S1B corresponds to this reduced model where some nodes and edges have been removed but the relationship between the signals and the output is still the same as in Fig. S1A. Though a simplification, this model is still large and we would like to find a further simplification. In the network representation of equations (S9)-(S13) shown in Fig. S1B the path-length from ABA to AtrbohF is 3 while the path-length from ethylene to AtrbohF is 1; however, in our experiments we observe a negligible difference between the ROS produced by ABA and the ROS produced by ethylene five minutes after treatment. We also note that the edge between ethylene and AtrbohF in Fig. S1B is positive because ETR1 inactivates the inactivator of AtrbohF in the absence of ethylene. These observations suggest that the number of steps between ethylene perception and ROS production is about the same as with ABA signals (the immediate events after ethylene binding by ETR1 in guard cells are still unknown), or if the number of events is different then the timescale of the reactions is approximately the same in both cases. Furthermore, in both cases the maximum rate of ROS production is limited by the total amount of AtrbohF ([AtrbohF T ]), which means that the response to either signal has the same theoretical maximum rate. Given that data are unavailable for the receptors, OST1, and ABI1, and that ABA and ethylene have similar timescales for producing ROS we further simplify our model so that AtrbohF becomes active directly from the ABA and ethylene signals, as shown in Fig. S1C.
To determine the equations that represent Fig. S1C, we must have a hypothesis of how ABA and ethylene signals activate AtrbohF. One possibility is that the signals activate AtrbohF through the same pathway, i.e., there is a bottleneck for both signals upstream of AtrbohF. In this case the signals are essentially interchangeable. This assumption would imply that, for example, a 2 µM dose of ABA and a combined 1 µM ABA plus 1 µM ethylene dose are the same:

Figure S2
: Antioxidant production at the end of linear activation cascades. A: Cascade of length n 1 whose input has a logic or function, the cascade responds to ABA or ethylene treatment (equation (S16)). B: Cascade of length n 2 whose input has a logic and function, the cascade becomes active only when ABA and ethylene are present simultaneously (equation (S17)).
The identity of the signalling "bottleneck" remains unknown. Another hypothesis that has a similar result but that does not require a common node in the pathways of the signals would be to assume that the signals converge for the first time at AtrbohF, and activate it independently of each other: where α 1 is a product of [AtrbohF T ] and other rate-limiting parameters in the ABA pathways, and κ 1 is the ABA-specific Michaelis constant, and likewise for ethylene α 2 is the rate-limiting parameter and κ 2 the Michaelis constant. We emphasise that this is not a rigorous derivation of the ROS-activation dynamics but a deduction guided by our current knowledge of the system and assumptions deemed reasonable. In Sec. 4 we derive the compound Michaelis-Menten term in equation (S14) using the QSSA (see Sec. 4). Though both hypotheses of ROS production can produce a similar response the latter one is better suited for modelling stomatal closure because it does not require the assumption of additional interactions, includes the former as a special case (when κ 1 = κ 2 and α 1 = α 2 = [AtrbohF T ]), and its derivation is more clear. Now we turn our attention to the antioxidant pool AOX; it is clear that a homogeneous antioxidant pool responsive only to the concentration of ROS as described by equation (S13) is not compatible with our experimental observations. We consider the possibility of two different antioxidant mechanisms described by the variables [AOX 1 ] and [AOX 2 ] which lie at the end of linear activation cascades driven by [ABA] and [ACC] (Fig. S2). As discussed in the main text, there is evidence to suggest that two distinct antioxidant mechanisms might be at work during stomatal closure. The first of these mechanisms (AOX 1 ) would be activated by either an ABA or an ethylene signal and includes endogenous antioxidant production (to maintain unstimulated equilibrium levels), whereas the second antioxidant (AOX 2 ) is only active when ABA and ethylene signals are present simultaneously.
We place the response of the antioxidant to the signals at the end of cascades of abstract variables to emulate the delay observed in the removal of ROS. Each cascade has a constant input (depending on the ABA and ethylene signals, see below) and has a solution proportional to the normalised lower-incomplete Gamma function P(n i , h i t) (see Sec. 3 where n i represents the length of the cascade and h i is related to the deactivation rates of the cascade, i = 1, 2.
The input of the cascade culminating in AOX 1 must follow a boolean or logic and saturate, therefore we use the compound Michaelis-Menten form we used to describe ROS-production previously as the input for the ABA and ethylene signals, so the equation The cascade leading to the activation of AOX 2 must operate as a logic and gate, ie becoming active only if [ABA] > 0 and [ACC] > 0. We consider that the input of the cascade is downstream of the ABA and ethylene receptors and the response of the cascade must also exhibit saturation, so we model the input as the product of two Michaelis-Menten forms: which enforce a logical and operation of the signals. The expression for AOX 2 is Now we have a complete reduced model of a ROS production-module in guard cells:

NO production
As discussed in the main text, ROS induces NO production in guard cells treated with ABA via the enzyme NIA1 (Fig. S3A); the possibility that NO is also produced in guard cells treated with ethylene was also discussed. An ODE describing endogenous and enzymatic NO production in guard cells dependent on ROS synthesis is where α 30 is a constant rate of NO production by other processes, the Michaelis-Menten term is ROS-induced NO production via NIA1, and the last term is NO decay and removal. Note that in unstimulated guard cells [ROS] > 0 and as measurements of NO to distinguish between ROS and non-ROS induced production are not available, we gather the ROS-dependent and ROS-independent NO production in a single term. We include a second term that describes further enzymatic NO production from ethylene (Fig. S3B), which could be either from NIA1 or another, yet unidentified source (see main text). The new expression for NO production becomes Figure S3: NO production models in guard cells. A: NO-production unit. NO is produced in a ROS-dependent way by NIA1. B: Model of NO-production given in equation (S18). The key assumption of this model is the existence of a ROS-independent pathway of NO production in response to ethylene.

Ca 2+ increase, cytosolic alkalinisation, and ion efflux
The dynamics of Ca 2+ -release and action in guard cells are complex and not yet fully understood [10]. Though the importance of Ca 2+ in guard cell signalling (and cell viability in general) is beyond doubt, in this work we do not include an equation describing its behaviour for three reasons: i. Reports of Ca 2+ -behaviour after ABA treatments in the literature describe both oscillations and rises in cytosolic levels, and experimental data-sets encompassing both single and combined ABA and ethylene treatments do not exist.
ii. As mentioned in Ref. [10], the way in which a cytosolic Ca 2+ rise (or oscillations) transmits signals during stomatal closure is not yet clear. Current hypotheses state that ABA "primes" receptors of Ca 2+ , making the rise in cytosolic levels helpful but not essential for successful closure.
iii. Our experiments suggest that ABA-ethylene cross-talk occurs at the ROS-level, upstream of Ca 2+ in the guard cell signal transduction network (Fig. 2 of the main text), so we direct most of our efforts to understanding signal transduction at this level.
During ABA-induced stomatal closure the pH in the cytosol of guard cells increases from 7.0 to 7.5, and the pH in the tonoplast decreases from 5.5 to 5.0 after treatment with ABA. The concentration of H + determines pH through the Henderson-Hasselbach equation: pH = − log([H + ]) [19]. During stomatal closure the membrane H + -ATPases are inactivated, which means that changes in pH are the result of the transport of protons from the cytosol into the tonoplast by the vacuolar proton pumps (V-ATPases), activated by OST1 [24]. The model of cytosolic alkalinisation and ion-efflux shown in Fig. S4B has one equation for the potassium ion concentration [K + ], and one equation for the outwards potassium channels [K + out ]: Figure S4: Events in stomatal closure downstream of ROS. A: Cytosolic alkalinisation following treatment with ABA, Ca 2+ -increase, membrane depolarisation, and ion efflux. B: Model of the late events in stomatal closure presented in equations (S19) and (S20).
Equation (S19) shows the change in [K + out ], the active outwards K + channels. The first and last terms represent the constant flux of channels between the active and inactive states, respectively. The second term represents the extra number of channels made available by the increase in cytosolic pH following an ABA stimulus, mediated by OST1. The third term (α 42 [N O]) is the increase in K + out activity as a result of membrane depolarisation, possibly via NO-induced Ca 2+ release (ie via the path NO → Ca 2+ ⊣ H + -ATPase → Polarity ⊣ K + out in Fig. S4A). We assume that NO does not target K + out , as ion efflux is required for stomatal closure (we note that although NO has been shown to block K + out in Vicia faba guard cells [23], the authors of the study are unsure whether NO action is specifically targeted to K + out .) Equation (S20) shows the change in [K + ]. The first term represents the increase of ions that enter through the inwards-rectifying channels (K + in ), which are inactivated by NO. The second term is the ion efflux through the outwards channels that is proportional to the active channels [K + out ] and the ion concentration itself. We include an equation for K + out but not for K + in because the alkalinisation of the cytosol has the effect of increasing the number of available channels to extrude ions, whereas the inactivation of K + in is only represented by a term in the equation for K + .

Loss of turgor
The relationship of this model to stomatal aperture is via the last variable [K + ]. Cell volume (and hence stomatal aperture) is determined by the ion and solute concentration in the cell relative to the external concentration [16,25]. Therefore, we take ions and solutes (K + in particular) as a simple proxy for aperture: We assume for this work that the relationship between ionic concentration and aperture is linear; as more data become available (K + and solute concentration) this relationship is likely to change in order to become more realistic, possibly becoming nonlinear.

A model of signal transduction for stomatal closure
We use equations (S15)-(S20) to construct the model of ABA and ethylene-induced stomatal closure represented graphically in Fig. 4 of the main text. We normalise the variables in the model by their non-stimulated equilibrium levels (ie the initial conditions) so they represent percentage of control: With the normalised variables we can transform equation (S21) to so we take the normalised potassium ion concentration as equivalent to the normalised aperture. Dropping the hat notation, the equations of the model become (note that the parameters have been renamed) which is the model we present in equations (1)-(6) of the main text.

Allostery hypothesis for ABA-ethylene cross-talk in ROS production
We have explored an alternative hypothesis to explain the decrease in ROS production in guard cells under a combined ABA-ACC treatment. Motivated by the findings in Ref. [22] which show that the NADPH-oxidase AtrbohF has two phosphorylation sites, we explore the idea that ROS can only be produced when AtrbohF has a single phosphorylation, and that double phosphorylation returns the enzyme to a non ROS-producing state (Fig. S5A). In this scenario, an ABA signal can only lead to the phosphorylation of only one of the two active sites, and likewise an ethylene signal can only lead to the phosphorylation of the other active site only. When the ABA and ethylene signals are present simultaneously time then AtrbohF may become doubly-phosphorylated: once by ABA and once by ethylene, rendering it unable to produce ROS. We use a simple mass-action model to test this hypothesis:  Figure S5B shows the time-courses of the model above to the ROS production along with data obtained from our experiments. The single stimulus solutions of the equations (to ABA or ACC stimuli alone) have a reasonable fit to the data but under a combined stimulus, the equations are not able to reproduce the peak in ROS production five minutes after treatment.

Modelling of activation through linear cascades
We model [AOX 2 ] and the production of [AOX 1 ] as the result of a linear activation cascade with a constant input λ. A model for such cascade is given by the following system of ODEs [9,5]: . . .
with initial conditions x i (0) = 0, ∀i. If we assume that all the deactivation rates of all components in the cascade are identical (i.e., β i = β, ∀i) then the cascade provides optimal amplification [5]  and the exact solution for the last component (x n (t)) is [2]: where P(n, βt) is the normalised lower incomplete Gamma function [17]. where we condense the reactivity, activation and degradation rates in: with the timescale of the cascade β 13 = β.

Compound Michaelis-Menten forms
In this section we derive the compound Michaelis-Menten term we use in our model of ABA and ethylene-induced stomatal closure. Suppose an enzyme E catalyses the production of a substance P from two different substrates S 1 and S 2 in the following chemical reactions: , E T ∈ R + . We eliminate equation (S36) using this conservation relation. Now equations (S37) and (S38) become: The quasi-steady-state approximation (QSSA) states that before any meaningful amounts of P are produced, the enzymes and complexes reach an equilibrium [14,21]. Thus, equations (S40) and (S41) are both equal to zero, because they do not depend on equation (S39) and we can analyse them in isolation: If k a = (k 2 + k 3 )/k 1 and k b = (k 5 + k 6 )/k 4 then The production of P can be expressed as which is a compound Michaelis-Menten form. When [S 1 ] > 0 and [S 2 ] = 0 the expression in equation (S42) becomes proportional to the standard Michaelis-Menten form: likewise when [S 1 ] = 0 and [S 2 ] > 0 we have

Parameter fitting
Equations (1)-(6) of the main text have 28 parameters whose values must be determined. (Note that the variables are rescaled dividing them by 100, so that control levels and initial conditions are 1, to improve numerical stability of the fitting process.) We use our experimental time-course measurements of ROS and aperture to fit the parameters in Eqs.
With these conditions, the number of unknown parameters has been reduced from 28 to 23 and we define the vector θ ∈ R 23 in parameter space: θ = [α 11 , α 12 , k 11 , k 12 , β 11 , β 12 , β 13 , n 2 , α 22 , α 23 , k 21 , k 22 , α 24 , n 1 , β 20 , α 31 , k 31 , α 32 , α 41 , α 42 , α 43 , α 51 , k 51 ] , whose components are all non-negative. We define the set of treatments as T = {T 1 , T 2 , T 3 , T 4 , T 5 } where: where ||·|| 2 is the euclidean norm. That is, we measure the distance between our ROS and aperture measurements and the model for a given θ in the parameter space. The global optimisation problem is to find θ ‡ where We use the Squeeze-and-breathe optimisation method [3] to find θ ‡ . The method requires an initial probability distribution for each parameter (called a prior). In this work, we used a uniform distribution U (0, 10) for all parameters. On each iteration 500 points in the parameter space (in R 23 + ) are sampled from the prior. Each point is used as a starting point to minimise E D (θ) locally (using the Nelder-Mead simplex algorithm). The 50 local minima with the smallest errors are used to construct a posterior distribution of the parameters. The posterior is used as a prior for the next iteration where another 500 points are sampled and minimised until the convergence criteria has been met. Figure S7 shows the convergence of the method for fitting the parameters of our model. On Fig. S7A, we show the decrease in the difference (on a semilogarithmic scale) between the errors of the parameter sets found at the end of each iteration and the global minimum E D (θ ‡ ) ≈ 0.0215,  Figure S7B shows the convergence criterion defined in Ref. [3]. We stop the iterations of the method once the difference between mean of the errors of the 50 parameter sets from consecutive iterations (φ k , shown on a semilogarithmic scale) is smaller that 10 −5 . During the first 20 iterations of the method φ k appears to decrease exponentially and thereafter the trend still continues downwards albeit no longer exponentially. On Fig. S7C we show the mean of the cosines of the angles between all local minima from each iteration. This is to assert that the method converges to a single region of the parameter space. After iteration 43 the mean cosine is O(10 −4 ). Based on these metrics we conclude that θ ‡ is a good estimation of the model parameters, given the present data. See Ref. [3] for a detailed description of the method. Figure S8 shows the distribution of the best 50 parameters after 43 iterations of the algorithm. Red dots mark the mean of each parameter (values in Table 1). The behaviour of the model that we observe in Figs. 3, 5, and 6 of the main text is given by these parameters.

Model selection
We use Akaike's information criterion with a correction for finite sample sizes (AICc) to select the best out of the models presented here. The quantity measured by the AICc balances how well a model fits the data with the number of parameters to discourage overfitting and is a convenient way to perform model selection [4]. The AICc of a model i is given by: where n is the number of observations, p i is the number of parameters in model i, and RSS i is the residual sum of squares of the model. Given many models, the one with the lowest AICc i is considered preferable [4].
We computed the AICc for several versions of the models considered in order to obtain equations (S22)-(S27) and (S28) (which were constructed methodically from known interactions), under different assumptions concerning the value of some of the parameters which reflect different hypotheses about the way in which ABA and ethylene signals are processed in stomata (all values are summarised in Table S2):    Table S2: Number of parameters and values of Akaike's Information Criterion with correction for small sample size (AICc) of the 7 models described in Sec. 6. The number of data points in all the models is n = 36.
• Model 9: 21 parameters, AICc 9 = 268.6. No saturation in AOX 2 activation, and fixed cascade length of 2, no saturation in ethylene-induced production of NO, and no NO effect on K + channels.
• Model 10: 21 parameters, AICc 10 = 142.3. Same as Model 9 except that we assume that the events that activate AOX 1 after ABA or ethylene stimuli have the same speed (equal Michaelis constants), and NO has an effect on K + .
• Model 12: 22 parameters, AICc 12 = 158.6. Length of cascade of AOX 2 is a parameter (n 2 ) , ABA-induced pH rise follows Michaelis-Menten kinetics, and the Michaelis constant in ethylene-induced NO is k 12 , the same as in the ethylene receptors.
Based on the AICc, the best models are 15, 13 and 10. Model 10 is discarded because it exhibits stomatal reopening with increased doses of ABA and ethylene, something that has been disproved experimentally. The error of models 13 and 15 are so close that although AICc 15 < AICc 13 , we consider them equally plausible until further experiments shed light into the behaviour of NO and K + .

Sensitivity analysis
We perform a sensitivity analysis of the model in equations (S22)-(S27) (models 13 and 15 from Sec. 6) using the extended Fourier amplitude sensitivity test (eFAST) [13,26] computed in MAT-LAB with the Systems Biology Toolbox [20]. The eFAST method computes sensitivity scores between 0 and 1 based on how fast the Fourier decomposition frequency of each parameter propagates from variations in inputs [13].  Figure S9: (Colour) Global sensitivities of the parameters in model (S22)-(S27) computed using the extended Fourier amplitude sensitivity test (eFAST). For each parameter, the mean sensitivity is given by the blue line, the bars cover the range between the minimum and the maximum sensitivity observed (1000 simulations, see Refs. [20,13]). The parameter sensitivities in the model are low with the exception of those related to the antioxidant cascades (β 12 , β 13 , α 23 ) which, as described in the main text, are crucial to the observed behaviour of guard cells.
In Fig. S9 we show the sensitivities of the parameters in model (S22)-(S27) with fixed n 1 = 9 and n 2 = 3. The most sensitive parameter in the model is α 23 (mean sensitivity 0.89), the deactivation rate in the cascade that produces AOX 1 in response to either ABA or ethylene. The next two most sensitive parameters are β 12 and β 13 (mean scores 0.63 and 0.68, respectively) which are the scavenging rate of ROS by AOX 2 , and the degradation of the cascade components leading to the activation of AOX 2 . Other sensitive parameters are related to ethylene-induced NO production and its effect on potassium ion extrusion.
In Fig. S10 we examine the sensitivities of the parameters in the model with the alternative hypothesis of pH rise induced by ethylene, discussed in Sec. 2.6. As in the main model, the most sensitive parameters are α 23 and β 13 (mean scores 0.71 and 0.59), related with the timescales of both antioxidants (the sensitivity of β 12 is slightly less than in the previous model). The parameter α 42 in equation (S28) which relates NO to K + extrusion shows increased variability in its sensitivity range. The parameter α 43 , which is involved with the hypothetical effect of ethylene on pH (and thus on the number of K + out channels available) also appears to be sensitive.