Abstract

Motivation: Hormone pathway interactions are crucial in shaping plant development, such as synergism between the auxin and brassinosteroid pathways in cell elongation. Both hormone pathways have been characterized in detail, revealing several feedback loops. The complexity of this network, combined with a shortage of kinetic data, renders its quantitative analysis virtually impossible at present.

Results: As a first step towards overcoming these obstacles, we analyzed the network using a Boolean logic approach to build models of auxin and brassinosteroid signaling, and their interaction. To compare these discrete dynamic models across conditions, we transformed them into qualitative continuous systems, which predict network component states more accurately and can accommodate kinetic data as they become available. To this end, we developed an extension for the SQUAD software, allowing semi-quantitative analysis of network states. Contrasting the developmental output depending on cell type-specific modulators enabled us to identify a most parsimonious model, which explains initially paradoxical mutant phenotypes and revealed a novel physiological feature.

Availability: The package SQUADD is freely available via the Bioconductor repository at http://www.bioconductor.org/help/bioc-views/release/bioc/html/SQUADD.html.

Contact:martial.sankar@unil.ch; christian.hardtke@unil.ch

Supplementary information:Supplementary data are available at Bioinformatics online.

1 INTRODUCTION

The ontogenesis of plants is characterized by an intrinsic developmental plasticity, which reflects their capacity to adapt to environmental conditions and is frequently conveyed by modulation of plant hormone pathways (Nemhauser et al., 2006). Thus, plant hormones are not only essential for various endogenous developmental programs, but also for the perception of, and adaptation to, environmental change. Several plant hormone pathways are known, and their often context-specific, overlapping activities suggest that they influence each other, sometimes in a hierarchical, sometimes in a synergistic manner. Auxin and brassinosteroids are among the best characterized hormones with well-defined signaling pathways and cellular effects. Both hormones are limiting for cell elongation, a process in which they act synergistically because of crosstalk between their pathways (Hardtke, 2007; Kuppusamy et al., 2008).

A distinct feature that sets auxin apart from other plant hormones is the fact that it can be moved through the plant body by specialized molecular machinery. This regulated and directional transport, called polar auxin transport (PAT), contributes significantly to the local control of auxin action (Leyser, 2005). Although other hormones, including brassinosteroids (Savaldi-Goldstein et al., 2007), might well be mobile within a certain range from their site of production, specialized transport systems apart from PAT have not been found so far. Together with auxin biosynthesis, PAT determines cellular auxin concentration, which is translated into a transcriptional response via the canonical auxin signaling pathway. This involves nuclear auxin receptors, transport inhibitor response 1 (TIR1) and homologs, a class of SCF-type E3 ubiquitin ligase F-box proteins that are activated by binding auxin, hence targeting transcriptional co-repressors of the auxin/indole-3-acetic acid (AUX/IAA) family for proteasome-mediated degradation (Dharmasiri et al., 2005; Kepinski and Leyser, 2005). Since AUX/IAAs inhibit the activation potential of auxin response factors (ARFs), this releases ARFs to activate transcriptional targets of auxin signaling through auxin-responsive promoter elements (AuxREs) (Benjamins and Scheres, 2008).

Importantly, auxin signaling involves a negative feedback loop, as AUX/IAA genes are among the most prominent, early auxin signaling targets. Moreover, auxin signaling also affects PAT, because it controls the expression of several of the integral plasma membrane auxin efflux carriers, the PIN-FORMED (PIN) proteins (Sauer et al., 2006; Vieten et al., 2005; Wisniewska et al., 2006). Therefore, auxin transport and signaling are intimately linked, and local auxin activity is conveyed by their interplay (Benjamins and Scheres, 2008; Leyser, 2005).

The brassinosteroid pathway represents a more classic signaling paradigm, where perception of the hormone at the plasma membrane modulates the activity of nuclear targets to eventually trigger gene expression changes (Kim et al., 2009; Vert et al., 2005). Brassinosteroids are perceived by integral plasma membrane proteins with a leucine-rich repeat receptor-like kinase topology, the brassinosteroid insensitive 1 (BRI1) and BRI1-like proteins. Binding of brassinosteroid promotes association of the receptors with another membrane-bound kinase, BRI1-associated receptor kinase 1, triggering a series of phosphorylation events. This signaling cascade, via several intermediates, inhibits the activity of BRASSINOSTEROID-INSENSITIVE 2 (BIN2), a glycogen synthase kinase 3/shaggy-like kinase. In the absence of brassinosteroid, BIN2 phosphorylates the highly homologous transcription factors bri1-EMS-SUPPRESSOR 1 (BES1) and BRASSINAZOLE RESISTANT 1 (BZR1). Thus, upon brassinosteroid signaling, aided by the phosphatase bri1-SUPPRESSOR 1 (BSU1), BES1 and BZR1 accumulate in a dephosphorylated state, which results in altered nucleo-cytoplasmic partitioning and DNA binding affinity towards brassinosteroid-responsive promoter elements (BRREs), and thus target gene activation. This signaling pathway involves a homeostatic feedback loop, since it directly impinges on brassinosteroid biosynthesis.

The pathways described above appear to be universally expressed throughout the plant, and have been implicated in many developmental or physiological processes. This is somewhat paradoxical, given that many of these processes do not appear to be related. However, while less variability has been observed in the brassinosteroid signaling pathway, (tissue-)specificity of auxin effects could be explained by the complexity of the auxin signaling components, notably the numerous AUX/IAA and ARF genes. The respective proteins are likely not fully redundant, thus their differential expression could introduce a cell type-specificity to auxin response (Weijers et al., 2005). Moreover, the ARF proteins do not only encompass transcriptional activators, but also transcriptional repressors. Another explanation could be that auxin, and brassinosteroid, action depends in part on molecular pre-patterns that could involve context-specific modulators and read-outs of auxin action (Badescu and Napier, 2006; Benjamins and Scheres, 2008; Braun et al., 2008; Strader et al., 2008; Tromas et al., 2009). One such modulator could be the BREVIS RADIX (BRX) gene, which is specifically expressed in the vasculature and is rate-limiting for auxin action, likely by impinging on brassinosteroid biosynthesis (Beuchat et al., 2010; Mouchel et al., 2006). As BRX activity is controlled by auxin at both the transcriptional and post-translational level (Scacchi et al., 2009), this suggests that BRX mediates crosstalk between the two hormone pathways, adding to accumulating evidence for a rate-limiting role of the brassinosteroid pathway in auxin response (Hardtke, 2007; Kuppusamy et al., 2008; Nemhauser et al., 2004; Vert et al., 2008). Another factor involved in auxin–brassinosteroid crosstalk is the repressive auxin response factor, ARF2, which has been shown to be a substrate of the BIN2 kinase (Vert et al., 2008) and thus a direct point of crosstalk between the two hormones.

