Mathematical Modelling and Analysis of the Brassinosteroid and Gibberellin Signalling Pathways and their Interactions

The plant hormones brassinosteroid (BR) and gibberellin (GA) have important roles in a wide range of processes involved in plant growth and development. In this paper we derive and analyse new mathematical models for the BR signalling pathway and for the crosstalk between the BR and GA signalling pathways. To analyse the effects of spatial heterogeneity of the signalling processes, along with spatially-homogeneous ODE models we derive coupled PDE-ODE systems modelling the temporal and spatial dynamics of molecules involved in the signalling pathways. The values of the parameters in the model for the BR signalling pathway are determined using experimental data on the gene expression of BR biosynthetic enzymes. The stability of steady state solutions of our mathematical model, shown for a wide range of parameters, can be related to the BR homeostasis. Our results for the crosstalk model suggest that the interaction between transcription factors BZR and DELLA exerts more influence on the dynamics of the signalling pathways than BZR-mediated biosynthesis of GA, suggesting that the interaction between transcription factors may constitute the principal mechanism of the crosstalk between the BR and GA signalling pathways. In general, perturbations in the GA signalling pathway have larger effects on the dynamics of components of the BR signalling pathway than perturbations in the BR signalling pathway on the dynamics of the GA pathway. The perturbation in the crosstalk mechanism also has a larger effect on the dynamics of the BR pathway than of the GA pathway. Large changes in the dynamics of the GA signalling processes can be observed only in the case where there are disturbances in both the BR signalling pathway and the crosstalk mechanism.


a b s t r a c t
The plant hormones brassinosteroid (BR) and gibberellin (GA) have important roles in a wide range of processes involved in plant growth and development. In this paper we derive and analyse new mathematical models for the BR signalling pathway and for the crosstalk between the BR and GA signalling pathways. To analyse the effects of spatial heterogeneity of the signalling processes, along with spatiallyhomogeneous ODE models we derive coupled PDE-ODE systems modelling the temporal and spatial dynamics of molecules involved in the signalling pathways. The values of the parameters in the model for the BR signalling pathway are determined using experimental data on the gene expression of BR biosynthetic enzymes. The stability of steady state solutions of our mathematical model, shown for a wide range of parameters, can be related to the BR homeostasis which is essential for proper function of plant cells. Solutions of the mathematical model for the BR signalling pathway can exhibit oscillatory behaviour only for relatively large values of parameters associated with transcription factor brassinazole-resistant1's (BZR) phosphorylation state, suggesting that this process may be important in governing the stability of signalling processes. Comparison between ODE and PDE-ODE models demonstrates distinct spatial distribution in the level of BR in the cell cytoplasm, however the spatial heterogeneity has significant effect on the dynamics of the averaged solutions only in the case when we have oscillations in solutions for at least one of the models, i.e. for possibly biologically not relevant parameter values. Our results for the crosstalk model suggest that the interaction between transcription factors BZR and DELLA exerts more influence on the dynamics of the signalling pathways than BZR-mediated biosynthesis of GA, suggesting that the interaction between transcription factors may constitute the principal mechanism of the crosstalk between the BR and GA signalling pathways. In general, perturbations in the GA signalling pathway have larger effects on the dynamics of components of the BR signalling pathway than perturbations in the BR signalling pathway on the dynamics of the GA pathway. The perturbation in the crosstalk mechanism also has a larger effect on the dynamics of the BR pathway than of the GA pathway. Large changes in the dynamics of the GA signalling processes can be observed only in the case where there are disturbances in both the BR signalling pathway and the crosstalk mechanism. Those results highlight the robustness of the GA signalling pathway.

Introduction
The sessile nature of plant life highlights the importance of efficient regulatory mechanisms allowing plants to respond to environmental stimuli and to adapt to changing environmental conditions. Plants have developed a set of highly integrated signalling pathways. Plant hormones, e.g. auxin, gibberellin, cytokinin, brassi-nosteroids, ethylene, are key signalling molecules and their activities depend on the cellular context and interactions between them.
The family of steroidal plant hormones brassinosteroids (BRs) is responsible for the regulation and control of a wide range of essential processes including responses to stresses ( Bajguz and Hayat, 2009;Gruszka, 2013 ), photomorphogenesis ( Belkhadir and Jaillais, 2015;Zhu et al., 2013 ), root growth ( Müssig et al., 2003 ), and stomatal development ( Kim et al., 2012 ). Over the last 47 years, the effects of brassinosteroids on plant cells and plants as a whole, as well as their signalling pathways have been studied in detail ( Clouse, 2015 ). In particular the signalling pathway of brassinolide (BL), the most biologically active of the discovered BRs, has http://dx.doi.org/10.1016/j.jtbi.2017.08.013 0022-5193/© 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ ) been examined in great detail, and is now one of the most understood pathways in plant biology ( Belkhadir and Jaillais, 2015;Clouse, 2011;Kim and Wang, 2010;Li and Deng, 2005;She et al., 2011;Wang et al., 2014c;Yang et al., 2011;Zhu et al., 2013 ). BR signalling functions by controlling the expression of various genes regulating developmental processes, of which 10 0 0s have been identified ( Sun et al., 2010 ), including several genes which regulate the production of proteins that act as enzymes during BR biosynthesis . To ensure controlled growth, homeostasis of BRs in plant tissue is carefully maintained by negative feedback in the BR signalling pathway . Absence of BRs and/or perturbed BR signalling have also been linked to many growth defects, including dwarfism and male sterility ( Clouse, 1996;Clouse and Sasse, 1998 ).
Gibberellins (GAs) are another family of plant hormones involved in many developmental processes in plants, including seed germination, stem elongation, leaf expansion, trichrome development, pollen maturation and the induction of flowering ( Achard et al., 2008;Davière and Achard, 2013 ). There are over 130 categorized gibberellins, a few of which are bioactive molecules, the most common being GA 1 , GA 3 and GA 4 ( Yamaguchi, 2008 ).
Along with detailed information on individual plant hormone signal transduction pathways, their target genes, and their effects, it is now also known that interactions between various molecules involved in different signalling pathways have an effect on physiological phenomena in plants ( Bai et al., 2012;Belkhadir and Jaillais, 2015 ). For example both auxin and brassinosteroids play a role in the patterning of vascular shoot bundles ( Ibañes et al., 0 0 0 0 ), as well as exhibiting some cross-regulation of biosynthesis via the BRX gene ( Sankar et al., 2011 ). The GA signalling pathway exhibits a high level of interaction with the BR signalling pathway regulating growth ( Clouse and Sasse, 1998;Úbeda Tomás et al., 2009 ) and responses to stresses ( Achard et al., 2008;Bajguz and Hayat, 2009 ) among others. The interactions between different signalling pathways are called crosstalks, and the investigation of their mechanisms is key to better understand plant hormonal responses to external stimuli ( Albrecht et al., 2012;Belkhadir and Jaillais, 2015;Gruszka, 2013;Yang et al., 2011 ). Despite the need for better understanding of the mechanisms of crosstalk between hormone signalling pathways, it is hard to obtain experimentally quantitative data on the dynamics of all molecules involved in the signalling pathways and interactions between them. For example, it is very difficult to measure the dynamics of hormones such as BR in real time, in part due to their occurring naturally at extremely low levels. Therefore development of accurate mathematical models of hormone activity can help to analyse and better understand interactions between signalling pathways and their impact on plant growth and development. Hence, the main aim of this paper is to derive and analyse novel mathematical models for the BR signalling pathway, and the crosstalk between the BR and GA signalling pathways.
Along with many modelling results for the gibberellin and auxin signalling pathways ( Gordon et al., 2009;Liu et al., 2010;Middleton et al., 2010;2012;Muraro et al., 2011 ), only a few models can be found for BR-related signalling processes. A logic model for the activation states of the components of the BR and auxin signalling pathways and their interactions was derived in Sankar et al. (2011) , and provided a qualitative description of the dynamics of the BR pathway. In a simple model for BR-mediated root growth, proposed in van Esse et al. (2012) , the growth dynamics is assumed to be dependent on the quantity of receptor-bound BL, which was considered to be constant. A system of ordinary differential equations was considered to describe the dynamics of BRregulated transcription factor BES1 (BRI1-EMS SUPRESSOR 1) and its interactions with R2R3-MYB transcription factor BRAVO, related to division of plant stem cells ( Frigola et al., 2017;Vilarrasa-Blasi et al., 2014 ). To our knowledge there are no previous results on mathematical modelling of the crosstalk between BR and GA signalling pathways. In our mathematical model for the BR signalling pathway, we consider the dynamics of BR, free and bound receptors, inhibitors, and phosphorylated and dephosphorylated transcription factors, not considered in previous models. The spatially homogeneous dynamics of the molecules involved in the signalling pathway that we consider are modelled by a system of six ordinary differential equations (ODEs). To analyse the effect of spatial heterogeneity of the signalling processes on the dynamics of BR, we derive a coupled model composed of partial differential equations (PDEs) for BR, inhibitor, and phosphorylated transcription factor, ODEs for receptors, defined on the cell membrane, and the ODE for dephosphorlyated transcription factors, localised in the cell nucleus. Along with spatial distribution of the concentration of BR in the cell cytoplasm, we observe similar dynamics for solutions of the ODE and averaged solutions of the PDE-ODE models when those solutions converge to a steady state as t → ∞ . However spatial heterogeneity has significant effect in the case when at least one of the two models has periodic solutions, which is determined for possibly biologically not relevant parameter values.
To model the crosstalk between BR and GA signalling pathways we first rigorously derive a reduced model for the GA signalling pathway from the full GA signalling pathway model proposed in Middleton et al. (2012) . Then we couple the model for the BR signalling pathway with the reduced model for GA signalling pathway by considering three different interaction mechanisms between BR and GA pathways. By analysing the effect of different interaction mechanisms on the dynamics of molecules involved in the signalling processes, we determine that one of these mechanisms has a more significant effect on the dynamics of the signalling pathways than the other mechanisms. Similar to the BR signalling pathway model, we also consider the influence of spatial heterogeneity of the signalling processes on the dynamics of solutions of the BR-GA crosstalk model. Using the mathematical models developed here, we analyse how interactions between the BR and GA signalling pathways depend on the model parameters and the strength of interaction mechanisms. We observed that, in general, parameter changes in both pathways have a stronger effect on the components of the BR signalling pathway than on the components of the GA signalling pathway. Our results also suggest that the interaction between transcription factors exerts more influence on the dynamics of the signalling pathways than BR signallingmediated GA biosynthesis. Further, our results suggest that perturbations in the GA signalling pathway have larger effects on the dynamics of components of the BR signalling pathway than the perturbations in the BR signalling pathway on the dynamics of components of the GA signalling pathway, apart from in the case when we have disturbances in both the BR signalling pathway and the crosstalk mechanism.
The structure of this paper is as follows. In Section 2 a biological overview of the BR and GA signalling pathways, and their interactions is given. In Section 3 we derive the mathematical model for the BR signalling pathway and estimate the values for the model parameters using experimental data from Tanaka et al. (2005) . In Section 4 we perform qualitative analysis of the model for the BR signalling pathway, examining how the behaviour of solutions of the mathematical model depends on the values of the model parameters. We also define the set of parameters for which the system of ODEs has stable stationary solutions and the set of parameters for which it undergoes Hopf bifurcation. In Section 5 we extend our model for the BR signalling pathway to examine the effects of spatial heterogeneity in the signalling processes. In Section 6 we consider the reduction of the mathematical model for the GA signalling pathway, proposed in Middleton et al. (2012) , derive a new model for the crosstalk between BR and GA signalling pathways, and analyse the influence of different interaction mechanisms and changes in the dynamics of one of the signalling pathways on the qualitative and quantitative behaviours of the coupled system. We also analyse the influence of spatial heterogeneity of the signalling processes on the interactions between the BR and GA signalling pathways. We summarise and discuss our results in Section 7 .

