Responses to auxin signals: an operating principle for dynamical sensitivity yet high resilience

Plants depend on the signalling of the phytohormone auxin for their development and for responding to environmental perturbations. The associated biomolecular signalling network involves a negative feedback on Aux/IAA proteins which mediate the influence of auxin (the signal) on the auxin response factor (ARF) transcription factors (the drivers of the response). To probe the role of this feedback, we consider alternative in silico signalling networks implementing different operating principles. By a comparative analysis, we find that the presence of a negative feedback allows the system to have a far larger sensitivity in its dynamical response to auxin and that this sensitivity does not prevent the system from being highly resilient. Given this insight, we build a new biomolecular signalling model for quantitatively describing such Aux/IAA and ARF responses.

The species and the dynamical equations of the minimal model were specified in the Main part of the paper (see Sect. 2). The interactions between the species in this model are based on simplified mass-action reactions. Specifically, the kinetic parameter α (respectively β) quantifies how the presence of auxin (respectively IAA) "degrades" IAA (respectively ARF) but this is an effective reaction which does not consume the active agent, here auxin (respectively IAA). This simplification, which enforces the feed-forward nature of the system, facilitates the mathematical analyses that we now present. Nevertheless, in the last sub-section of this part we shall show that our conclusions hold even if this simplification is dropped.

A. Steady states
In the steady state, the right hand sides of Eqs.1-3 are set to 0, leading to the following conditions on the steadystate concentrations: [auxin] ss = S auxin τ auxin , (S1) [IAA] ss = S IAA /(α [auxin] ss + 1/τ IAA ), [ARF ] ss = S ARF /(β [IAA] ss + 1/τ ARF ), where S auxin , S IAA and S ARF are necessarily time-independent and the "ss" subscript on the concentration of each species denotes that it is taken in the steady-state. In the last equation, one may substitute [IAA] ss by its expression in terms of S auxin or [auxin] ss if desired. If S IAA is fixed (unregulated case), all the equations are explicit, showing that there is a unique steady-state solution. In the presence of regulation, S IAA is not a priori known and must also be determined via the steady-state conditions.

B. Ensuring identical steady states for the regulated and unregulated minimal models
In the auxin signaling networks found in various organisms, there is a regulatory feedback [1], i.e., S IAA depends on [IAA], [ARF ] and even on other molecular species. A hypothetical justification could be that feedback allows for multi-stationarity, so Eqs. S1,S2,S3 could have more than one steady-state level of IAA and ARF for a given value of S auxin . However that would require a positive feedback, whereas the feedbacks seen in plants correspond to negative retroactions, i.e., when IAA is increased it lowers the value of S IAA . As a result IAA is considered to act as a self-inhibitor. As mentioned in the main part of the article, even in the absence of regulation, it is possible to push the saturation curves in the steady state to large values of the input flux S auxin . This can be seen by considering Eqs. S1-S3 in the unregulated model: IAA will become small (compared to its maximum value occurring when S auxin = 0) when the incoming auxin flux becomes much larger than S * auxin = 1/(ατ auxin τ IAA ). By decreasing α, the range in S auxin over which the steady state values of [IAA] and [ARF ] avoid saturation can be increased at will, no regulatory feedback is necessary. But there is a drawback: any dynamical response will be weak whenever the coupling (α) is small.
What form of regulation will ensure both large α and that the steady-state curves saturate only for large S auxin ? To make the comparison between the two models (without and with regulation) as fair as possible, we shall impose their input-output relation curves (cf. Figures 1:B1-B3) to be the exactly the same. Consequently, they will have identical static properties, and in particular the static response functions will coincide for the two models. To mathematically define the model with regulation, we take its differential equations to be those of the model without regulation but introduce two exceptions. First, α in the case of regulation has a larger value, α (reg) α (no−reg) . Second, the source term S IAA for production of IAA is modified to compensate that change in α, enforcing the whole steady-state input-output relation to be the same in the regulated and unregulated models. Such a modification can be thought of as introducing a regulation in the production of IAA, for instance at the transcriptional level. There is still a lot of freedom for how to make such a regulatory change depend on different molecular species while maintaining exactly the same steady-state values. For both biological and mathematical reasons, we shall take S IAA to be a function of only IAA, making IAA a (negative) self-regulator. Comparing the equations with and without regulation in the steady state, we set the change in IAA production to be equal to the change in IAA consumption, i.e., This form indicates that at low [IAA] (high auxin-influx rate) the IAA production rate in this regulated model is enhanced by a factor α (reg) /α (no−reg) compared to the unregulated case. Inversely, when auxin influx vanishes, the production rate of IAA is the same in the regulated and unregulated cases; this could have been anticipated since without such influx, α plays no role so one must then have S . In addition, it is easy to see that because the production rate decreases as [IAA] increases (as expected since IAA acts as a self-inhibitor), there is a unique solution to Eq. S2, and so the regulated system does not exhibit any multi-stationarity either. Lastly, plugging Eq. S5 into Eq. 2.2, we see that the regulated model is obtained from the unregulated model solely by multiplying the expression for d[IAA]/dt by the overall factor α (reg) /α (no−reg) .