Despite the detailed knowledge about the composition of the signaling pathways, their discrete output in a given context still remains largely obscure due to the quantitative nature of hormone signaling and the technical difficulty in measuring it. The best existing read-outs so far are (artificial) reporter genes, which allow some generic quantification of hormone pathway activity by response to endogenous transcription factors (Muller and Sheen, 2008; Sabatini et al., 1999). However, because of the multiple feedback loops in each pathway and a lack of quantitative data on the activity of individual components, e.g. the level and stability of transcripts and proteins, steady state signaling amplitude in a given condition and its developmental consequences are still nearly impossible to predict.

In an attempt to overcome these limitations, we sought to apply network modeling, which has proven useful in deciphering, sometimes paradoxical, experimental observations in Arabidopsis (e.g. Digiuni et al., 2008; Liu et al., 2010; Locke et al., 2006). We applied a Boolean logic approach (e.g. Pandey et al., 2010) to assess and analyze the auxin and brassinosteroid signaling networks, and their interaction. To compare the behavior of the network components across one or several conditions, we transformed a discrete dynamic model into a qualitative continuous system using the SQUAD software (Di Cara et al., 2007). This approach has been successfully applied previously to model gene regulatory networks (Philippi et al., 2009; Sanchez-Corrales et al., 2010). A comparison to alternative strategies …(e.g., Diaz and Alvarez-Buylla, 2009) can be found in (Wittmann et al., 2009). Our model was able to explain various paradoxical experimental observations, thereby helping us to correctly place the BRX gene in the respective network and approximate how the cell type-specific expression of BRX could contribute to differential cell fate.

2 METHODS

2.1 Plant growth, genotyping and hormone measurements

The Arabidopsis brx-2 and arf2-8 null mutants used to create the double mutant have been described (Scacchi et al., 2009; Vert et al., 2008). Plant tissue culture, genotyping and auxin measurements were performed according to standard methods as described (Beuchat, et al., 2010; Sibout, et al., 2006).

2.2 Signaling network representation

To reconstruct the auxin and brassinosteroid pathways, literature information was used to assemble the signaling circuits with respect to the logical formalism and signaling network framework. As described (Klamt et al., 2006), signaling network models are structured by input, intermediate and output layers. The input represents the starting points of a signaling circuit, which formally are nodes without incoming arrows. Inversely, output nodes are circuit end points that can depict developmental or physiological output processes (Hyduke and Palsson, 2010). Together, the input and output layers define the boundaries for the intermediate layer, which represents the core signal transduction cascade. Due to their linear nature, signaling networks are prone to node reduction, which permits to decrease the complexity of the discrete analysis and avoid any delay issue during the analysis of the continuous form (Naldi et al., 2010). The network models were implemented using cellDesigner v4.0.1 (www.celldesigner.org). The model descriptions are available in XML format in the Supplementary Material.

2.3 Logical model simulation

Model simulations were performed using SQUAD v2.0 (Di Cara et al., 2007), which relies on a standardized qualitative dynamical systems method (Mendoza and Xenarios, 2006) to provide an ordinary differential equation (ODE) to each model component (Equation 1). This leads to the transformation of the original discrete step function into a sigmoid response curve. Each node's ODE relies on two parameters: the gain of the sigmoid, h, and the decay gi. In the absence of kinetic data, in this study we used the default values of h = 10 and gi = 1.

2.4 SQUAD add-on

An add-on to SQUAD developed for this study is provided as an R package named SQUADD (SQUad ADD-on) in the Supplementary Material. This add-on permits the generation of simulation matrices and prediction heat maps. Simulation matrices are useful to assess the individual node simulation profiles across several conditions (e.g. Supplementary Figure S1), whereas prediction heat maps are useful to analyze the node activation change between, e.g. a perturbed condition and a ‘wild-type’ condition (e.g. Fig. 5). The method takes the SQUAD simulation datasets as an input and interpolates the data points with a locally weighted smoothing line or least square fitting line. At a user-defined time value, the ratio of the interpolate values between two conditions can then be calculated. Finally, the color scale is defined according to the ratio range. In this study, we applied a lowess interpolation at t = 50. The SQUADD software including a tutorial can be downloaded as a Bioconductor package at http://www.bioconductor.org/help/bioc-views/release/bioc/html/SQUADD.html.

3 RESULTS

3.1 Modeling the core of cellular auxin perception

The goal of our study was to develop a cellular model for auxin–brassinosteroid interaction in cell elongation. Roots grow at their tips, where an apical meristem that harbors stem cells continuously produces new cells that undergo a stage of proliferation, followed by elongation and eventually differentiation (Osmont et al., 2007). Both auxin and brassinosteroids are limiting in this process. To build up our core models of the auxin and brassinosteroid pathways, we relied on solid, well-established experimental data from the literature as summarized in recent reviews (Benjamins and Scheres, 2008; Hardtke, 2007; Kieffer et al., 2010; Kim and Wang, 2010; Scheres and Xu, 2006; Vert et al., 2005). An overview of these regulatory interactions can be found in Supplementary Table S1. We started by creating a model of the core auxin signaling loops, which includes the AUX/IAA transcriptional repressors, the TIR1-type auxin receptors and the ARF transcriptional activators. As a generic developmental read-out of the model we defined auxin-responsive elements and thus genes that are controlled by ARF activity, including AUX/IAA genes as well as PIN genes. Finally, auxin itself was added as a continuous stimulus in the model, based on the observation that shoot-derived auxin is transported towards the root tip from cell to cell via PAT during seedling development. Evacuation of auxin from the cell by PIN proteins was integrated into the model as a negative effect of PINs on cellular auxin activity.