Biological background
The mathematical modelling and analysis of the BR signalling pathway and of the interactions between the BR and GA signalling pathways is the main aim of this paper. In this section we present an overview of the BR and GA signalling pathways, and the interactions between them.

The BR signalling pathway
The signalling process starts at the cell plasma-membrane with the perception of BR by the receptor BRASSINOSTEROID INSENSI-TIVE1 (BRI1) ( Li and Chory, 1997 ). Upon BR binding to BRI1, two main events then take place, association of BRI1 to a co-receptor, BRI1-ASSOCIATED RECEPTOR KINASE1 (BAK1), and dissociation of the inhibitor protein BRI1 KINASE INHIBITOR1 (BKI1). This triggers a transphosphorylation cascade between BRI1 and BAK1, leading further to phosphorylation of BRASSINOSTEROID-SIGNALLING KI-NASE1 (BSK1), another membrane-bound kinase. Next, BSK1 phosphorylates the protein phosphatase BRI1-SUPPRESSOR1 (BSU1), which dephosphorylates a protein kinase BRASSINOSTEROID IN-SENSITIVE2 (BIN2), eventually leading to its degradation ( Ryu et al., 2010 ). In the absence of BR, phosphorylated BIN2 has a role in phosphorylating the two transcription factors, BRASSINAZOLE RESISTANT1 (BZR1) and BRI1-EMS-SUPPRESSOR1 (BES1) ( Li and Deng, 2005 ), also known as BZR2 (for the purposes of this paper it is unnecessary to distinguish between the two, so we refer to them jointly as BZR). When phosphorylated, BZR is less stable and thus more unlikely to activate or repress any of the 10 0 0s of genes associated with BR signalling ( Ryu et al., 2007 ), it is also thought that association of phosphorylated BZR to a 14-3-3 protein inhibits its entry to the nucleus. PROTEIN PHOSPHOTASE 2A (PP2A) is responsible for the de-phosphorylation of BZR, which allows its entry into the nucleus and then its activation of BR responsive genes.
BZR functions as a repressor of certain genes associated with the biosynthesis of BR, notably for example CONSTITUTIVE PHO-TOMORPHOGENESIS AND DWARFISM (CPD), DWARF4 (DWF4), ROTUNDIFOLIA3 (ROT3) and BRASSINOSTEROID-6-OXIDASE 1 (BR6ox1) . That is, active, de-phosphorylated BZR inhibits production of BR. So, high levels of BR cause low levels of phosphorylated BZR, leading to inhibition of BR biosynthesis and decreasing levels of BR. Conversely, low levels of BR lead to high levels of phosphorylated BZR and activation of BR biosynthesis, increasing the levels of BR. This completes the negative feedback loop of the BR signalling pathway.

The GA signalling pathway
Gibberellin Signalling is achieved by enhancing the degradation of DELLA proteins, which influence the expression of GAresponsive genes ( Achard and Genschik, 2009 ). GA molecules are perceived by the GA receptor, GIBBERELLIN INSENSITIVE DWARF1 (GID1), a nuclear-localised protein ( Ueguchi-Tanaka et al., 2005 ). Analysis of GID1's structure revealed that it has a GA-binding pocket, with a flexible extension adjacent ( Shimada et al., 2008 ). When GA binds to GID1, this extension undergoes conformational change, and covers the GA-binding pocket. When closed, the upper surface of this lid binds to DELLA proteins to form the GA.GID1.DELLA complex. The formation of the GA.GID1.DELLA complex enhances the degradation of DELLA proteins by mediating proteasome-dependent destabilization of DELLA proteins.
The GA signalling pathway exhibits negative feedback due to the influence of DELLAs on the expression of several genes which code components of the signalling pathway ( Yamaguchi, 2008 ). First, DELLA activates the GID1-encoding gene, leading to an increase in the translation of the GID1 protein. This means that in the absence of DELLAs, GID1 concentration also decreases which will slow down the proteasome-induced DELLA degradation, and that an abundance of DELLA leads to the production of more GID1 and enhances the DELLA degradation. Next, DELLA activates the transcription of genes encoding the enzymes GA 20-oxidase (GA20ox) and GA 3-oxidase (GA3ox). These enzymes catalyse several reaction steps in the GA biosynthesis pathway, meaning an abundance of DELLA increases both GA and GID1, leading to degradation of DELLA. Lastly, DELLA represses its own gene transcription.