C. Stability of the steady states and effect of regulation
For a given value of S auxin , the (unique) steady-state concentrations are determined by Eqs. S1-S3. Let In the absence of regulation, this Jacobian J is: In these expressions, all concentrations are taken at their steady-state values. In the presence of regulation, the only change is that the second line is multiplied by the factor α (reg) /α (no−reg) : Because the Jacobian is tridiagonal, its eigenvalues are given by the diagonal entries. These are all negative, showing that the steady state is always linearly stable. Furthermore, the steady state is all the more stable that these eigenvalues are large in absolute value. (Note that the associated characteristic times, generally referred to as the relaxation times, are given by minus the inverse of these eigenvalues. ) We then see from the second eigenvalue that regulation with α (reg) > α (no−reg) leads to enhanced stability.

D. Regulation's role in the linear dynamical response functions
The linearized dynamics determine not only the time-dependent behaviour of infinitesimal deviations from the steady state but also their response to infinitesimal perturbations of the input. Using the same notation as before, suppose that the auxin source term is allowed to have infinitesimal variations in time: S auxin (t) = S auxin + ∆S auxin (t). (S9) The system will respond to these variations in the input via where for what follows we shall take −→ ∆S(t) = {∆S auxin (t), 0, 0} i.e., we focus on perturbations corresponding to injecting auxin into the system according to an arbitrary time-dependent flux ∆S auxin (t). Assuming the system is in its steady state before the perturbation is applied, these equations lead to: where the 3 × 3 matrix for t ≥ t is called the linear dynamical response function. (By definition, χ(t − t ) is taken to be 0 for t < t .) Taking t ≥ t , the (ij)'th entry of this matrix, χ ij (t − t ), can be interpreted as the value of − − → ∆C i (t) arising if one applies to the system in its steady state a delta function pulse for the source −→ ∆S j at time t . This dynamical response function (a matrix function of time) is to be distinguished from the static response function R which is a time-independent vector whose components R 1 , R 2 , R 3 are simply the derivatives of the steady-state concentrations with respect to the intensity of the steady-state input: It is not difficult to show that: In Figures 1:C2-C3 we display the components χ i1 (t − t ) (i = 1, 2, 3 and t = 0) hereafter denoted χ auxin (t), χ IAA (t) and χ ARF (t), for the models with and without regulation. The case i = 1 provides the time-dependent response of auxin (which is not affected by regulation in this model). Examining Eq. 2.1, we see that ∆[auxin] rises instantaneously by one unit when a perturbation is applied in the form of an instantaneous (delta function at t = 0) pulse of unit integral, and thereafter it simply decays exponentially at rate −J 11 = 1/τ auxin . If one integrates over time to have what is called the total response to the pulse, we see that the total response is given by the ratio: total (S15) As a result, we see that whether there is regulation or not, the integral over time of the dynamical response function is the same. (Recall that regulation rescales J 21 and J 22 by the same factor α (reg) /α (no−reg) , see expressions S7 and S8). Interestingly, the generality of this result follows from the fact that this total response is given by Eq. S14 which itself depends only on the steady-state concentrations (cf. Eq. S13). Since by construction the input-output curves are the same with and without regulation, we see that the total response is also the same. The case of the response displayed by ARF (the third molecular species) can be treated similarly. After an influx pulse of auxin, ∆[ARF ] will initially rise quadratically in time and then it will reach a maximum. If the relaxation time of ARF is large compared to the time during which the IAA signal (the perturbation of IAA concentration) is significant, then this maximum value is approximately −β[ARF ] ss times the total response of ∆[IAA], itself given by Eq. S15. As a result, this maximum will be insensitive to regulation. If on the contrary the relaxation time of ARF is comparable to that of IAA or smaller, then the maximum value with regulation will be amplified approximately by a factor α (reg) /α (no−reg) , just as it is for IAA. After reaching its maximum, ∆[ARF ] will decay back to 0. And just as for ∆[IAA], the total ARF response is the same whether there is regulation or not.