3.1.1 Modeling strategy

Due to the absence of quantitative data for most of the signaling components and the inherent difficulty to obtain such, ideally kinetic data at spatio-temporal resolution in a developmental context, we analyzed the signaling network by a Boolean formalism approach (i.e. a logical model). In a Boolean formalism, each node (i.e. model component) is defined by a variable x, which represents its state of activation (x = 1) or inactivation (x = 0). The state of each node depends on its interaction or regulatory relationship with other nodes. By extension, the activation state of the system (i.e. of the totality of nodes) is given by the vector of size n (n = number of nodes) of the node activation states. Starting from given initial conditions and in response to a stimulus (in our models auxin and/or brassinosteroids), this vector is then updated synchronously or asynchronously, until the system reaches a stable, or oscillating (complex) steady state.

Boolean formalism is powerful but not very intuitive as the range of values is discrete and kinetics are not considered, i.e. the time scale is arbitrary. Components can be updated synchronously or asynchronously, however synchronous update of the vector is most of the time unrealistic, because it would assume that the components are activated or inhibited simultaneously. This is not the case in biological systems, where components are activated or inhibited on different time scales that depend on the kinetics of individual interactions (e.g. affinity in protein–protein interactions) or the speed of the activation/inhibition process (e.g. transcription versus translation or protein degradation). To overcome the lack of time scales in the Boolean approach, computational tools have been developed, such as the SQUAD software, in which vector updates are entirely asynchronous (Di Cara et al., 2007). For example, in a given sequence of A activates B activates C, a change of state in A only triggers a change of state in B with a certain delay. A change of state in C will then again only follow with a delay after the state of change in B. Moreover, SQUAD permits assessment of the behavior of components across an arbitrary time scale. This is because SQUAD transforms a classic dynamic Boolean model into a continuous one by using ordinary differential equations, one for each node, given by the following equation:  

formula
With this equation, intermediate values between 0 and 1 and therefore sigmoid curves can be obtained, which are useful to assess the behavior of components. The h parameter accounts for the sigmoid gain, i.e. it controls the steepness of the curve (high values will give curves close to a step function). gi represents the decay of the node. wi is the continuous form of the truth table of each node, i.e. the combinatorial effect of all the inputs into a given node. Due to lack of kinetic information, the computations in our models were performed with the default values of h = 10 and gi = 1. It is important to note however, that kinetic data that should eventually become available could be easily integrated into our models by setting h and gi parameters accordingly, allowing model refining. The logic rules and equations describing the nodes of all models presented in this study can be found in the Supplementary Material.

3.1.2 Behavior of the core SCFTir1-AUX/IAA-ARF auxin perception module

Based on the experimental data and strategy outlined above, we created a model for the core auxin perception machinery (Fig. 1A), which can be handled like a circuit. A signal is injected into the system via an input node and its route can be followed through the OR/AND gates and the node state switches until it is caught by the terminal output node. The biological signal for our model is auxin, which was thus set as the initial stimulus. This auxin stimulation node represents the events of both passive and active auxin influx. As the ultimate output of the model, we defined a generic ‘elongation growth’ node, which represents the developmental output of altered auxin-responsive gene expression on cell elongation. Finally, the tendency of the activation pattern of each node can be represented by a linear regression, which permits a numerical estimation of node behavior (see below).

Fig. 1.

The auxin and brassinosteroid signaling network models. (A) Network representation of the auxin signaling core. (B) Simulation output of the network presented in (A) after transformation into a qualitative continuous system using SQUADD (initial continuous auxin stimulation xIAA = 1, initial node states x = 0). (C) Network representation of the brassinosteroid signaling core. (D) Simulation output of the network presented in (C) after transformation into a qualitative continuous system using SQUADD (initial node states x = 0).

Fig. 1.

The auxin and brassinosteroid signaling network models. (A) Network representation of the auxin signaling core. (B) Simulation output of the network presented in (A) after transformation into a qualitative continuous system using SQUADD (initial continuous auxin stimulation xIAA = 1, initial node states x = 0). (C) Network representation of the brassinosteroid signaling core. (D) Simulation output of the network presented in (C) after transformation into a qualitative continuous system using SQUADD (initial node states x = 0).

The SQUAD analysis indicates that the network components exhibit oscillatory behavior (Fig. 1B), including the elongation output, which is revealed by the sigmoid feature of the program. This result is consistent with the idea of inherent buffering in signaling networks, which prevents a signaling lock-in as well as exaggerated response to transient stimulus and thus provides developmental flexibility. As SCFTIR1 is activated by auxin, Aux/IAAs are inhibited, allowing ARFs to reach their maximum activity, while decreasing SCFTIR1 activity stabilizes Aux/IAAs and thus inhibits ARFs. Importantly, to avoid a lock-in steady state of the network, negative feedback on auxin through PAT is required, which is mediated by the positive effect of auxin perception on PIN activity. In summary, the model correctly simulates the behavior of the core components in auxin perception. Due to the solid experimental evidence for the functioning of the auxin perception machinery, we considered the behavior of its components as a control in the analyses of subsequent, more complex models.

3.2 Modeling of brassinosteroid signaling and its interaction with auxin signaling

In the next step, we sought to extend our model by integrating the brassinosteroid signaling pathway and thus account for the observed synergism between the two hormone pathways (Hardtke, 2007; Kuppusamy et al., 2008; Mouchel et al., 2006; Nemhauser et al., 2004). First, similar to the auxin model, we built a model of brassinosteroid signaling based on experimental data from the literature (Kim et al., 2009; Kim and Wang, 2010; Vert et al., 2005) (Fig. 1C). Beyond the signal transduction from BRI1 to the BES1/BZR1 transcription factors and the eventual modulation of brassinosteroid-responsive gene expression, we also included an explicit feedback on brassinosteroid biosynthesis. This feedback is not only an important component of brassinosteroid homeostasis, but also the major dampening factor of the signaling pathway, conceptually similar to AUX/IAA feedback on auxin signaling. In response to a brassinosteroid stimulus, this model again swiftly reaches a quasi-steady state, including an oscillating developmental output (Fig. 1D). Thus, both of our models appear to correctly capture the signaling features described in the literature.