Crosstalk between the BR and GA signalling pathways
The interaction of BRs and GAs has been receiving much attention, due to their shared nature as critical plant growth regulators, combined with the fact that they share many overlapping functions such as regulation of cell elongation ( Catterou et al., 2001;Úbeda Tomás et al., 2009 ) and plant responses to abiotic stress ( Ahammed et al., 2015;Colebrook et al., 2014 ). However despite the interest, the exact mechanisms of these interactions have remained largely unclear, save from the fact that they control expression of several genes ( Bouquin et al., 2001 ). There has been much evidence that the signalling processes of BR and GA converge at the level of BZR and DELLA interaction. Direct crosstalk in this fashion was shown in Li et al. (2012) , where it was shown that overexpression of DELLA proteins reduced both the abundance and transcriptional activity of BZR. This was found to be due to the formation of a complex of DELLA and BZR, which removed BZR's transcriptional ability. There is also evidence of BRs regulating the biosynthesis of GAs, the so called "GA Synthesis" model of BR-GA crosstalk. Two main proposals have been made for the existence of this type of interaction. The effects of BR mutants on GA synthesis were examined in Tong et al. (2014) , and it was concluded that BZR enhances GA synthesis by activating synthesis of the GA3ox enzyme. In contrast to this the findings in Unterholzner et al. (2015) describe a much larger role for BZR in regulating GA synthesis. They provide evidence for a model where BZR activates the synthesis of the GA20ox enzyme, in addition to the effects described in Tong et al. (2014) . Thus BZR would influence the biosynthesis of GA in exactly the same manner as DELLAs for other interactions, however the significance of this mechanism of crosstalk is not yet established ( Ross and Quittenden, 2016;Tong and Chu, 2016;Unterholzner et al., 2016 ).

Derivation of a mathematical model for the BR signalling pathway
In this section we derive a mathematical model of the BR signalling pathway. In order to build a simple, yet sufficiently accurate and efficient model incorporating BR biosynthesis negative feedback, we first construct a reduced reaction schematic that describes the pathway mechanism. This reduction is achieved via a simplification of two principal parts of the signalling pathway: the complex BR biosynthesis network, and the cytoplasm localized phosphorylation cascade. Hence, we build a model focussing on three key components, hormone (BR), inhibitor (BKI1) and transcription factor (BZR) ( Fig. 1 ). In the mathematical model we consider the binding of free BR molecules to the BRI1 receptors leading to the dissociation of BKI1. This is modelled as an almost instantaneous reaction, with BR + BRI1.BKI1 interacting and resulting into BR.BRI1 + BKI1. In order to model the effects of the signalling cascade induced by the membrane-bound receptors and subsequent effects on the phosphorylation state of BZR, we assume that the effects of active receptors in triggering the cascade may be approximated by the free BKI1 that is released upon this activation. We further assume that the free BKI1 catalyses dephosphorylation of BZR-p and destabilises the BIN2 proteins, thus reducing the phosphorylation of BZR.
We denote by b the concentration of hormone BR, by k the concentration of inhibitor BKI1, by r k the concentration of receptorinhibitor complex BRI1.BKI1, by r b the concentration of receptorhormone complex BR.BRI1, by z the concentration of (dephosphorylated) transcription factor BZR, and by z p the concentration of (phosphorylated) transcription factor BZR-p. Then assuming spatial homogeneity of the signalling processes, the interactions between b, k, r k , r b , z , and z p are described by the system of six ordinary differential equations db dt Here β b is the binding rate of b to r k , and β k is the binding rate of k to r b . We model this as only two reactions by assuming that when either BR or BKI1 are bound to BRI1 either dissociates sufficiently fast that the levels of the BR.BRI1.BKI1 remains roughly zero. We assume the reaction to occur in some finite closed volume, so the loss of BR is only described by the degradation coeffi- The phosphorylation state of BZR is modelled as being dependent on the levels of free BKI1. We justify this by noting that upon signalling activation, the receptor phosphorylates BSU1 which governs the phosphorylation of BZR via BIN2. Thus since BKI1 is also released upon BR binding, we may model these effects by assuming that free BKI1 activates (or catalyses) the dephosphorylation of BZR-p at rate δ z . In BKI1's absence BZR is phosphorylated at rate ρ z , and when BKI1 is present it inhibits the phosphorylation of BZR such that when k = 1 /θ z , the rate of phosphorylation is halved.
Finally we model BR biosynthesis as being directly inhibited by BZR. BZR represses the expression of several genes encoding enzymes, namely CPD, DWF4, ROT3 and BR6ox1, that are required for the conversion of many of the precursors involved in BR biosynthesis. Hence we use a Hill function with exponent h b > 1 to model the cumulative inhibitory effect of BZR on these genes. We estimate the parameter h b by considering the detailed BR biosynthesis pathway(s) presented in Chung and Choe (2013) . The BR-mediated enzymes that are involved in the biosynthesis are CPD, DWF4, ROT3, and BR6ox1 which mediate four, five, six, and five steps in the biosynthetic pathway respectively. Since BR biosynthetic enzymes act multiplicatively at different steps of the reaction network, the expressions modelling their actions can be approximated by the product of these expressions, which would mean that the exponents of these functions would be summed, thus we consider h b = 20 .
The model equations (1) imply that the total concentrations of BKI1, BRI1 and BZR are conserved, thus we consider k + r k = K tot , r b + r k = R tot , and z + z p = Z tot , and derive a reduced model: Various values for R tot were reported in van Esse et al. (2011) , and we chose the value for WT seedling roots. We further assume that K tot = R tot in order that the receptor should have the ability to be completely inactive, but not be saturated by BKI1. We also use physiological values reported in the literature in order to write β k and μ b in terms of other parameters, for full calculations see Appendix A . As such we are left with 8 parameters for which we have no direct estimate, namely β b , α b , θ b , δ z , Z tot , ρ z , θ z and h z . These parameters were estimated indirectly by validating the numerical solutions of the mathematical model (2) against experimental results, using numerical optimisation algorithms.
By deriving the steady state concentration of BR, denoted [ BR ] 0 , from the level of endogenous 24-epiBL reported in Wang et al.
(2014b ), we were able to write the rate of BR degradation μ b in terms of α b , θ b , Z tot , δ z , ρ z , θ z , h z and h b as follows We can write β k in terms of β b and other known parameters in two ways. First, using [ BR ] These two constraints (4) and (5) (2) , plotted against experimental data from Tanaka et al. (2005) . For the WT plants grown under control conditions the parameters given in Table 1 were used, parameter sets for other cases can be found in Table A.6 .  Table 2 were used.
(2) was fitted to experimental data using both conditions (4) and (5) in order to compare their effects, see Figs. 2 and 3 . The experimental data from Tanaka et al. (2005) , used to determine model parameters, give the relative BR biosynthetic gene expression of CPD, DWF4, ROT3, and BR6ox1, and were measured by RT-PCR analysis, then converted to give relative values with the initial values equal to one. Values were measured for three independent experiments (three replicates), and the data presented by points and (where available) error bars correspond to the mean and standard error respectively. Gene expression was measured in both Wild-Type (WT) and bri1-401 mutant (where perception of BR by BRI1 is inhibited) plants, and this was accounted for in our parameter estimation by assuming that the parameter β b was greater for the WT than for the mutant. β k was allowed to vary freely for the mutant case since constraints (4) and (5) are not def-initely valid in this case. Both of these phenotypes were grown under control conditions as well as two other cases: one where plants were grown in a medium containing 5μM of Brassinazole (BRZ), a BR-specific biosynthesis inhibitor, and one grown in a medium containing 0.1 μM of Brassinolide (BL) having first been grown in the medium containing 5 μM BRZ for two days. Data comparing the control case with the case of growth in the 5 μM BRZ medium were recorded for five days. Data comparing further growth after two days of the BRZ medium case with the case of addition of 0.1 μM BL to the medium were recorded for further 24 h, and as such for these cases the initial conditions were taken to be the values of the numerical solution to the model for the 5 μM BRZ medium at time t = 2 days. For the control conditions we made no amendments to the model, for the case of plants growing in the BRZ-medium we imposed bounds upon the parameters such that α b should be smaller in this case since addition of BRZ reduces the biosynthesis of BR. In order to examine the case where BL was added to the growth medium, an extra term governing influx of exogenous BL and efflux of endogenous BR was added to the equation describing BR dynamics db dt where φ b is the relative permeability of the cell membrane to BR and was one of the optimised parameters for the relevant cases, and ω b is assumed to be 0.1 μM in accordance with the experimental procedure.
Optimisation was achieved by comparing the biosynthetic expression defined by the numerical solutions of model (2) , Tanaka et al. (2005) . In order to compare experimental data with the output from our model we first normalised the simulation data by their initial values such that they took values comparable to the experimental data. We then took the mean of the four gene data sets, weighted by the number of times the respective proteins appear in the biosynthetic pathway. The optimisation was carried out in Python using the curve_fit function in the SciPy module ( Jones et al., 2001 ). curve_fit applies nonlinear least squares minimization using the trust region reflective algorithm as default, with a default tolerance of 10 −8 . The model was fitted to the data set for each case sequentially, starting with the WT under control conditions since this data set was the largest. The parameters μ b and β k were replaced by the expressions (3) and (4) or (5) , respectively, for the WT data under control conditions since this is the only case where such parameter constraints are definitely valid. The parameters generated from the fitting for WT under control conditions were then used as the initial guesses for all other cases, where μ b and β k were also allowed to be fitted. Parameters that were not expected to vary under the different growth conditions were allowed very small variations to account for error in the first fitting, whereas parameters that were expected to vary had much wider bounds. Numerical simulations of model (2) using the optimised parameters, given in Tables 1 and 2 for β k determined by (4) and (5) respectively, are plotted against the experimental data in Fig. 2 and 3 and show good agreement with experimental data, having R 2 values of 0.89 and 0.92 respectively.