E. Modifying the minimal model by using true mass action kinetics
In the minimal model considered so far, auxin degrades IAA but this degradation has no consequence on auxin. In the same vein, that model lets IAA degrade ARF but without the IAA dynamics being affected. In reality, the actions arise through the formation of complexes which lead to changing the dynamics of all the molecular actors involved. It is thus natural to ask whether the conclusions reached in that minimal (and feed-forward) model still hold if one uses a more realistic description of the dynamics that correctly includes mass-action in the reactions. That is the purpose of this section. For this extended model to be specified, we have auxin act on IAA by first forming an effective auxin-IAA complex which mimics the auxin-TIR1-IAA complex. This effective auxin-IAA complex then degrades IAA and releases auxin, in direct analogy with what happens biochemically. Similarly, to mimic the role of IAA on ARF, we have them form heterodimers according to mass action kinetics. These heterodimers act in effect to sequester ARF, just as in the non-minimal models. These changes make a bit more complex the minimal model but for the most part one can nevertheless study this modified model analytically.
The equations of the minimal model have to be modified to take into account the formation of the auxin-IAA and ARF-IAA complexes. Furthermore, in the spirit of our new model, we consider that ARF is long lived so that its total concentration (free plus in the ARF-IAA complex) is fixed. The use of mass action for the formation of complexes then changes Eqs. 1 -3: where we have also taken into account the choice of no production or degradation of ARF. The new parameters γ and δ are the dissociation rates of the auxin-IAA and ARF-IAA complexes. The dynamics of the concentration of auxin-IAA and ARF-IAA complexes adds one additional differential equation: and also one constraint: where ARF T is the total concentration of ARF (free or in the ARF-IAA complex). Because of this constraint, we do not need to follow explicitly the dynamics of [ARF-IAA] in this extended model as it is not an independent quantity. The steady-state concentrations of the four independent molecular species are readily obtained: [auxin] ss = S auxin τ auxin , (S21) Note that each equation involves only terms determined from previous equations so the whole system can be solved by using these equations in the order of appearance. So far we have implicitly considered no regulation. To include regulation, just as in the feed-forward minimal model, we rescale α and adjust S IAA so that the steady states remain the same. The procedure used in the feed-forward model works exactly in the same way, and in fact the expressions for S IAA in the two minimal models are identical and given in Eq. S4.
Having now defined the sequestrating minimal model with and without regulation, let us consider the steady-state behaviour for increasing S auxin . Supplementary Figures S1:A1-A3 shows that the inclusion of the mass action dynamics (and associated sequestration) does not qualitatively change the dependence of concentrations of auxin, IAA or ARF on S auxin (compare these curves to the ones in Figure 1:B1-B3).
What about the linear dynamical response function? To investigate that, we compute the Jacobian matrix giving the linearized dynamics about the steady state. In contrast to the feed-forward minimal model where the Jacobian was a 3 × 3 matrix, here there are 4 molecular species to consider as dynamical variables so the Jacobian is a 4 × 4 matrix. Except for that greater complexity, the framework for computing the response functions is the same. In Supplementary Figures S1:B1-B3 we show these functions giving the time dependence of the variations in concentrations of auxin, IAA and ARF when one introduces an infinitesimal pulse of auxin into the system at t = 0.
The conclusion is that again the minimal model with full mass action behaves very much like the feed-forward minimal model, the negative feedback showing an amplification in the dynamical response compared to the unregulated case.
From a biological point of view, note that the release from sequestration allows for a fast response to an auxin signal, bypassing any (long) times associated with transcription or translation. This feature is shared with the Vernoux model and also our new model of Sect. 4 in the Main text.

II. MATHEMATICAL ANALYSIS OF THE MODEL OF VERNOUX ET AL.
In the last sub-section, we show how the dynamics of the model provide a decomposition of the network into modules, without the need for any other information.

A. Rate of transcription
The rate equation of transcription for producing the IAA mRNAs (Eq. 3.5 in the main text) follows from the thermodynamical framework of Bintu et al. [3]. It involves multiple rates of transcription because Vernoux et al. assume that there are two AuxREs (auxin response elements) in the regulatory region of the gene coding for IAA. They allow for a first (basal) level of transcription when that region is free, another one when it is occupied by one ARF, and a last one when it is occupied by two ARFs. The binding of the heterodimer ARF-IAA to either of the AuxREs is assumed to shut off transcription, in line with the idea that this heterodimer acts as an inhibitor of transcription. The final expression used by Vernoux et al. for Eq. 3.5 is: The term κ − A is motivated by the presence of additional species competing for the binding to the regulatory region but not leading to any transcription. These could be for instance "repressor" ARFs which bind to the AuxRE but do not recruit the RNA polymerase (ARF repressors are known to exist). Furthermore, f and f A quantify the transcriptional amplification due to, respectively, one ARF activator and two ARF activators being bound to the regulatory region. The coefficients ω A , ω I and ω D indicate cooperativity effects stemming from the binding to the DNA of two ARF activators (ω A ) or the formation of dimers (ω I and ω D ); K d and B d are the dissociation constants for the ARF-IAA dimerization and the ARF binding to DNA reaction.