3.2.1 Auxin–brassinosteroid interaction—the ARF2 connection

To integrate the two models, auxin- and brassinosteroid-responsive gene expression was linked to a common developmental output node, reflecting the observed synergism of the two hormones in controlling cell elongation (Hardtke, 2007; Nemhauser et al., 2004; Wang et al., 2005). Thus, in the integrated model, the developmental output reports the impact of both auxin- and brassinosteroid-responsive promoter elements, which are indeed frequently found together in genes that respond to both hormones.

Beyond the common occurrence of auxin- and brassinosteroid-responsive elements in promoters, only few components that might mediate more direct, non-genomic auxin–brassinosteroid interaction have been described, such as the inhibitory ARF2. ARF2 activity is directly modulated through the BIN2 kinase, which can phosphorylate ARF2 to prevent it from DNA binding (Vert et al., 2008). Thereby brassinosteroids would activate an inhibitory component of auxin signaling, since all ARFs can bind to the same auxin-responsive promoter elements (Guilfoyle and Hagen, 2007). Moreover, a significant effect of ARF2 on genes involved in auxin biosynthesis has been reported, suggesting that brassinosteroid might modulate auxin biosynthesis (Vert et al., 2008). However, the contribution of cellular auxin biosynthesis to our model is negligible for the developmental output as long as auxin flow through the cell is maintained by PAT, which is consistent with established models and experimental findings (Grieneisen et al., 2007). Based on these reported results, we placed ARF2 as a central node connecting auxin and brassinosteroid signaling in our model (Fig. 2A). Moreover, we also integrated ARF2 into the auxin signaling pathway, because inhibitory ARFs, just like activating ones, do interact with AUX/IAA proteins, although at seemingly lower affinity (Tiwari et al., 2003).

Fig. 2.

Interaction between auxin and brassinosteroid signaling through ARF2. (A) Network representation of the integrated pathways connecting auxin and brassinosteroid signaling through the crosstalk component ARF2. (B) Simulation output of the network presented in (A) after transformation into a qualitative continuous system using SQUAD (continuous auxin stimulation xIAA = 0, initial node states x = 0).

Fig. 2.

Interaction between auxin and brassinosteroid signaling through ARF2. (A) Network representation of the integrated pathways connecting auxin and brassinosteroid signaling through the crosstalk component ARF2. (B) Simulation output of the network presented in (A) after transformation into a qualitative continuous system using SQUAD (continuous auxin stimulation xIAA = 0, initial node states x = 0).

The exact role of inhibitory ARFs (ARFi) is not clear, and it has been suggested that they could act by competing with activating ARFs (ARFa) for binding sites, or by directly inhibiting them upon dimerization (Guilfoyle and Hagen, 2007). For our logical model, this is directly not relevant as either mode of interaction impinges on kinetic parameters, but not the logic. In our model, we placed ARF2 under AUX/IAA control, such that AUX/IAA degradation in response to auxin stimulus would remove one level of inhibition from ARFi. Notably, however, in our logical model the exact mode of ARF2 inhibition does not matter for the working of the auxin core perception module, i.e. the quasi-steady state behavior is similar, although the (unknown) signaling amplitude might differ with or without ARFi. For the integrated model, again an oscillating quasi-steady state developmental output was observed, however with a slightly more restricted amplitude (Fig. 2B).

3.2.2 Auxin–brassinosteroid interaction—integration of BRX

A subject of special interest to us was to integrate the BRX gene into our model, since its position with respect to auxin and brassinosteroid signaling remains obscure despite the experimental evidence collected so far (see above). In particular, we aimed to evaluate a suspected role of BRX in an accessory auxin signaling pathway, which could serve to modulate auxin-response in a context-dependent manner (Scacchi et al., 2009).

To place BRX in our model, we first considered the regulatory influence of auxin on BRX activity. The BRX promoter contains prototypical AuxREs, and consistently, BRX expression is highly auxin-inducible, presumably through the canonical auxin signaling core (Mouchel et al., 2006; Scacchi et al., 2009). At the same time, auxin has a more direct effect on BRX protein, by promoting its transfer into the nucleus, where it presumably regulates transcription in conjunction with NGATHA (NGA) B3-type transcription factors (Scacchi et al., 2009). Thus, we placed BRX directly under auxin control. This regulation was considered as activating, although the nuclear transfer of BRX eventually results in BRX degradation (Scacchi et al., 2009). The latter would only restrict the duration of nuclear BRX activity, however, which can be neglected in a logical model.

The positive regulation of BRX by auxin, both at the transcriptional and post-translational level, would not influence the developmental output unless BRX activity feeds back into the network. This feedback is assumed to pass through the brassinosteroid pathway, possibly by promoting brassinosteroid biosynthesis, since brx phenotypes can be largely rescued by brassinosteroid application (Mouchel et al., 2006), and since gain-of-function lines constitutively and ectopically over-expressing full length BRX or a stabilized, partially active C-terminal fragment contained significantly higher levels of the major active brassinosteroids, brassinolide and its precursor, castasterone (Beuchat et al., 2010). Therefore, a positive influence of BRX on brassinosteroid biosynthesis was included in our model (Fig. 4A).

Finally, we also included the interaction of activated (nuclear localized) BRX with NGA transcription factors into our model, since the latter have been proposed to positively regulate auxin biosynthesis (Trigueros, et al., 2009). Whether BRX has an inhibitory or an activating influence on NGA remains to be experimentally determined. However, for our model this was not relevant, because as stated above, as long as auxin is provided by PAT, cellular auxin biosynthesis is negligible for the developmental output.

Different variants of the full model (i.e. including ARF2, BRX and NGA1) were created (Fig. 4A; Table 1) to account for the unknowns described above, such as the regulatory relation between BRX and NGA1. Moreover, we also created variants that assumed a potential direct influence of BRX on auxin-responsive gene expression (Scacchi et al., 2009). Only variants with a valid (i.e. steady state) developmental output were retained for further analysis (variants 4.1.1, 4.2.1, 5.1.1 and 5.2.1; Table 1) (Supplementary Fig. 1). Developmental output that could not reach a steady state any longer and fluctuated in a stochastic manner was observed if the positive link between BRX and brassinosteroid biosynthesis was removed, reinforcing its significance.