Qualitative analysis of the mathematical model for the BR signalling pathway
In this section we consider well-posedness and qualitative analysis of model (2) . We start by non-dimensionalising our model, transforming the variables as t = 1 and introducing the dimensionless parameters Table 1 The model parameters associated with the BR signalling pathway model (2) for WT when fitted using expression (4) for β k , together with their sources (Parameter set 1).  Tanaka et al. (2005) Table 2 The model parameters associated with the BR signalling pathway model (2) for WT when fitted using expression (5) For simplicity of presentation we denote by P ⊂ [1 , ∞ ) × R 7 + × N 2 the parameter space for system (6) , where for each p ∈ P , p = (κ, β k , β b , θ b , , δ z , ρ z , θ z , h b , h z ) . We assume that κ has a minimum value of 1 since κ < 1 would imply saturation of receptor by inhibitor (i.e. K tot > R tot ), leading to BR signalling being permanently switched on.

Theorem 1. The system (6) has a unique, global solution
and any p ∈ P.
Since f is locally Lipschitz-continuous, the Picard-Lindelöf theorem ensures local existence of a unique solution of (6) , see e.g. ( Amann, 1990 ). To obtain global existence and uniqueness we prove boundedness of solutions by demonstrating the existence of a positive-invariant region for system (6) , i.e. showing that for a solution u of (6) starting in M = always be contained within this region. To show that a region M is positive invariant under the flow of system (6) , we show that ( Amann, 1990 ), where n is the inward normal vector on ∂M , see B.1 for more details. This implies uniform boundedness of solutions of (6) with initial values in M , and continuous differentiability of f ensures global existence and uniqueness.
Theorem 2. For any parameter set p ∈ P, there exists a unique steady Proof. Considering equations for a steady state solution ( b * , k * , z * ) of (6) and employing simple algebraic manipulation, k * is defined as a root of the following non-linear function and b * and z * are determined as functions of k * , for more details see ( B.2 ). We may immediately see, since g(0) = −β b and g(1) = at least one root in [0, 1] for any p ∈ P , and hence system (6) must contain at least one steady state solution in M . In order to find the number of roots of g(k * ) = 0 in [0, 1], consider the derivative of g which is positive for all p ∈ P and k * ∈ [0, 1], see ( B.2 ) for the formula for g ( k * ). Thus g is monotonically increasing. Strict monotonicity of g coupled with existence of at least one root in [0, 1] implies that g has a unique root in [0, 1], and hence (6) has a unique steady state in M .

Linearised stability and bifurcation analysis
To study the qualitative behaviour of solutions of mathematical model for BR signalling pathway, we performed linearised stability analysis for system (6) and analyse the impact of variations in values of model parameters on the behaviour of solutions of (6) . For the parameters obtained via validation of mathematical model by experimental data, see Tables 1 and 2 , steady state solutions of (6) are linearly stable, with eigenvalues  (2) and can induce oscillatory behaviour, but only in the case when all other parameters are as in Table 1 and not as in Table 2 . We consider δ z and ρ z as bifurcation parameters because these parameters directly correspond to processes in the BR signalling pathway, whereas is only an approximation for the dynamics of the cytoplasmic phosphorylation cascade. Biologically, increase of δ z could potentially correspond to faster phosphorylation of BSU1 by BAK1, or decrease of δ z corresponding to reduced action of PP2A in dephosphorylating BZR. Further, decrease of ρ z could correspond to BIN2-deficient or insensitive mutants e.g. bes1-D, and increase of ρ z could correspond to BIN2-overexpressing mutants, e.g. bin2.
In the bifurcation analysis of system (2) we considered increased value for the dimensional parameter θ z compared to the standard value, Table 1 (i.e. θ z = 41 . 2 μM −1 ), which was essential to determine the region for parameters δ z and ρ z where system (6) undergoes bifurcation. For the value of θ z = 3 . 95 μM −1 obtained through fitting model solutions to experimental data, the steady state solution of (6) is linearly stable for a wide range of values of δ z and ρ z , i.e. δ z , ρ z ∈ (0, 50).
Proof. We performed linearised stability analysis to determine the parameter subspace for which the stationary solution of (6) is linearly stable, as well as the range of parameters for which we have periodic solutions for the model (6) .
Using the Jacobian (B.1) of system (6) , evaluated at the steadystate ( b * , k * , z * ), we calculate the characteristic equation for the system to be a cubic polynomial of the form λ 3 + a 2 λ 2 + a 1 λ + a 0 = 0 , where a 2 , a 1 and a 0 are all positive, real constants (see B.3 ). Since the characteristic equation is a cubic polynomial we obtain that there are only 2 possible sets of eigenvalues, either that they are all real or that there is one real eigenvalue λ 1 and two complex conjugate eigenvalues λ 2 and λ 3 . Further, since all coefficients have the same sign any real eigenvalue must be negative, specifically zero cannot be an eigenvalue of J for any parameters of system (6) in P . Together these two facts tell us that in the case where the eigenvalues are all real, or that the complex conjugate eigenvalues have negative real part the steady state is stable, and that there is a possible bifurcation point when the two complex conjugate eigenvalues cross the imaginary axis, corresponding to a Hopf bifurcation.
We showed numerically that for a closed loop L in ( δ z , ρ z ) the system has complex conjugate eigenvalues with zero real parts, Fig. 4 , and in the region D enclosed by the loop L the complex conjugate eigenvalues have positive real part, whereas for (δ z , ρ z ) ∈ (0 , 50) 2 \ D the real part of the complex conjugate eigenvalues is negative. Hence the points of the loop L correspond to the bifurcation points where the stability of stationary solutions of (6) changes. We also showed that at such points the eigenvalues have non-zero imaginary parts and hence do not pass through the origin, Fig. 5 a), which supports the proof of the fact that zero cannot be an eigenvalue of J , presented above. For this we designed a scheme in MATLAB to calculate the eigenvalues of the Jacobian J in (B.1) , for values of ( δ z , ρ z ) ∈ (0, 50) 2 , with dimensional parameter values θ z = 41 . 2 μM −1 and all other parameters as in Table 1 . For each δ z ∈ (0, 50), we determined the values of ρ z ∈ (0, 50) for which J has a pair of non-zero purely imaginary eigenvalues. The derivatives of the real part of the eigenvalues w.r.t. both δ z and ρ z were also calculated numerically, Fig. 5 b). The values of d dδz Re (λ 2 , 3 ) and d dρz Re (λ 2 , 3 ) at the critical points are non-zero, apart from exactly four points on the curve in Fig. 5 b), where the derivative w.r.t. one of the parameters will be zero. Those four points correspond to the points where δ z or ρ z are at their extreme values, i.e. when δ z takes an extreme value we have that d dρz Re (λ 2 , 3 ) = 0 , and when ρ z takes an extreme value that d dδz Re (λ 2 , 3 ) = 0 . When δ z is fixed at an extreme point, varying ρ z will not cause the point to enter the region bounded by the curve and will not correspond to a bifurcation point w.r.t. ρ z , similar for ρ z fixed at an extreme point.
Hence zero derivative at this point does not break the transversality condition. Thus, in conjunction with Theorems 1 and 2 we have shown that system (6) satisfies all conditions for the existence of a local Hopf bifurcation ( Hassard et al., 1981 ).
Therefore, for all parameter sets such that (δ z , ρ z ) ∈ (0 , 50) 2 \ D , θ z = 41 . 2 μM −1 , and all other dimensional parameters are defined as in Table 1 , we have that the steady state solution of system (6) is linearly stable. At the points (δ z , ρ z ) ∈ L the system (6) undergoes a Hopf bifurcation, and for (δ z , ρ z ) ∈ D we have periodic solutions for the BR signalling pathway model.