B. Steady States
To solve for the steady state(s), we set to 0 the left hand side of Eqs. (4)- (8) and (10) of the main text. A convenient strategy to solve the resulting equations is to first express all quantities in terms of the concentration of IAA protein. Based on Eqs. (6) and (7), at steady state one has: and It is useful to rewrite these relations as: . Then, by plugging Eq. S29 into Eq. 3.2, one obtains at steady state: Thus the steady-state concentrations of ARF, ARF-IAA, and IAA 2 are all known in terms of [IAA] ss . As a result, the rate of production of IAA messenger RNA is also known: .
(S31) Eq. 3.5 then provides the concentration of IAA messenger RNA. As a result, knowing the concentration of IAA protein determines the steady-state concentrations of all the other species, as promised.
To obtain now the steady-state concentration of IAA protein, we enforce equality of IAA production rate (given by π I [IAA m ] ss where [IAA m ] ss is known in terms of [IAA] ss ) and IAA total decay rate. This total decay rate is the sum of the degradation rates of IAA in all of its different forms. Using Eq. 3.1 and the previously derived equations, this total decay rate is: where again all quantities are to be re-expressed in terms of [IAA] ss . The associated self-consistent equation can be described graphically by having [IAA] ss on the x-axis and by putting on the y-axis (i) the rate of production of IAA protein and (ii) the total decay rate of IAA protein. Where these two curves cross gives the value of [IAA] in the steady state. The solution is unique in the Vernoux model because the first curve is strictly decreasing (the feedback loop is negative) while the second curve is strictly increasing (the more IAA is present, the more it gets degraded).

C. The dynamics provide an a posteriori decomposition of the network into modules
To define their dynamical equations, the authors of the Vernoux model [2] used the quasi-steady-state approximation to replace the ubiquitination pathway by effective rates for the degradation of IAA (cf. in particular Eqs. 4 and ??). Might it be possible to use this kind of approximation to further reduce their model? The linearized dynamics provides a systematic way to do so since each eigenmode of the Jacobian has a characteristic relaxation time. If we focus on the fastest relaxing mode in the Vernoux model, we find that it mainly involves ARF and ARF-IAA. This means that the kinetic coefficients for the formation and dissociation of ARF-IAA are large so one may use the quasi-equilibrium approximation for that reaction. One could thus replace the Vernoux model by a simpler one where ARF and ARF-IAA are no longer dynamical quantities but are given by Eqs. S30 and S27. To reduce still further the complexity of the model, we can consider the next mode, associated with the second shortest relaxation time. This eigenmode involves mainly IAA and IAA 2 . Thus using the same approach one may decide to have IAA 2 be in quasi-equilibrium with IAA and remove it as a dynamical variable. This leads us to ask what is the role of IAA 2 in the Vernoux model in the next sub-section.

D. Role of IAA homodimer
A sensible guess for the role of the IAA homodimer is that it may act as a kind of reservoir to buffer or delay the response to changes in auxin concentrations. To investigate this question, we have taken the Vernoux model and modified it so that there is no longer any formation of IAA 2 (this can be done for instance by simply setting k II = 0). Just as when we compared the Vernoux model with and without regulation, it is necessary here to adjust some of the rates for the steady states to be nearly the same with and without IAA homodimer. It is possible to do so as shown in the Supplementary Figure S2. Given that property of the steady-states, we can of course look at the dynamical properties. Supplementary Figure S2 shows that one consequence of removing the formation of IAA 2 is to quicken the response to an auxin perturbation. A second consequence is to increase the amplitude of the associated dynamical response. Both of these features seem desirable in a signaling network, justifying why in our new model, we have opted for no homodimerization of IAA.

E. Resilience in the Vernoux model
Characteristic resilience times in the Vernoux model were obtained by applying specific perturbations to the system and seeing when the amplitudes of the responses were down by a given factor. We did these measurements in the Vernoux model, both with and without the negative feedback. Illustrative results are given in Supplemental Table S1.