Table 1.

Hierarchical summary view of the different models generated and tested in this study

Model version Description No. of nodes No. of edges Feedbacks ±(−) Developmental output valid? 
Auxin core model 2 (1) Yes 
BR core model 1 (0) – 
Integrated auxin-BR model 18 24 4 (1) – 
 Add BRX 
IAA activates BRX 20 28 5 (1) – 
4.1 BRX activates NGA1 20 29 6 (1) – 
4.1.1 BRX activates BR biosynthesis 20 29 6 (1) Yes 
4.1.2 BRX inhibits BR biosynthesis 20 29 6 (1) No 
4.2 BRX inhibits NGA1 20 29 6 (1) – 
4.2.1 BRX activates BR biosynthesis 20 29 6 (1) Yes 
IAA activates BRX activates ARF(a) 20 29 5 (1) – 
5.1 BRX activates NGA1 20 30 6 (1) – 
5.1.1 BRX activates BR biosynthesis 20 30 6 (1) Yes 
5.1.2 BRX inhibits BR biosynthesis 20 30 6 (1) No 
5.2 BRX inhibits NGA1 20 30 6 (1) – 
5.2.1 BRX activates BR biosynthesis 20 30 6 (1) Yes 
BRX activates AuxRE 20 30 6 (1) No 
BRX inhibits Aux/IAA 20 30 6 (1) No 
Model version Description No. of nodes No. of edges Feedbacks ±(−) Developmental output valid? 
Auxin core model 2 (1) Yes 
BR core model 1 (0) – 
Integrated auxin-BR model 18 24 4 (1) – 
 Add BRX 
IAA activates BRX 20 28 5 (1) – 
4.1 BRX activates NGA1 20 29 6 (1) – 
4.1.1 BRX activates BR biosynthesis 20 29 6 (1) Yes 
4.1.2 BRX inhibits BR biosynthesis 20 29 6 (1) No 
4.2 BRX inhibits NGA1 20 29 6 (1) – 
4.2.1 BRX activates BR biosynthesis 20 29 6 (1) Yes 
IAA activates BRX activates ARF(a) 20 29 5 (1) – 
5.1 BRX activates NGA1 20 30 6 (1) – 
5.1.1 BRX activates BR biosynthesis 20 30 6 (1) Yes 
5.1.2 BRX inhibits BR biosynthesis 20 30 6 (1) No 
5.2 BRX inhibits NGA1 20 30 6 (1) – 
5.2.1 BRX activates BR biosynthesis 20 30 6 (1) Yes 
BRX activates AuxRE 20 30 6 (1) No 
BRX inhibits Aux/IAA 20 30 6 (1) No 

3.3 Behavior and testing of the integrated models

Similar to the previous model, the extended versions produced an oscillating quasi-steady state developmental output in most variants (Fig. 4B). This output displayed different trajectories however, consistent with the idea that accessory factors can introduce differential cell fate (Fig. 4B). We would like to note that the cyclic attractors (i.e. oscillating steady states) are consistent with various experimentally observed feedbacks, for instance the transcriptional re-activation of AUX/IAAs in response to their auxin-mediated degradation and the associated oscillation of auxin response. Whether other component states predicted to oscillate by our models in fact do so, for instance hormone levels, is not known and experimentally challenging to determine. Dynamic reporters that can be followed by live imaging at high temporal and cellular resolution would be required for this purpose. We would like to highlight however that modified logical rules can be found to eliminate cyclic attractors in our model if justified, and that our model describes an abstracted signaling network rather than, for instance, a gene regulatory network. Even if the steady states emerging from our models are oscillatory, they clearly describe the difference between the perturbated conditions.

In the next step, we aimed to test whether the models accurately reflect observed experimental findings. First, we perturbed the models by removing BRX a short time after stimulation (Fig. 4C). This resulted in dramatically altered developmental output, which even reached a linear steady state close to inactivity as compared to the full model in two variants (5.1.1 and 5.2.1). By contrast, removing ARF2 from the model did not have any major effect on the developmental output. Both results match with the severity of the root growth phenotypes of the respective mutants, i.e. the strong short root phenotype of brx versus absence of a morphological root phenotype in arf2 (Mouchel et al., 2004; Vert et al., 2008).

Finally, we tested whether the developmental output of our model could be rescued by a brassinosteroid stimulus after a loss of BRX function, simulating the phenotypic rescue of brx mutants by brassinosteroid application (Mouchel et al., 2006). This was indeed the case, i.e. adding brassinosteroid stimulus to the model after elimination of BRX activity and subsequent breakdown of developmental output swiftly restored the output to a pattern that was largely similar to the unperturbed model (Fig. 4D). Interestingly, this was only possible if BRX was not positioned within the brassinosteroid signal transduction module.

3.3.1 Testing the model - genetic interaction between ARF2 and BRX

Experimental support for the notion that BRX does not fit into the brassinosteroid signaling chain was obtained from the analysis of arf2 brx double mutants. Previously, it was found that the long hypocotyl 5 (hy5) mutation significantly suppresses the brx root growth phenotype in brx hy5 double mutants (Scacchi et al., 2009). Supposedly, this is due to constitutively increased auxin-responsive transcription as conferred by hy5 loss of function (Sibout et al., 2006), which can in turn offset diminished auxin-responsiveness in brx mutants to some degree. Along these lines, it was speculated that quantitatively impaired auxin-response in brx mutants might be compensated by arf2 mutation. To test this hypothesis, we generated an arf2 brx double mutant from respective null alleles. Assessment of the root growth phenotype of this double mutant revealed that it was not significantly different from the brx control line (Fig. 3A). Thus, the arf2 mutation cannot suppress the effect of brx, which is consistent with perturbation of our models, i.e. the developmental output after removal of BRX and ARF2 activity from the model was identical to the output observed after BRX removal only (Fig. 4E). This finding, together with the brassinosteroid rescue of brx phenotypes, would also suggest that the main point of auxin–brassinosteroid synergism at the transcriptional level are genes that are activated by both AuxREs and BRREs in their promoters.