Spatially heterogeneous model for the BR signalling pathway
The BR signalling pathway is a process which has distinct functions at different spatial locations in the cell ( Fig. 6 ). Thus it is important to extend the ODE model (1) and analyse the dependence of the dynamics of the pathway components on the spatial distribution.
We consider a spatially heterogeneous model for the BR signalling pathway in the one-dimensional domain c = (0 , l c ) representing a part of the plant cell cytoplasm, where l c denotes the length of the cell segment we consider. The boundaries of c are denoted by n modelling the cell nucleus, and c representing the cell membrane ( Fig. 7 ). We assume the diffusion of BR ( b ), BKI1 ( k ) and BZR-p ( z p ) in the cytoplasm: The only reaction that takes place in c is the degradation of BR since we assume that the phosphorylation status of BZR is modulated only in the nucleus. The dynamics occurring on the plasma membrane c are the interactions between b, k , and receptor BRI1 ( r k , r b ). Since we assume that the receptors are membrane-bound a system of ODEs, similar to the corresponding equations in system (1) , is considered to model the dynamics of r k and r b . The effect of the interactions between b, k, r k , and r b on the dynamics of b and k is defined by Robin boundary conditions for b and k on c . Finally, we assume that the BZR-p cannot diffuse out of the cell, which we model by a zero-flux boundary condition on  (7) -(9) . BR (red circles), BKI1 (blue squares) and BZR-p (yellow diamonds with black dots) diffuse freely in the cytoplasm, where BR is also degraded. At the membrane, both BR and BKI1 are perceived by BRI1 (black y-shapes) and form complexes with it. In the nucleus (brown) BKI1 activates dephosphorylation of BZR-p to BZR (yellow diamonds), and inhibits phosphorylation of BZR to BZR-p. BR is synthesised in the endoplasmic reticulum (ER, grey crescent) which is continuous with the nuclear membrane. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7.
A diagram of the one-dimensional domain in which system (7) -(9) was solved. c = (0 , l c ) represents the cytoplasm, c the plasma membrane, and n the nucleus.
c . Thus on c we have In the nucleus n we consider both the phosphorylation and dephosphorylation of BZR. Although the exact subcellular location of BR biosynthesis has not been experimentally demonstrated, the likely location is the endoplasmic reticulum ( Shimada et al., 2001 ). The endoplasmic reticulum is continuous with the nuclear membrane, hence we model the BR biosynthesis as occurring on n . We assume that BKI1 cannot enter the nucleus and consider zeroflux boundary conditions for k on n .
Since ˜ r k , ˜ r b and ˜ z are confined to the boundary, they have units of mol / m 2 . Thus we have the following scaled relationships between variables ˜ r k , ˜ r b and ˜ z in (7) -(9) and the corresponding variables r k , r b and z in (1) In order to preserve the balance of units some parameters from model (1) also had to be rescaled, specifically˜ α b = l c α b ,˜ δ z = l c δ z and ˜ θ b = θ b /l c . We estimated the diffusion constant D b = 60 μm 2 min −1 by taking the value reported for Progesterone (a steroidal hormone with a similar structure to BL) in physiological solution from Sieminska et al. (1997) . We also took D k = D z = 0 . 125 μm 2 min −1 for the diffusion constant for proteins from Sturrock et al. (2011) . l c was taken to be 7.43 μm from measurements of root cell sizes from unpublished data. All other parameters were taken to be the same as in Table 1 or Table 2 . Model (7) -(9) was solved numerically to analyse the changes of the model solutions due to the spatial heterogeneity of signalling processes. In the fitted parameter regimes, averaged solutions of the model (7) -(9) have similar behaviour to solutions of (2) ( Fig. 8 ), however a distinct spatial distribution of BR is characteristic for (7) -(9) , see Fig. 9 . In the oscillatory parameter regime discussed in Section 4.1 , the PDE-ODE model (7) -(9) was found to have increased amplitude and period of oscillations, see Fig. 8 . For solutions of model (7) -(9) we also observe oscillatory behaviour for a much wider range of values of δ z and ρ z than for model (2) .

Derivation of a mathematical model for crosstalk between the BR and GA signalling pathways
Here we consider a mathematical model for the crosstalk between the BR and GA signalling pathways. BRs and GAs can play similar independent roles in development both having important roles in growth ( Clouse and Sasse, 1998;Úbeda Tomás et al., 2009 ), as well as acting together via shared gene expression, and through interactions between their signalling pathways ( Bouquin et al., 2001;Tanaka et al., 2003 ). We aim to use our model to analyse the three different mechanisms of crosstalk between the BR and GA signalling pathways proposed in Li and He (2013) , Tong et al. (2014) and Unterholzner et al. (2015) and to attempt to establish which mechanism has the most significant influence on the dynamics of the pathways.

Mathematical modelling of GA signalling
A detailed model of the GA signalling pathway was derived and examined in Middleton et al. (2012) , consisting of a system of 21 ODEs and 42 parameters. It considers both the interaction of the GA 4 , GID1 and DELLA, and the GA 4 biosynthesis pathway, where GA 12 is converted to GA 15 is converted to GA 24 is converted to GA 9 , all catalysed by enzyme GA20ox, and GA 9 is converted to GA 4 , catalysed by GA3ox.
Since a model for crosstalk must include variables for both pathways involved, as well as potentially introducing new variables specific to their interactions, a system of ODEs describing  (6) , and the PDE-ODE model (7) -(9) for various parameter values. For comparison, the solutions to the PDE-ODE model have been averaged over the space, and the initial conditions are such that they are equal for the ODEs and the averaged PDEs. a) Both models were solved with the fitted parameters (parameters in Tables 1 and 2 correspond to set 1 and set 2 respectively), and show similar behaviour. b) Models were solved for δz = 1 . 02 × 10 −2 , ρz = 1 . 33 × 10 −3 , θz = 41 . 2 and all other parameters as in Tables 1 and 2 respectively. The ODE model tends quickly to steady state, whereas the PDE-ODE model exhibits damped oscillations. c) Models were solved for δz = ρz = 4 , θz = 41 . 2 , and all other parameters as in Table 1 . Both systems exhibit periodic solutions, but the PDE-ODE model has a much reduced frequency and much increased amplitude as compared to the ODE model. d) δz = 14 , ρz = 35 , θz = 41 . 2 , and all other parameters as in Table 1 , ODE model has moved outside of the oscillatory parameter regime, however the same is not true for the PDE-ODE model, which continues to have oscillatory behaviour for values of δz , ρz in excess of 10 0 0. Fig. 9. The spatial distribution of BR in the cell segment c , solution of (7) -(9) , at different times, where 'parameter set1' and 'parameter set2' correspond to parameter values in Table 1 and Table 2 respectively. The steady state concentration of BR is greater for solutions with parameter set 1, Table 1 , than for solutions with parameter set 2, Table 2 . such a model may be very large. This being the case, we first reduce the size of the GA signalling model that was derived in Middleton et al. (2012) . We reduce the model by assuming that the dynamics of the molecules involved in the GA biosynthesis pathway are much faster compared to the dynamics of other processes involved in the signalling pathway that their levels remain at a steady state. This enables us to write a system of algebraic equations governing GA 12 , GA 15 , GA 24 , GA 9 , GA20ox and GA3ox, as well as their relevant mRNAs and complexes which may then be solved such that they may therefore be substituted out of the main system. These substitutions generate a new term governing the biosynthesis of GA 4 , where biosynthesis is directly dependent on DELLA concentration, which is then simplified to a fourth order Hill function, signifying the 4 steps in the biosynthesis pathway. We also assume that only one configuration of DELLA.GID1 c .GA 4 may be formed. Overall this removes 27 parameters and 13 variables from the system in Middleton et al. (2012) , and introduces a new term for the GA 4 biosynthesis with new parameters α g and ϑ g . The interested reader may examine both the reduction process as well as explicit calculations of α g and ϑ g in Appendix D . Then for r , r o g , r c g , r d , r m and d m , the concentrations of GID1, GID1 o .GA 4 , GID1 c .GA 4 , DELLA.GID1 c .GA 4 , GID1 mRNA and DELLA mRNA respectively, we obtain dr dt = −β g rg + γ g r o g + α r r m − μ r r,  and for d l and g , the concentrations of DELLA and GA 4 respectively, we have All parameter values are taken from those for the full model in Middleton et al. (2012) ( Table 3 ). For all components present in both the full and the reduced model, the initial conditions are taken to be the same as in Middleton et al. (2012) . Both the full and reduced models were solved numerically, and these solutions are plotted in Fig. 10 . The two models show excellent agreement being almost indistinguishable for six terms, and only minor short-term discrepancies for DELLA.GID1 c .GA 4 and DELLA. Having thus successfully reduced the GA signalling pathway model while retaining its core behaviour, we derive a model of the crosstalk between the BR and GA signalling pathways by coupling models (2) and (10), (11) .