III. QUALITATIVE AND QUANTITATIVE ASPECTS OF OUR NEW CALIBRATED MODEL OF AUXIN SIGNALING
As mentioned in the Main part of the paper, the model is available on the BioModels repository. The first subsection explains the bioinformatic search for AuxREs motifs in Arabidopsis IAA regulatory regions. The second and third sub-section cover the mathematical aspects of our model. The forth provides details on the calibration of the resulting model while the actual numerical values of the 17 parameters are given in Supplemental Table S3.
A. Scanning the Arabidopsis IAA regulatory regions for AuxREs and analysis of their clustering Auxin Response Factors are transcription factors and as such they regulate other genes. Interestingly, they also are involved in IAA transcriptional regulation. In Arabidopsis, 23 different ARFs have been identified in the past few decades among which 5 seem to act as activators of IAA transcription (ARF5-8 and 19) while others are believed to act as repressors [4]. As previously discussed, ARFs can directly interact with IAA to form a hetero-dimer; that hetero-dimer probably competes with ARF in binding to AuxRE, and may act as an inhibitor of transcription. The putative static network underlying ARF and IAA interactions is given in ref. [2].
In this system, activators normally initiate transcription by binding to the Auxin Response Elements (AuxRE) upstream of the region coding for IAA. ARFs generally consist of a DNA Binding Domain (DBD) that recognises the consensus sequence TGTCTC along the DNA [5]. This binding motif has been further refined during the past years by both directed mutation experiments [6] and bioinformatic searches [7], leading to the construction of specific Position Weight Matrices (PWMs). ARF23 is the only exception and it is made of a truncated DBD [4]. All the other ARFs show furthermore a domain associated with their function as activators or repressors of transcription [8]. In ARF 1-12, 14-15 and 18-22, that domain is connected to two other ones, called III and IV, which allow binding to IAA and some other ARFs [4]; those two domains are the key domains for ARFs to form hetero-and homo-dimers. As we previously discussed, IAA binding ARFs prevents transcriptional activation. Activator ARFs can initiate transcription and are expected to bind DNA as monomers or dimers [6]. Binding by the dimer should be preferred to binding by the monomer if the AuxRE allows it [6]. However, a scan of the upstream regions of the 21 Arabidopsis IAA genes, searching for the AuxRE consensus sequences, did not show the existence of any pairs of positions allowing for binding via ARF homo-dimerization, i.e., consensus motifs were not found to be separated by about 10 nucleotides (one helix turn), other than in one case (Table S2). A refinement of such an analysis by using PWMs instead of the consensus sequence overcomes the constraint of requiring the presence of the exact consensus sequence. This can be done, e.g., via available open source programs which, given a PWM, allow one to score sequences and return a p-value for each position detected as significant [9]. In our specific case, we wrote a shell script to scan PWMs against the 21 Arabidopsis IAA upstream sequences. For completeness, we compared the obtained scores with those generated by scanning the same PWMs against reshuffled sequences: the program selects then as "true AuxRE" all those positions whose scores were higher than the highest random one. This allows one to get a set of positions labeled as AuxREs and potentially corresponding to true binding sites. In order to probe the possibility of binding by homo-dimers, we determined whether the distance between two subsequent motifs is close to 10 nucleotides (corresponding to a full turn of the DNA). Motivated by what is seen from the crystal structures [6], one can also expect the homodimers to have mirror symmetry in which case the paired AuxREs should be oppositely oriented (Figure 3:B). In our analysis of the sequences, we used two different PWMs, respectively from [6] and [7], from which we identified 16 and 38 hits (Figures 3:C1-C2 and C3-C4 respectively). A study of the distances between hits did not show any signal indicating presence of adjacent binding sites and thus we found no evidence in favour of homodimerization. This is quite in agreement with the results previously found by [10]. Such a divergence between that work and the indications from crystal structures leads to a puzzle in ARFs homo-dimer vs monomer DNA binding that might be sorted out with further experiments in the next few years. For the moment, in our model, we shall take into account both choices, where binding arises via a monomer or a dimer of ARF.
It remains an open question still which is the way ARFs bind to the AuxREs of the regulatory regions of IAA genes. In this regard, we shall consider in our modeling both possibilities, having transcription activation via ARF monomers or dimers.
Along with homo-dimerization, it has been recently proposed that ARFs may be able to form higher order complexes, e.g. in the form of oligomers [11]. In this scenario, both activators and repressors might participate together in binding the DNA; it would be interesting to study which are the consequences of oligomerization on the functioning of the system, either focusing on the ARF and IAA species or more generally within a full network including auxin signaling components.