Fig. 3.

Auxin content in brx mutants and root growth phenotypes. (A) Primary root growth in wild-type, brx mutant, arf2 mutant and brx arf2 double mutant seedlings in tissue culture. (B) Free auxin (IAA), auxin conjugate (IAA-Asp, IAA-Glu) and auxin precursor (IAM) content of brx mutant seedlings (brxS), their wild-type background (Sav-0) and gain-of-function lines obtained from ectopic over-expression of the BRX C-terminus (BRXCT) or full length BRX.

Fig. 3.

Auxin content in brx mutants and root growth phenotypes. (A) Primary root growth in wild-type, brx mutant, arf2 mutant and brx arf2 double mutant seedlings in tissue culture. (B) Free auxin (IAA), auxin conjugate (IAA-Asp, IAA-Glu) and auxin precursor (IAM) content of brx mutant seedlings (brxS), their wild-type background (Sav-0) and gain-of-function lines obtained from ectopic over-expression of the BRX C-terminus (BRXCT) or full length BRX.

Fig. 4.

Alternative integration of BRX into the network model. (A) Extended model of auxin–brassinosteroid interaction after integration of BRX and NGA1. Model variants that result in valid (i.e. coherent) developmental output throughout perturbed conditions are indicated. (B) Developmental output response curve of the different model variants over an arbitrary time scale (black), with the least square regression (red) indicated (continuous auxin stimulus x = 1). (C) Developmental output response curve after inactivation of BRX (xBRX coerced to 0 at t = 13). (D) Developmental output response curve after inactivation of BRX (xBRX coerced to 0 at t = 13) and rescue by brassinosteroid stimulus (continuous BR stimulus xBR = 0.4 at t = 25). (E) Developmental output response curve after inactivation of BRX (xBRX coerced to 0 at t = 13), followed by inactivation of ARF2 (xARFi coerced to 0 at t = 25).

Fig. 4.

Alternative integration of BRX into the network model. (A) Extended model of auxin–brassinosteroid interaction after integration of BRX and NGA1. Model variants that result in valid (i.e. coherent) developmental output throughout perturbed conditions are indicated. (B) Developmental output response curve of the different model variants over an arbitrary time scale (black), with the least square regression (red) indicated (continuous auxin stimulus x = 1). (C) Developmental output response curve after inactivation of BRX (xBRX coerced to 0 at t = 13). (D) Developmental output response curve after inactivation of BRX (xBRX coerced to 0 at t = 13) and rescue by brassinosteroid stimulus (continuous BR stimulus xBR = 0.4 at t = 25). (E) Developmental output response curve after inactivation of BRX (xBRX coerced to 0 at t = 13), followed by inactivation of ARF2 (xARFi coerced to 0 at t = 25).

3.3.2 Semi-quantitative presentation of perturbation consequences by predictive heat maps

In the next step, we aimed to determine which of the integrated models best accommodate additional experimental data. To facilitate this task, we generated predictive heat maps to assess the activity state change of components between a perturbation and the normal condition. As SQUAD uses an exponential interpolation of Boolean states to obtain continuous activation values, and since activation states do not really represent quantitative experimental data, an absolute quantitative prediction of component activities is not possible. However, a relative, semi-quantitative analysis can be performed to reflect a qualitative change in behavior of components between any two conditions. For such an evaluation, we implemented an R-based script to generate predictive heat maps using a simple algorithm (see Supplementary Material, SQUAD add-on R package). This tool takes the SQUAD simulation datasets as an input and interpolates the various data points for components with a locally weighted smoothing line at a user-defined time value. To compare two model variants, it then calculates the ratio of the interpolated values for a given component between a perturbation and the normal condition. Integrated across the time scale, this then indicates a relative change in the activity status of a node, i.e. an overall increase or decrease in activity. Finally, this relative change is then presented in a color scale grid, permitting intuitive assessment of the consequences of model perturbations and variations (Fig. 5).

Fig. 5.

Semi-quantitative evaluation of model output by predictive heat maps. (AD) Predictive heat maps of relative changes of node activities in the model variants and perturbations as described in Figure 4, in comparison to the respective unperturbed models.

Fig. 5.

Semi-quantitative evaluation of model output by predictive heat maps. (AD) Predictive heat maps of relative changes of node activities in the model variants and perturbations as described in Figure 4, in comparison to the respective unperturbed models.

Starting from our set of model variants, we generated predictive heat maps to compare the relative node states for the perturbations described above. We then determined for which variant model the experimental observations in the brx mutant fit best, focusing on two key results: (i) the observation that the expression of an AuxRE reporter gene is reduced in brx roots (Mouchel et al., 2006); (ii) the observation that brx roots mimic wild-type roots treated with auxin transport inhibitors (Scacchi et al., 2009), consistent with a downregulation of certain PIN genes in brx (Mouchel et al., 2006); Only models 5.1.1 and 5.2.1 matched these observations.

The principal difference between these two models is in the regulatory relationship between BRX and NGA1, which is reflected in a predicted decrease (5.1.1) or increase (5.2.1) in auxin biosynthesis and thus free auxin levels upon BRX loss of function. Surprisingly, direct measurement indicated that free auxin levels are indeed significantly increased in brx seedlings (Fig. 3B), suggesting that model 5.2.1 is most parsimonious with experimental observations.

To further evaluate the robustness of our model, we also ran it in Boolean, synchronous mode using the GinSim software (Gonzalez et al., 2006), which confirmed the cyclic attractors observed in the continuous form (see Supplementary Material). Finally, we also conducted a stability analysis by sampling a range of parameter values and time points. This analysis corroborated earlier findings about sensible parameter ranges (Mendoza and Pardo, 2010) and revealed very little variation in the qualitative behavior of our model, confirming that it is robust (see heat maps in the Supplementary Material).

3.4 Conclusions