The crosstalk model
We consider direct crosstalk between the BR and GA signalling pathways. Besides shared gene expression, two main mechanisms for BR-GA crosstalk have been suggested: direct interaction at the level of transcription factors BZR and DELLA ( Li and He, 2013 ), and BZR-mediated biosynthesis of GA ( Unterholzner et al., 2015 ). Here we derive a mathematical model to compare the behaviour of the system under the three different mechanisms of crosstalk (interaction between BZR and DELLA, BZR-mediated biosynthesis of GA, and a combination of both of them) and to analyse their effects on the BR and GA signalling processes.
The existence of an interaction between BZR and DELLA constituting a crosstalk between the BR and GA signalling pathways was established in Bai et al. (2012) , Gallego-Bartolomé et al. (2012) and Li et al. (2012) . Further to this, in Li and He (2013) it was found that the formation of a complex BZR.DELLA inhibits the transcriptional activities of both BZR and DELLA ( Fig. 11 ).
A second mechanism of BR-GA crosstalk has been proposed where BZR also has a direct influence on the GA 4 biosynthesis pathway ( Fig. 12 ). In Tong et al. (2014) it was discovered that BR induces the expression of GA3ox-2, leading to the accumulation of GA 1 , the most bioactive GA for rice. It was also discovered that the external application of BR induces GA20ox expression in Arabidopsis ( Lilley et al., 2013 ). Later it was shown that BRs both regulate GA20ox expression, and are required for the transcription of GA3ox  ( Unterholzner et al., 2015 ). Despite good evidence for the coexistence of these two mechanisms, it is still unclear to what extent each mechanism operates to effect changes in BR and GA signalling processes ( Ross and Quittenden, 2016;Tong and Chu, 2016;Unterholzner et al., 2016 ).
To derive the model for BR-GA crosstalk we coupled model (2) and model (10), (11) by adding new interaction terms describing the crosstalk mechanisms. Two terms were added, a term corresponding to the interaction of BZR and DELLA, which required the introduction of a new variable denoting the concentration of the complex BZR.DELLA, and a term corresponding to the BZRmediated biosynthesis of GA. By varying the parameters of these new terms, we were then able to analyse the influence that each mechanism exert over the corresponding signalling processes.
BZR.DELLA complex formation is modelled as a reversible reaction between BZR and DELLA. The concentration of BZR.DELLA is denoted by z d , and the parameters β z and γ z denote the binding rate of BZR and DELLA and the dissociation rate of BZR.DELLA, respectively, From this we modify system (2) to take account of the dynamics of BZR.DELLA db dt System (11) is modified by taking account of the dynamics of BZR.DELLA, and extending the Hill function describing GA biosynthesis to include BZR-mediated gene expression. For the crosstalk model we do not assume that exogenous GA is entering the system, hence the term φ g (ω g − g) describing this process in (11) is removed. BZR is assumed to control gene expression in the same manner as DELLA, but with relative activity described by φ z , where φ z is the ratio between the binding thresholds of DELLA and BZR dd l dt = −β d d l r c The values for the parameters β z , γ z governing the direct interaction between DELLA and BZR are estimates from values for similar complex formation and separation from Middleton et al. (2012) . We first assumed that the binding thresholds of BZR and DELLA were equal i.e. φ z = 1 , and then analysed the influence of variation in φ z on the dynamics of the solutions of (10) , (12) -(14) .
In order to analyse the relative effects of the different crosstalk mechanisms, we examined the behaviour of the model for different values of β z , γ z and φ z . The model was solved numerically for four different conditions: no crosstalk ( β z , γ z , φ z = 0 ), crosstalk via BZR activated GA biosynthesis only ( β z , γ z = 0 ), crosstalk via BZR.DELLA complex formation only ( φ z = 0 ), and crosstalk via both mechanisms ( Fig. 13 ). Only for variables including GA were there any differences when including BZR-mediated biosynthesis of GA, and for all other variables there was close agreement between the cases of no crosstalk and crosstalk via biosynthesis only, and close agreement between the cases of crosstalk via complex formation only and crosstalk via both mechanisms. To examine the influences of both mechanisms on the BR and GA signalling pathways, the relevant parameters were varied. Variation in β z and γ z led to longterm behaviour changes for the components of the BR signalling pathway, and short-term changes for the components of the GA signalling pathway ( Fig. 14 ). Variations in φ z had very little effect with differences in behaviour in variables containing GA for the cases φ z = 0 and φ z > 0 (results omitted). From this we concluded that direct interaction between BZR and DELLA is the predominant form of crosstalk between the BR and GA signalling pathways.
In order to examine the effects of mutations in the BR signalling pathway upon BR-GA crosstalk, changes in parameters δ z , ρ z and θ z are introduced into model (10) , (12) -(14) , Fig. 15 . Damped oscillations in the dynamics of the components of the BR signalling pathway induced no such behaviour in the dynamics of the GA signalling pathway, and despite causing some changes in the early behaviour of the dynamics of the GA signalling pathway their long term behaviour remains similar. Alongside the induced perturbation in the BR signalling pathway, we also varied β z and γ z to examine how crosstalk affects the dynamics of both pathways in this case. We found in the case that β z / γ z is large, the dynamics of the BR signalling pathway are virtually unchanged from the cases with smaller values of β z / γ z , however the long-term behaviour of the GA signalling pathway was significantly effected, with all components except GA 4 and DELLA m having substantially different stationary solutions. Given the observed large effect of crosstalk on the BR pathway as compared with the GA pathway, see Figs. 13 and 14 , we analysed how directly altering each pathway would affect the other ( Fig. 16 ). We modelled the overexpression of BR and GA hormones by increasing parameters α b and α g . Overexpression of BR (increase of α b ) resulted in major changes in the BR signalling pathway with increases in the concentrations of BR, BKI1 and BZR, but almost negligible change in the GA signalling pathway. Overexpression of GA had more widespread effects, causing increases in the concentrations of BZR, GID1 o .GA 4 , GID1 c .GA 4 , DELLA.GID1 c .GA 4 and GA 4 , and decreases in the concentrations of BR, BKI1, GID1 o , GID1 m and DELLA.BZR.
In all cases of BR-GA interactions discussed here we considered the parameter values for the BR signalling pathway model as in Table 1 . However similar behaviours are observed also if considering the parameter values as in Table 2 , hence we do not present the simulation results for those cases.