B. The dynamical equations for each molecular species
The dynamical equations of our model are as follows. IAA transcription is described via the reaction: where F 1 is a function of [ARF ] and [ARF − IAA] (the two species which can bind to the AuxRE) and λ 1 is the rate of IAA mRNA production when the AuxRE is bound by an ARF (monomer) molecule. F 1 is in fact the probability that ahe AuxRE is bound by an ARF transcription factor; its functional form reflects an underlying thermodynamic equilibrium. Explicitly, we take the choice made in [12] and adapt it to our setting so that: Transcription requires that the AuxRE be bound by ARF. Note that ARF-IAA acts as a competitive binder and that when it is bound there is no transcription. These IAA messenger RNAs have a finite lifetime: This reaction proceeds at the rate µ IAAm times the concentration of the IAA messengers. The messengers are also used as templates for translation, leading to the production of IAA proteins. For the purpose of the modeling, the multiple steps involved in translation are simply coarse grained into one bulk reaction, i.e.: with a rate δ times the concentration of the IAA messenger. IAA degradation can arise via two routes: spontaneous, corresponding to a natural lifetime of the protein, or active, catalyzed by auxin-dependent biochemical processes. The spontaneous route is relevant mainly for very low concentrations of auxin and the corresponding rate parameter is µ IAA : The active degradation process depends on the formation of complexes containing auxin and its receptor TIR1. Auxin can bind to TIR1 proteins to form an auxin-TIR1 complex, which in turn may bind IAA. In this complex, IAA is then ubiquitinated, that is some of the protein's residues become tagged, changing IAA into IAA * . The set of these reactions reads as in ref. [12]: In the spirit of keeping the model relatively simple, we use mass action to describe all the kinetics of formation and dissociation of complexes. The other complexes in our model are the ARF-IAA and ARF 2 dimers, corresponding to the processes: where p a and p d are respectively the association and the dissociation rates of this heterodimer. Similarly, we have: where q a and q d are respectively the association and the dissociation rates of the ARF homodimer.
As mentioned in the previous section, the total concentration of ARF is considered to be fixed to some given Lastly, we assume that auxin is pumped into the system at a rate S auxin : Since auxin is introduced into the system, it must also be either evacuated or degraded. We choose the latter option, introducing an effective lifetime (independent of the cell's state): where µ auxin = 1/τ auxin . The union of all these processes allows us now to completely specify the model mathematically based on a set of ordinary differential equations for the time dependence of the concentration of all 10 molecular species: Our goal is to study how the concentration of IAA and of the driving factor ARF are affected by changes in auxin influx. Having specified the model in a formal way, it still has to be calibrated, i.e., its 17 parameters must be set, 15 being associated with rates, the two remaining being respectively the total amount of ARF and the total amount of TIR1.

C. Determining the steady states
We follow the strategy previously used to compute the steady states in the simpler models: we first express all concentrations in terms of that of IAA and then we derive the self-consistent equation for [IAA] ss . Although the expressions of all the different quantities in terms of [IAA] ss are complicated, at the end one obtains one selfconsistent equation for [IAA] ss which is easily solved numerically. To illustrate the successive steps, begin by imposing the steady-state condition on the complexes involving ARF. This leads to: where P = p a /p d and Q = q a /q d . By using the conservation law for total amount of ARF protein, one gets: where we use the notation of ref. [12]: K = k a /k d and L = l a /(l d + l m ). Note that we also used the steady-state expression for auxin concentration, [auxin] ss = S auxin /µ aux . As TIR1 total concentration is fixed too, one can use the same logic as before, leading to: whose solution for TIR1 can be expressed in terms of the parameters and [IAA] ss . We will denote this relation by Now that the concentrations of all species have been determined in terms of that of IAA, the self-consistent equation for [IAA] ss is obtained by imposing that the rate of production of IAA protein (translation) is equal to its (passive and active) decay rate. This gives immediately: The graphical representation of this self-consistent equation is provided in Figure 4:A. Just as for the Vernoux model, by monotonicity of the two curves there exist a unique solution for the steady-state concentration of IAA. We can compute it numerically for any given set of parameters and auxin influx rate S auxin . Then, by plugging the solution obtained for [IAA] ss into the other equations, one obtains all steady-state concentrations in the network (Figures 4:B1-B3). The behaviors of auxin, IAA and ARF concentrations as a function of S auxin are qualitatively the same as in the models analyzed in the previous sections.

D. Model calibration and use of diffusion limited reactions
The values of all 17 parameters of our new calibrated model are provided in Table S3. They are also given in the model definition in the BioModels repository. The setting of these values rested on estimates from published works such as [13,14], constraints coming from measurements of equilibrium constants, and theoretical bounds for "on" rates. This last type of constraint is not commonly used so we now explain in detail how it arises and can be used.
The rate k on in any reaction cannot be arbitrarily large because it is necessarily bounded from above by the so called "diffusion limit". To understand where this limit comes from, consider a reaction A + B which produces C. Motivated by the case where B is an enzyme, for pedagogical purposes consider that B is at a given point and that the A molecules diffuse throughout the cell volume. A necessary condition for forming a C molecule is the encounter or "collision" of an A molecule with B. The number of collisions per unit time felt by B, k S , can be computed [15] and it is referred to as the Smoluchowski encounter rate: where D A is the diffusion constant of the A molecules, R t is the radius of the target zone "offered" by B for the reaction and ρ A is the density of A molecules in the cell volume. Note that k S does not have the same dimensions as k on : indeed k S gives the number of collision per unit of time while k on has dimensions 1/(M · s). To relate k S to k on , one must first use the definition of k on : and then convert from numbers of molecules per unit volume to concentrations in moles per liter, giving where n A is the number of molecules of A in the cell, N A is Avogadro's number and V cell is the cell volume expressed in litres. Analogously, since we took one molecule of B in the cell, one obtains [B] = 1/N A V cell . Plugging these expressions into the differential equation for [C] and using k S one finally gets: The Einstein-Smoluchowski relation and Stokes' equation for mobilities of spheres in viscous liquids allow one to give an estimate of the diffusion constant, D A = k B T /6πηR A , where η is the viscosity of the liquid and R A is the radius of the molecule. Inserting this relation in the previous one for k on one gets: Recall that R t is the radius of the contact region offered by B for the reaction. For most reactions it is believed to be typically R t 1 − 2 nm. For our derivation, we assumed B did not diffuse. If instead both A and B molecules diffuse, the relation is modified by considering as the diffusion constant the sum of both diffusion constants, D A and D B . This leads to a final expression for k on given by: For our purposes, we considered T = 298.15K (as in [16]), η = 1.5 · 10 −2 P for the viscosity of the cytoplasm [17], and the other known constants also being given by the literature. To determine the characteristic sizes R A and R B when A and B are ARF and IAA, we first studied their domains using published crystal structures, in particular the DBD and the III/IV domains for ARFs and the III/IV domains for IAA. We visualized the structures of these domains using UCSF Chimera [18] to obtain the characteristic size of these molecules. Since each protein is multidomain, the estimates of R ARF and R IAA require some care. For ARFs, the crystal structures do not furnish any information about the junction domain between the DBD and the III/IV. We thus assumed this domain to have the same dimensions as the smallest one, i.e., III/IV. Our approach was to consider the multi-domain molecules as rods formed by aligning the domains, and to use the result for diffusion of rigid rods obtained computationally in ref. [19]. This approach leads to an effective radius of the protein using Stokes' law, leading to R ARF = 8.75 nm. In the case of IAA, we assumed the four domains to be all similar to the III/IV ones. The reasoning then led to an equivalent radius for IAA equal to R IAA 3.2nm. Plugging these values and R t = 2nm into the equation for k on and using the measurements of the dissociation constants, K D , in ref. [16] allowed us to determine k of f . The corresponding estimates of q a , q d , p a and p d are given in Tab. S3. Setting the other parameters in our new model described in Sect. 4 was possible using estimates from the literature (see in particular ref. [13]) so for instance the value of λ 1 was typical of genes constitutively expressed. Furthermore we used theoretical reasonings, for instance based on the bounds from the diffusion limit. With all this information it was possible to obtain values for all 17 parameters of this new model.

E. Complexity reduction via time-scale separation
In this subsection, we ask whether a decomposition of our new model into modules can be obtained automatically without any a priori biases. To do so we use the separation of time scales approach along with Principal Component Analysis [20]. Principal Component Analysis is generally used on sets of points in high dimensional spaces to obtain a lower dimensional projection. In effect, it first builds the multi-Gaussian distribution having the same mean and covariance matrix as the considered set of points and then it projects the high-dimensional set of points onto the main principal axes, that is the directions corresponding to the leading eigenvectors of that covariance matrix so as to maintain the maximum amount of variance. In the method of time-scale separation of a dynamical system, one can analyze the Jacobian matrix J describing the linearized dynamics about a steady state. When treating the dynamics of the deviations from the steady-state concentrations, the most general form of these linearized dynamics can be written as: where − −−− → ∆C(t) is the variation with respect to the steady-state of the different concentrations in the system. Whenever the Jacobian matrix is invertible, the previous equation can be integrated in terms of its eigenvalues and eigenvectors, namely λ k and |w k , where k goes from 1 to the total number of independent molecular species. In this space of eigenvectors the Jacobian is diagonal and so each eigenvector simply decays with time as exp(λ k t).
In general λ k can be complex, but for the decay rate of that eigenvector what matters is the real part of λ k (by stability of our system, this value is negative). One then defines the relaxation time of that eigenvector as the inverse of the opposite of that real part. The higher the rate of relaxation, the shorter the corresponding relaxation time. In the separation of time-scales approach, one wants to identify the fast parts of the dynamics. Formally, the quasisteady state approximation can be obtained by taking to infinity the eigenvalues of the fast eigenvectors. If one does that, the amplitude of those fast eigenmodes are in effect set to 0, corresponding to enforcing a (quasi-equilibrium) constraint between the involved molecular species. To interpret the meaning of those fast modes, one can use Principal Component Analysis. For illustration, consider the two fastest modes of our model which was defined in Sect. 4. We represent in Supplementary Figure S3, for each species, its components on those two eigenvectors based on the steady state when S auxin = 0.02 nM min −1 . The separation of time scales is seen to be quite simple: (1) the fastest mode (x-axis of the figure) consists mainly of auxin-TIR1, IAA and auxin-TIR1-IAA, meaning again that the association and dissociation rates of the triple complex are high so one has quasi-equilibrium between those species; (2) the second fastest mode (y-axis of the figure) consists mainly of ARF, IAA and ARF-IAA, meaning that the association and dissociation rates are high so one has quasi-equilibrium between those species. Were one to take the on and off rates of these two reactions to infinity, then formally the species ARF-IAA and auxin-TIR1-IAA could be removed from the dynamics and their dynamical equation replaced by a constraint. (Note that in practice such a change is not particularly helpful from a numerical point of view.) In Supplementary Figure S3 the "slow" species occur near the origin, e.g., IAA m and ARF 2 . The first is easy to justify: the relaxation time of IAA m is intrinsically long because the lifetime of the messenger RNAs are on the order of an hour, which is much longer than the time scale of any enzymatic biochemical reaction. The slowness of ARF 2 on the other hand was not necessarily expected. Indeed, since the association and dissociation rates of the heterodimer ARF-IAA are high, one might have anticipated them to be high also for the ARF homodimer. We expect that the separation of scales between these dimers reflects the role of ARF-IAA in the sequestration while the function of ARF 2 is associated with ARF's downstream targets (see last sub-sections).     The consensus sequence used is TGTCTC. We defined as clusters the groups of motifs whose distance is compatible with one turn of the DNA helix.  [2,13] and the insights obtained from our dynamical analyses in this paper. For the toggle switch, we used both results from FLAIR and the analysis reported in ref. [21].

Toy model without mass-ac�on law
Here we define the toy model without mass-action law: we define the variables' list, their derivatives and the steady-state variables. We fix the parameters and define the ODEs for the regulated and unregulated case.

Linear Response Analysis
We hereby implement the Linear Response analysis. We first define the jacobian matrices for the unregulated and regulated case respectively and we then plot the components of the Jacobian matrix's exponential for ARF, IAA and auxin describing their relaxation with respect to auxin perturbations.

Toy Model with Mass-Ac�on Law
Here we define the toy model with mass-action law: we define the variables' list, their derivatives and the steady-state variables. We fix the parameters and define the ODEs for the regulated and unregulated case.

Linear Response Analysis
We hereby implement the Linear Response analysis. We first define the jacobian matrices for the unregulated and regulated case respectively and we then plot the components of the Jacobian matrix's exponential for ARF, IAA and auxin describing their relaxation with respect to auxin perturbations.    Definition of the unregulated case: we take the rate of transcription at the basal level of auxin

Plots of time courses with regulation
We use the steady states as initial conditions, with an additive perturbation in auxin as previously Linear and nonlinear part of the dynamics and inclusion of parameters

Jacobian of the dynamics
We hereby linearize the dynamics around the steady state and define the Jacobian matrices for the negative feedback case, the unregulated case and the model without IAA-homodimers. We define the range of variation of auxin influx We initialize the steady states by solving the self-consistent equation for IAAss Hereby we define the lists of the steady states for IAA, ARF and auxin for the negative feedback and the unregulated case, respectively We start varying the auxin influx and find the solution to the self-consistency equation. We use this result to recompute all the steady states at each step

Linear Response analysis
We hereby consider S auxin = 0.02 and compute the corresponding steady states We hereby verify that the eigenvalues of the Jacobian are mostly real and they are all negative (this applies to the real part when complex), i.e. the time evolution is given by a decay

Linear Response Analysis
Hereby we define the value of the amplitude of auxin influx perturbation we are going to apply and compute the corresponding steady-state value of IAA.
We compute the Jacobian matrices with this steady states for the negative feedback and the unregulated case.

PCA on the Jacobian
We implement the Principal Component Analysis described into the Supplementary Material to understand whether some species may be clustered according to their relaxation time after auxin perturbation. This will consist in computing the eigenvectors of the two Jacobian matrices, i.e., negative feedback and unregulated case, and plot the components corresponding to the modes associated to the first two highest eigenvalues.

Plots of the temporal responses
We compute first the exponential of the jacobian matrices.
We plot the components corresponding to the response of IAA, ARF and auxin with respect to variations in auxin.

Toggle-switch implementa�on
Hereby we define the step-size perturbation for the auxin influx and we implement the equations for the auxin network with this perturbation. We define hereby the dynamical equations for the toggle switch and set the steady-state values as initial conditions for the integration with a step-size perturbation. We prepare the final big system of ODEs by merging together the equations of the auxin network and the toggle switch along with all initial conditions. eqssysneg = Joineqspertneg, toggleswitch, icsneg, icstoggleneg; eqssysconst = Joineqspertconst, toggleswitch, icsconst, icstoggleconst; We integrate the equations for the negative feedback and the case without regulation.

Self-consistency equa�ons and plot
Here we define the steady-state concentrations for the species we need to plot the self-consistency equation. We then define the two sides of the equations and finally plot them for different values of auxin influx.