In summary, we have built qualitative and continuous logical models of auxin and brassinosteroid signaling, and their interaction. The models account for well established as well as novel experimental observations and demonstrate how a cell type-specific modulator, BRX, could be integrated. Taking into account novel experimental data also allowed us to choose between alternative models to correctly place BRX in the network. This most parsimonious model predicts an inhibition of NGA1 by BRX and an accessory role of BRX in auxin perception and signaling, guiding future experimentation. Notably, no valid developmental output was observed when a direct or indirect (e.g. through NGA1), but ARF-independent impact of BRX on AuxREs was assumed. The recent experimental confirmation that BRX can interact directly with ARFs was in fact motivated by this prediction (Scacchi et al., 2010). Finally, our models can serve as a basis for the work of others and, thanks to the SQUAD approach, can be extended for quantitative, kinetic data as they become available.

ACKNOWLEDGEMENTS

We would like to thank Dr J. Nemhauser for arf2 seeds. M.S. performed the modeling and the development of SQUAD add-ons. K.S.O. and B.G. generated and analyzed arf2 brx mutants. J.R. and D.T. performed hormone measurements. M.S., I.X. and C.S.H. designed the research and wrote the article together with M.S.

Funding: Swiss National Science Foundation SystemsX ‘Plant growth in a changing environment’ project.

Conflict of Interest: none declared.

REFERENCES

Badescu
G.O.
Napier
R.M.
Receptors for auxin: will it all end in TIRs?
Trends Plant Sci.
 , 
2006
, vol. 
11
 (pg. 
217
-
223
)
Benjamins
R.
Scheres
B.
Auxin: the looping star in plant development
Annu. Rev. Plant Biol.
 , 
2008
, vol. 
59
 (pg. 
443
-
465
)
Beuchat
J.
, et al.  . 
BRX promotes Arabidopsis shoot growth
New Phytol.
 , 
2010
, vol. 
188
 (pg. 
23
-
29
)
Braun
N.
, et al.  . 
Conditional repression of AUXIN BINDING PROTEIN1 reveals that it coordinates cell division and cell expansion during postembryonic shoot development in Arabidopsis and tobacco
Plant Cell
 , 
2008
, vol. 
20
 (pg. 
2746
-
2762
)
Dharmasiri
N.
, et al.  . 
The F-box protein TIR1 is an auxin receptor
Nature
 , 
2005
, vol. 
435
 (pg. 
441
-
445
)
Diaz
J.
Alvarez-Buylla
E.R.
Information flow during gene activation by signaling molecules: ethylene transduction in Arabidopsis cells as a study system
BMC Syst. Biol.
 , 
2009
, vol. 
3
 pg. 
48
 
Di Cara
A.
, et al.  . 
Dynamic simulation of regulatory networks using SQUAD
BMC Bioinformatics
 , 
2007
, vol. 
8
 pg. 
462
 
Digiuni
S.
, et al.  . 
A competitive complex formation mechanism underlies trichome patterning on Arabidopsis leaves
Mol. Syst. Biol.
 , 
2008
, vol. 
4
 pg. 
217
 
Gonzalez
A.G.
, et al.  . 
GINsim: a software suite for the qualitative modelling, simulation and analysis of regulatory networks
Biosystems
 , 
2006
, vol. 
84
 (pg. 
91
-
100
)
Grieneisen
V.A.
, et al.  . 
Auxin transport is sufficient to generate a maximum and gradient guiding root growth
Nature
 , 
2007
, vol. 
449
 (pg. 
1008
-
1013
)
Guilfoyle
T.J.
Hagen
G.
Auxin response factors
Curr. Opin. Plant. Biol.
 , 
2007
, vol. 
10
 (pg. 
453
-
460
)
Hardtke
C.S.
Transcriptional auxin-brassinosteroid crosstalk: who's talking?
Bioessays
 , 
2007
, vol. 
29
 (pg. 
1115
-
1123
)
Hyduke
D.R.
Palsson
B.O.
Towards genome-scale signalling-network reconstructions
Nat. Rev. Genet.
 , 
2010
, vol. 
11
 (pg. 
297
-
307
)
Kepinski
S.
Leyser
O.
The Arabidopsis F-box protein TIR1 is an auxin receptor
Nature
 , 
2005
, vol. 
435
 (pg. 
446
-
451
)
Kieffer
M.
, et al.  . 
Defining auxin response contexts in plant development
Curr. Opin. Plant. Biol.
 , 
2010
, vol. 
13
 (pg. 
12
-
20
)
Kim
T.W.
, et al.  . 
Brassinosteroid signal transduction from cell-surface receptor kinases to nuclear transcription factors
Nat. Cell. Biol.
 , 
2009
, vol. 
11
 (pg. 
1254
-
1260
)
Kim
T.W.
Wang
Z.Y.
Brassinosteroid signal transduction from receptor kinases to transcription factors
Annu. Rev. Plant Biol.
 , 
2010
, vol. 
61
 (pg. 
23.21
-
23.24
)
Klamt
S.
, et al.  . 
A methodology for the structural and functional analysis of signaling and regulatory networks
BMC Bioinformatics
 , 
2006
, vol. 
7
 pg. 
56
 
Kuppusamy
K.T.
, et al.  . 
Cross-regulatory mechanisms in hormone signaling
Plant Mol. Biol.
 , 
2008
, vol. 
69
 (pg. 
375
-
381
)
Leyser
O.
Auxin distribution and plant pattern formation: how many angels can dance on the point of PIN?
Cell
 , 
2005
, vol. 
121
 (pg. 
819
-
822
)
Liu
J.
, et al.  . 
Modelling and experimental analysis of hormonal crosstalk in Arabidopsis
Mol. Syst. Biol.
 , 
2010
, vol. 
6
 pg. 
373
 
Locke
J.C.
, et al.  . 
Experimental validation of a predicted feedback loop in the multi-oscillator clock of Arabidopsis thaliana
Mol. Syst. Biol.
 , 
2006
, vol. 
2
 pg. 
59
 
Mendoza
L.
Pardo
F.
A robust model to describe the differentiation of T-helper cells
Theory Biosci.
 , 
2010
, vol. 
129
 (pg. 
283
-
293
)
Mendoza
L.
Xenarios
I.
A method for the generation of standardized qualitative dynamical systems of regulatory networks
Theor. Biol. Med. Model.
 , 
2006
, vol. 
3
 pg. 
13
 
Mouchel
C.F.
, et al.  . 
Natural genetic variation in Arabidopsis identifies BREVIS RADIX, a novel regulator of cell proliferation and elongation in the root
Genes Dev.
 , 