Modelling spatial heterogeneity in the crosstalk signalling
We also considered the spatial heterogeneity in the crosstalk model. Here, the coupled PDE-ODE model for the BR signalling pathway (7) -(9) was modified to include the components of the GA signalling pathway. GA was assumed to be able to diffuse freely throughout the cell, and all other components of the GA signalling pathway were assumed to be nuclear-localized. We retain the same assumptions for the components of the BR signalling pathway as in  (14) ; ' no crosstalk ' refers to the case where βz = γz = φz = 0 ; ' biosynthesis ' refers to the case where the only crosstalk between the BR and GA signalling pathways is BZR-induced biosynthesis of GA ( βz = γz = 0 , φz = 0); ' complex ' refers to the case where the only crosstalk between the BR and GA signalling pathways is complex formation of BZR and DELLA ( βz , γ z = 0, φz = 0 ); ' full ' refers to the case where crosstalk between the BR and GA signalling pathways can be via both mechanisms BZR-mediated GA biosynthesis and complex formation of BZR and DELLA ( βz , γ z , φz = 0).

Fig. 14.
Variation of the values of βz and γ z suggests that the BR signalling pathway is more sensitive to the effects of crosstalk than the GA signalling pathway; ' no crosstalk ' corresponds to the case βz = γz = 0 , ' small ' the case where βz / γ z has order of magnitude zero, i.e. βz / γ z ≈ 0.75 , for ' standard ' βz / γ z has order of magnitude two, i.e. βz / γ z ≈ 75, and for ' large ' βz / γ z has order of magnitude four, i.e. βz / γ z ≈ 7500.
Section 5 . This led to a set of reaction-diffusion equations for BR, BKI1, BZR-p the same as in (7) , and for GA: We assume that receptor based interactions of BR, BKI1 and BRI1 remain the same as in (8) , and no influx of exogenous GA on the plasma membrane: The production of BR, change in phosphorylation status of BZR, and interactions between BZR and DELLA occur in the plant cell nucleus and are modelled as flux boundary conditions and ODEs on the boundary n : All processes of the GA signalling pathway apart from degradation of GA are localised to the nucleus and are modelled by a system of ODEs for GID1 o , GA 4 .GID1 o , GA 4 .GID1 c , DELLA.GA 4 .GID1 c , DELLA, GID1 m and DELLA m , and a flux boundary condition for GA on the boundary n : Fig. 15. Oscillations in the BR pathway lead to short term changes in the dynamics of the GA pathway, and increasing βz / γ z causes more long term effects on the GA signalling pathway. In the ' standard ' case all parameters are as in Tables 1 and 3 and βz = 10 μM −1 min −1 , γz = 0 . 133 min −1 . For the ' BR perturbed ' case parameters δz and ρz have had their orders of magnitude increased by three, θz = 41 . 2 μM −1 , and all other parameters are as in Tables 1 and 3 and βz = 10 , γz = 0 . 133 . The ' BR perturbed, βz / γ z large ' has the same parameter values as the ' BR perturbed ', however here βz / γ z has order of magnitude four, i.e. βz / γ z ≈ 7500, similar to Fig. 14 .   Fig. 16. Examining the effects of hormonal overexpression on the dynamics of solutions of model (10) , (12) -(14) for the crosstalk between the BR and GA signalling pathways.
For the ' standard ' case all parameters are as in Tables 1 and 3 and βz = 10 μM −1 min −1 , γz = 0 . 133 min −1 ; in the ' BR overexpression ' case, α b 's order of magnitude is increased by one; in the ' GA overexpression ' case αg 's order of magnitude is increased by one.
Finally crosstalk is modelled via formation of the BZR.DELLA complex dz d dt Similar to model (7) -(9) , to ensure appropriate units we obtain rescaled relations between boundary-localised variables in model (15) -(19) and the corresponding variables in model (10) , (12) - (14) , Since r m and d m were already relative quantities with no units there was no need for them to be scaled. To balance units certain parameters also had to be rescaled: parameters α b , δ z , α r , α d , α g , ϑ r and ϑ d from (10)  System (7), (8), (15) -(19) was solved numerically, and these solutions were averaged over the space and then compared to the solutions of system (10) , (12) -(14) . The behaviour of the two models was similar for most of the parameter sets considered in Section 6.2 , see e.g. standard cases in Figs. 15 and 17 , however spatially heterogeneous steady states are characteristic for model (7), (8), (15) -(19) , see Fig. 18 . Spatial heterogeneity does have a significant effect in the case where perturbations in the phosphorylation of BZR lead to damped oscillations in the BR signalling pathway, with much different behaviours of solutions of system (7), (8), (15) -(19) than for (10) , (12) -(14) , comparing Figs. 15 and 17 . When the BR signalling pathway is 'strongly perturbed', i.e. δ z and ρ z are increased sufficiently high, the oscillations in the components of the BR signalling pathway also cause oscillatory behaviour in the components of the GA signalling pathway, suggesting that spatial heterogeneity may contribute to instability of solutions in extreme cases.