2004
, vol. 
18
 (pg. 
700
-
714
)
Mouchel
C.F.
, et al.  . 
BRX mediates feedback between brassinosteroid levels and auxin signalling in root growth
Nature
 , 
2006
, vol. 
443
 (pg. 
458
-
461
)
Muller
B.
Sheen
J.
Cytokinin and auxin interaction in root stem-cell specification during early embryogenesis
Nature
 , 
2008
, vol. 
453
 (pg. 
1094
-
1097
)
Naldi
A.
, et al.  . 
Diversity and plasticity of Th cell types predicted from regulatory network modelling
PLoS Comput. Biol.
 , 
2010
, vol. 
6
 pg. 
e1000912
 
Nemhauser
J.L.
, et al.  . 
Different plant hormones regulate similar processes through largely nonoverlapping transcriptional responses
Cell
 , 
2006
, vol. 
126
 (pg. 
467
-
475
)
Nemhauser
J.L.
, et al.  . 
Interdependency of brassinosteroid and auxin signaling in Arabidopsis
PLoS Biol.
 , 
2004
, vol. 
2
 pg. 
E258
 
Osmont
K.S.
, et al.  . 
Hidden branches: developments in root system architecture
Annu. Rev. Plant Biol.
 , 
2007
, vol. 
58
 (pg. 
93
-
113
)
Pandey
S.
, et al.  . 
Boolean modeling of transcriptome data reveals novel modes of heterotrimeric G-protein action
Mol. Syst. Biol.
 , 
2010
, vol. 
6
 pg. 
372
 
Philippi
N.
, et al.  . 
Modeling system states in liver cells: survival, apoptosis and their modifications in response to viral infection
BMC Syst. Biol.
 , 
2009
, vol. 
3
 pg. 
97
 
Sabatini
S.
, et al.  . 
An auxin-dependent distal organizer of pattern and polarity in the Arabidopsis root
Cell
 , 
1999
, vol. 
99
 (pg. 
463
-
472
)
Sanchez-Corrales
Y.E.
, et al.  . 
The Arabidopsis thaliana flower organ specification gene regulatory network determines a robust differentiation process
J. Theor. Biol.
 , 
2010
, vol. 
264
 (pg. 
971
-
983
)
Sauer
M.
, et al.  . 
Canalization of auxin flow by Aux/IAA-ARF-dependent feedback regulation of PIN polarity
Genes Dev.
 , 
2006
, vol. 
20
 (pg. 
2902
-
2911
)
Savaldi-Goldstein
S.
, et al.  . 
The epidermis both drives and restricts plant shoot growth
Nature
 , 
2007
, vol. 
446
 (pg. 
199
-
202
)
Scacchi
E.
, et al.  . 
Dynamic, auxin-responsive plasma membrane-to-nucleus movement of Arabidopsis BRX
Development
 , 
2009
, vol. 
136
 (pg. 
2059
-
2067
)
Scacchi
E.
, et al.  . 
Spatio-temporal sequence of cross-regulatory events in root meristem growth
Proc. Natl Acad. Sci. USA
 , 
2010
, vol. 
107
 (pg. 
22734
-
22739
)
Scheres
B.
Xu
J.
Polar auxin transport and patterning: grow with the flow
Genes Dev.
 , 
2006
, vol. 
20
 (pg. 
922
-
926
)
Sibout
R.
, et al.  . 
Opposite root growth phenotypes of hy5 versus hy5 hyh mutants correlate with increased constitutive auxin signaling
PLoS Genet.
 , 
2006
, vol. 
2
 pg. 
e202
 
Strader
L.C.
, et al.  . 
The IBR5 phosphatase promotes Arabidopsis auxin responses through a novel mechanism distinct from TIR1-mediated repressor degradation
BMC Plant Biol.
 , 
2008
, vol. 
8
 pg. 
41
 
Tiwari
S.B.
, et al.  . 
The roles of auxin response factor domains in auxin-responsive transcription
Plant Cell
 , 
2003
, vol. 
15
 (pg. 
533
-
543
)
Trigueros
M.
, et al.  . 
The NGATHA genes direct style development in the Arabidopsis gynoecium
Plant Cell
 , 
2009
, vol. 
21
 (pg. 
1394
-
1409
)
Tromas
A.
, et al.  . 
The AUXIN BINDING PROTEIN 1 is required for differential auxin responses mediating root growth
PLoS ONE
 , 
2009
, vol. 
4
 pg. 
e6648
 
Vert
G.
, et al.  . 
Molecular mechanisms of steroid hormone signaling in plants
Annu. Rev. Cell. Dev. Biol.
 , 
2005
, vol. 
21
 (pg. 
177
-
201
)
Vert
G.
, et al.  . 
Integration of auxin and brassinosteroid pathways by Auxin Response Factor 2
Proc. Natl Acad. Sci. USA
 , 
2008
, vol. 
105
 (pg. 
9829
-
9834
)
Vieten
A.
, et al.  . 
Functional redundancy of PIN proteins is accompanied by auxin-dependent cross-regulation of PIN expression
Development
 , 
2005
, vol. 
132
 (pg. 
4521
-
4531
)
Wang
S.
, et al.  . 
AUXIN RESPONSE FACTOR7 restores the expression of auxin-responsive genes in mutant Arabidopsis leaf mesophyll protoplasts
Plant Cell
 , 
2005
, vol. 
17
 (pg. 
1979
-
1993
)
Weijers
D.
, et al.  . 
Developmental specificity of auxin response by pairs of ARF and Aux/IAA transcriptional regulators
EMBO J.
 , 
2005
, vol. 
24
 (pg. 
1874
-
1885
)
Wisniewska
J.
, et al.  . 
Polar PIN localization directs auxin flow in plants
Science
 , 
2006
, vol. 
312
 pg. 
883
 
Wittmann
D.M.
, et al.  . 
Transforming Boolean models to continuous models: methodology and application to T-cell receptor signaling
BMC Syst. Biol.
 , 
2009
, vol. 
3
 pg. 
98
 

Author notes

Associate Editor: Alfonso Valencia

Comments

0 Comments