Discussion and conclusion
Plants rely on complex integrated hormonal signalling pathways in order to respond and adapt to changing environmental conditions. The BR and GA signalling pathways have a diverse range of effects on plant growth and developmental processes, some of which overlap. To examine the behaviours of these signalling pathways and their interactions we developed mathematical models of the BR signalling pathway, and crosstalk between the BR and GA pathways under assumptions of both spatial homogeneity: models (1) and (10) , (12) -(14) , and spatial heterogeneity of signalling processes: models (7) -(9) and (7), (8), (15) - (19) .
The parameters in the model for the BR signalling pathway (2) were determined upon validating the model by comparing its numerical solutions to experimental data from Tanaka et al. (2005) . Using numerical optimisation techniques we were able to get a good fit for the diverse data sets corresponding to the different growth conditions considered in Tanaka et al. (2005) . Our calculations resulted in values for δ z and ρ z of a much lower order of magnitude than the rest of the parameters in the model, suggesting that the phosphorylation of BZR occurs on a much slower time scale than e.g. perception of BR by BRI1. The model parameters were optimised under two different constraints on the value of β k , (4) and (5) , corresponding to two different mechanism governing the interactions of BR, BRI1, and BKI1, where BR and BKI1 binding to and dissociation from BRI1 is instantaneous, or a complex BR.BRI1.BKI1 can be formed, respectively. Most parameters had very similar or identical values, with the only notable differences being for (naturally) β k , and the parameters governing the phosphorylation of BZR: δ z , ρ z . Overall our model fitted the data well for both mechanisms (4) and (5) , with R 2 values of 0.89 and 0.92 respectively, suggesting that both of these mechanisms may accurately predict the dynamics of the BR signalling pathway. The only growth condition for which the model had a less accurate fit was for when the growth medium was supplemented with BL. Qualitative analysis of system (6) revealed that it has only one steady state for the parameter space P , and that this steady state is stable in a neighbourhood of the parameter sets defined in Tables 1 and 2 . Bifurcation analysis showed that for a bounded set of values of ( δ z , ρ z ), with θ z = 41 . 2 μM −1 and all other parameters as in Table 1 , system (6) has periodic solutions, i.e. system (6) undergoes a Hopf bifurcation upon varying parameters δ z , ρ z , corresponding to the dephosphorylation rate and phosphorylation rate of BZR respectively. No such bifurcations were found when considering a large range of values for δ z , ρ z , and θ z = 41 . 2 μM −1 and all other parameters as in Table 2 . However, since the values of δ z , ρ z for which a bifurcation exists are likely to be biologically unrealistic as they are much greater than the values yielded by validation of the model by experimental data, this suggests that within some range of normal function the BR signalling pathway has good stability, which is of crucial importance for proper function of plant tissues. Results also suggest that if BR and BKI1 associate and dissociate from BRI1 instantaneously, as in mechanism (4) , stability of the pathway is principally dependent upon the phosphorylation state of BZR. Further if BR.BRI1.BKI1 can be formed, as in mechanism (5) , stability of solutions to model (6) for a wide range of all parameters ensures effective BR homeostasis.
Numerical solutions for the spatially heterogeneous model (7) -(9) averaged over space and ODE model (2) for the BR signalling pathway demonstrated similar behaviours for the parameter sets when the solutions of (7) -(9) and of (2) converge to the steady state as t → ∞ . In the oscillatory parameter regime discussed in Section 4.1 , the PDE-ODE model exhibited distinct behaviour to the ODE model ( Fig. 8 ), suggesting that under such conditions the spatial heterogeneity of the signalling pathway has a large influence on the dynamics of the molecules involved in the BR signalling pathway.
We investigated the effects of interactions between the BR and GA signalling pathways and the BR-GA signalling crosstalk mechanisms ( Li and He, 2013;Tong et al., 2014;Unterholzner et al., 2015 ) by coupling (2), (10) and (11) to obtain model (10) , (12) -(14) , in order to establish which mechanism is more biologically significant ( Ross and Quittenden, 2016;Tong and Chu, 2016;Unterholzner et al., 2016 ). We examined the effects of BZR-mediated biosynthesis of GA and of BZR.DELLA complex formation on the behaviour of the BR and GA signalling pathways. We modelled the cases where there was: no crosstalk, BZR-mediated biosynthesis of GA only, BZR.DELLA complex formation only, and both BZRmediated biosynthesis of GA and BZR.DELLA complex formation. The cases of no crosstalk and BZR-mediated biosynthesis exhibited similar behaviour to each other, and the cases of BZR.DELLA complex formation and both BZR-mediated biosynthesis of GA and   Tables 1 and 3 and βz = 10 μM −1 min −1 , γz = 0 . 133 min −1 . The concentration of BR starts of uniform, and is produced at x = 7 . 43 μm, the nucleus, leading to the greatest concentration at this point. The overall concentration increases over the first 8 h, but then decreases steadily, arriving at a spatially heterogeneous steady state after 120 h. The concentration of GA is also greater at x = 7 . 43 μm where it is produced, but GA attains its spatially heterogenous steady state after 24 h, much more quickly than BR.
BZR.DELLA complex formation exhibited similar behaviour to each other ( Fig. 13 ). Further, upon variation of the ratio between DELLA and BZR binding thresholds φ z the change in behaviour of most components of both pathways was negligible, whereas the variations in BZR.DELLA binding β z and dissociation γ z rates had a strong effect upon the BR signalling pathway ( Fig. 14 ). From this we concluded that BZR.DELLA complex formation has a greater effect on the behaviour of the BR and GA signalling pathways. This suggests that interactions between BZR and DELLA exert more influence over the dynamics of the pathways than BZR-mediated biosynthesis of GA and are more likely to be able to promote any effective change in the behaviour of the signalling processes. Numerical results also demonstrated that for all parameter sets where solutions tended to a steady state, the solutions of model (10) , (12) -(14) and the averaged over space solutions of model (7), (8), (15) -(19) exhibit similar dynamics, although spatially heterogenous steady-states are obtained for BR and GA concentrations ( Fig. 18 ).
To examine whether disturbance of one signalling pathway had a greater effect on the other, we modelled the effects of overexpression of both BR and GA. Overexpression of BR resulted in large changes in the dynamics of the components of the BR signalling pathway only, and overexpression of GA resulted in large changes in the dynamics of the components of the GA signalling pathway, and small changes in the dynamics of the components of the BR signalling pathway ( Fig. 16 ). In addition variation of β z and γ z led to large changes in the dynamics of the BR pathway but only short term changes in the dynamics of the GA pathway. From this we concluded that in general perturbations in the GA signalling pathway have greater influence on the BR signalling pathway than vice versa.
To examine the effects of mutations in the BR signalling pathway on the GA pathway, models (10) , (12) -(14) and (7), (8), (15) -(19) were solved numerically considering different values for the parameters δ z and ρ z governing BZR phosphorylation. For the components of the GA signalling pathway solutions to the ODE model showed different short-term behaviour but tended to the same steady state, whereas the components of the BR signalling pathway tended to different steady state. If varying β z and γ z such that β z / γ z was large in the case when we have damped oscillation in the components of the BR signalling pathway, the dynamics of some of the components of the GA signalling pathway were significantly altered but the dynamics of the components of the BR signalling pathway were unchanged ( Fig. 15 ). For the spatially heterogeneous model (7), (8), (15) -(19) , in the oscillatory parameter regime oscillations in the BR signalling pathway propagated into the GA signalling pathway ( Fig. 17 ). From this we conclude that only in the case where there is disturbance in both the BR signalling pathway and the mechanism of crosstalk between the BR and GA signalling pathways the dynamics of molecules in the BR signalling pathway have greater influence on the dynamics of molecules in the GA signalling pathway. This, together with the results discussed in the previous paragraph, suggests that under normal crosstalk the stability of individual signalling pathways is maintained even if there are disturbances in the other pathway.
To conclude, our analysis of new mathematical models for BR signalling pathway and the crosstalk between the BR and GA signalling pathways, derived here, provide a better understanding of the dynamics of signalling processes, dependent on model parameters and spatial heterogeneity. Our results for the BR signalling pathway highlight its stability, and suggest that this stability is particularly dependent on the mechanisms governing the phosphorylation state of BZR and the subcellular locations of these processes. Our results suggest that direct interaction between BZR and DELLA exerts a larger influence on the dynamics of the BR and GA signalling pathways than BZR-mediated biosynthesis of GA, and hence may be the primary mechanism of crosstalk between the two pathways. Our analysis indicates that during normal plant function, the GA signalling pathway exerts more influence over the BR signalling pathway than BR on GA, but mutations in the BR signalling pathway cause BR signalling to exert some short time influence over GA pathway and greater influence when coupled with disturbances in the crosstalk mechanism. Both BR signalling and crosstalk between BR and GA signalling are important for plant growth and development. Our modelling and analysis results also can be used to model the interactions between growth and signalling processes in order to better understand the influence of the BR and GA signalling pathways on growth and developmental processes in plants. Transformed parameters when fitting the non-dimensionalised model (A.8) against experimental data from Tanaka et al. (2005) , and using (3) and (4) . Due to the nature of the non-dimensionalisation, the dimensional parameters θ b and Z tot cannot be found.  Tanaka et al. (2005) , and using (3) and (5) . Due to the nature of the non-dimensionalisation, the dimensional parameters θ b and Z tot cannot be found.

Constant
which means that we may not only define β k and β b as but that we may further divide β k through by β b to show that , and hence A value of K m = 4 . 28 μM was reported in Wang et al. (2014a ), giving an expression β k = 2 . 62 × 10 −3 β b . The steady state solution of the second equation in (2) may thus be written as which may be solved for k * , and hence an expression for [ BKI 1] 0 is found In order to check whether a more accurate estimation of the parameters was given by fitting a non-dimensional model which does not scale time by one of the fitting variables: (A.8) was also fitted to the experimental data. The parameters whose dimensional values could be found had close agreement with those found fitting model (2) , Tables A.4 and A.5 , and so it is better to consider the fitting of the dimensional model (2) to the experimental data. To prove that M is positive invariant, we first consider M 1 := { (b, k, z) ∈ R 3 | b, k, z ≥ 0 } and show the non-negativity of solutions of (6) . Restricting b = 0 , we obtain that f · (1 , 0 , 0) = β k (κ − 1 + k ) k +

B2. Calculations for the proof of Theorem 2
Any steady state ( b * , k * , z * ) of (6) satisfies 0 = 1 Using simple algebraic manipulation on the system above, we find b * and z * in terms of k * b * = 1 Then k * is defined as a root of the following non-linear function which is positive for all p ∈ P and k * ∈ [0, 1].

B3. Characteristic equation
The Jacobian for system (6) evaluated at the steady state ( b * , k * , z * ) is given by The coefficients of the characteristic equation associated to Jacobian J (B.1) are defined as follows where ( b * , k * , z * ) is a steady state of the system (6) .

Appendix C. Non-dimensionalisation of model (7) -(9) of the BR signallng pathway
Since in the PDE-ODE model (7) -(9) some of the variables are defined on the boundaries we must consider them in units of mol / m 2 instead of mol / m 3 , compared to the ODE model ( 2 ). To adapt the units, we scale by the length of the cell segment, i.e. (D.9h) The BR signalling processes occurring in the nucleus The GA signalling processes occurring in the nucleus dr dt = −β g rg + γ g r o g + μ r (r m − r) (E.5)