The core mechanism: sequestration vs. active repression.

To begin with, we shall consider the two main mechanisms by which auxin can lead to transcription of its target genes. Indeed, the overall modus operandi of the system consists in releasing the ARF from IAA to induce transcription. This de-repression acts at two levels: ARF are released from ARF:IAA dimers, corresponding to our sequestration module, and the decrease of the IAA population also reduces transcriptional repression, corresponding to our active repression module.

As there is no clear evidence in the literature allowing to exclude one of these paths, we will start our investigation by exploring their relative influences on the input/output behaviour of the whole system. To do so, one introduces a new parameter λ ∈ [0, 1] such that λ = 0 and λ = 1 correspond to a situation without sequestration and a situation without active repression, respectively. Intermediate values of λ correspond to both mechanisms being operative. Specifically, one replaces in Equation (5):

• all terms α_AI IA−θ_AI D_AI by λ(α_AI IA−θ_AI D_AI),
• all terms α_{GA I} IG_A−θ_{GA I} G_AI by (1−λ)(α_{GA I} IG_A−θ_{GA I} G_AI),
• all terms α_{DAI G} GD_AI−θ_{DAI G} G_AI by (1−λ)(α_{DAI G} GD_AI−θ_{DAI G} G_AI).

Then, one wishes to consider the effect of varying λ on the different input/output functions (1)–(4). How these functions depend on both x and λ can be interpreted geometrically as a two dimensional landscape. Since other parameters of the system may in general have a strong effect on its input/output behaviour, we have in fact considered such landscapes for different values of the parameters. Because (5) comprises over 20 parameters, a complete exploration is excluded, even using advanced computational techniques and hardware. Since decay rates are the most well known parameters of the system, they were fixed close to experimentally observed values, as reported in Table 1 and Section “Models”. On the other hand, binding and unbinding rates are not known and could in principle vary by several orders of magnitude. We thus chose to vary these parameters, within plausible ranges based on the literature [32].

Even with this restriction the number of parameters to vary was still computationally daunting, leading us to additional simplifying assumptions: we distinguished between protein-protein and protein-promoter binding reactions, but within each two types we assumed that all association (resp. dissociation) rates of the form α_XY (resp. θ_XY) were equal.

The parameters α_XY were varied over 5 orders of magnitude. By default, from comparable cases found in the literature they took values (in nM⁻¹.min⁻¹ for protein-protein association and min⁻¹ for protein-promoter association) in the sample set whereas the parameters θ_XY were varied over 6 orders of magnitude, taking default values (in min⁻¹ for protein-protein association and nM.min⁻¹ for protein-promoter association) in

The value 0 was not included for α_XY because it corresponds to an extreme situation where no dimers can form, which leads to an absence of response to any auxin input. Even with the assumptions made at this point, exploring all combinations involves 900 cases, which would have led to an unreasonably long computational time. We thus restricted the exploration by imposing that protein-protein and protein-promoter constants vary by at most one order of magnitude.

We use a notational convention to describe all the considered cases. Denote α_P < α_AG a situation where the association constants for protein-promoter vary in S_α and association constants are one order of magnitude smaller, hence varying in S_α/10 = {0.0001, 0.001, 0.01, 0.1, 1.0}. Similarly, we denote θ_P < θ_AG when the dissociation constants for protein-protein binding reaction are one order of magnitude smaller than constants for protein-promoter dissociation (and thus vary in S_θ/10 instead of S_θ). Table 2 summarizes all the cases we have explored, with the natural generalization of the notation just described:

Since for each case the ratio between protein-protein and protein-promoter constants was fixed, this left us with two global parameters α and θ that were varied in S_α and S_θ respectively, used as ranges for the higher constants in cases where they differed between protein-protein and protein-promoter.

To chose a relevant range of values for x, we performed a preliminary exploration to ensure that a non-negligible response was obtained in all or most explored parameter regimes. This led us to empirically define the range x ∈ [0, 2000].

In summary, for each of the 9 cases above we computed numerically the values of the 4 functions (1)–(4) for a 200×200 regular grid on the domain (x,λ) ∈ [0, 2000]×[0, 1]. This operation was repeated for a coarse grid of values taken by the two global parameters α and θ. The complete computation was distributed on 8 cores (2 cores per machine on 4 recent stations having between 6Gb and 24Gb memory per core) and took over 597h ≈ 24.88 days to be complete. The results for the default situation, “α_P = α_AG” & “θ_P = θ_AG” are shown in Figs. 3. Remarkably, the corresponding figures of the 8 other cases A-H (S1 Additional Information), did not differ significantly from the default parameter choice, whose main features are discussed below.

The landscapes shown on Fig. 3 present several features. First at the coarse grid level it appears that except for the response time, the output functions are far from uniform on the coarse grid: some regions in the (α,θ) space are more favourable than others. Both the absolute and relative responses reach higher values for $\frac{\theta }{100}\le \alpha \le \theta$, corresponding to a dissociation constant $K_d\doteq \frac{\theta }{\alpha }$ in the range [1nM, 100nM], i.e. strong to moderately strong dimers [32]. Within this range, the landscapes corresponding to K_d = 100 in fact have a globally high absolute response but a lower relative response than for more stable dimers. Along the x axis, higher values of the output functions are increasing, reaching a maximal value at a critical auxin level after which they tend to plateau.

In the range K_d ∈ [1nM, 10nM] of stronger dimers notably, there seems to be a minimal λ value below which ρ_rel landscapes stay low and never reach their plateau. For θ ≤ 0.01 for instance, i.e. the 3 lower rows in Fig. 3(a), there clearly seems to be a restricted range of λ values close to λ = 1 leading to higher responses. Though less striking, this phenomenon also appears on ρ_abs landscapes, where one also notices (e.g. on the row θ = 0.001) a value of λ below which no plateau is reached.

Although the ranges of K_d values changes from the above, the qualitative features of the default parameter choice were found to be conserved for the cases A-H. This indicates that variations within the association and dissociation rates may have a tuning effect on the values at which the different outputs are maximal, but do not affect the following general feature: the relative and absolute response landscapes present an optimum plateau, which lies above a minimal value of λ > 0 and does not include λ = 1.

In the case of the sensitivity σ, there appears to be a specific and narrow range of λ and x values at which sensitivity is maximal, Fig. 3(c). This specific range is dependent on the (α,θ) values on the coarse grid, a fact reminiscent of the role of a combinatorics of dimers (and thus dimerization kinetic constants) as a means to modulate auxin response [7]. This fact will be discussed further in later sections. In the present case, it is also important to notice that the maximal values of σ reached on all grids do not exceed 1. In other words, for none of the parameter values considered in Fig. 3 is the system ultrasensitive.

The previous observation suggests a closer inspection near the limit values λ ∈ {0,1}. It appears that these limits systematically correspond to a drastic decrease in the output functions, as shown in Figs. 4 and 5. On these figures on can see sharp transitions to lower values for the three functions ρ_abs, ρ_rel and σ near λ = 0 and λ = 1. If the sharpness might look like the result of numerical inaccuracies, further increases of the resolution can show that they in fact correspond to a smooth, yet abrupt, decrease (data not shown). Here again, these features were found also in the cases A-H shown in S1 Additional Information.

Consequently, one global conclusion is that λ does not seem to be a crucial parameter for the system, insofar as it does not take one of the extreme values 0 or 1 which both strongly reduce all the considered output functions. In other words, the core mechanism underlying auxin signalling seems to be a combination were both sequestration and active repression are required for the response to be high and robust.

This led us to remove the λ parameter from later investigations, where we explore more finely the influence of other parameters on the different output functions, as detailed in the following. The results above, will be used to select values of α and θ rates in the following, with the aim to guarantee a fairly high level for all the output functions (provided x is not too low).

Because ARF:IAA dimers have already been considered, as they constitute the core of what we have termed sequestration, we now turn specifically to the influence of homodimers (in the sense that they are composed of proteins of the same family, IAA or ARF).

Role of ARF:ARF dimers.

As ARF:ARF dimers have been found on promoters [23], a potential role of these dimers one may expect on intuitive grounds is to increase the steepness, or in fact the sensitivity of the transcriptional response.

In the response landscapes shown in Fig. 3 the sensitivity output is never greater than 1, i.e. the system is not ultrasensitive for the parameters explored in these landscapes in apparent contradiction with results suggesting that systems where a molecule (here IAA) sequestrates an other (here ARF) are known to present ultrasensitivity [32, 33].

Probably the strongest restriction imposed on the parameters used in Fig. 3 is the fact that all association (resp. dissociation) reactions have the same kinetic constant α (resp. θ). Although not a requirement, it has been observed in [32] that a disparity in decay rates can strengthen ultrasensitivity. This led us to consider alternative parameter values, where different binding reactions have distinct rates. None of our attempts showed ultrasensitive behaviour. Although the immense range of possible choices prohibited an exhaustive exploration, this indicates at least that ultrasensitivity may not be a dominant behaviour.

Finally, we investigated an intriguing possibility, suggested by recent experimental results on the structure of the domain III-IV of ARF and IAA that mediate protein-protein interactions. Several recent papers show indeed that ARF protein have the ability to regulate transcription with cooperativity, and can agglomerate as multimers of higher order than 2 proteins [4, 26, 27]. Even though the function of these structural properties would be a modelling research topic in itself, we considered a very simplified representation of these properties in our model. Namely, we introduced an additional parameter n_A representing the number of ARF proteins agglomerated in a multimer. This parameter can also be interpreted as a measure of cooperativity. The following alterations of Equation (5) were performed:

• terms α_AA A²−θ_AA D_AA were replaced by ${\alpha }_{AA}{A}^{n_A}-{\theta }_{AA}{D}_{AA}$ in the equations $\dot{A}$ and ${\dot{D}}_{AA}$.
• terms α_AGA AG_A−θ_AGA G_AA were replaced by α_AGA A^n_A−1 G_A−θ_AGA G_AA in the equations ${\dot{G}}_A$ and ${\dot{G}}_{AA}$.

These alterations are clearly phenomenological and a more realistic account would require the introduction of new variables for each multimer involving p ≤ n_A ARF proteins, see [34] for a related discussion. In this new formulation, D_AA does not represent ARF:ARF dimers, but multimers composed of n_A ARF proteins.

As parameters like n_A are known to increase the steepness of the transcriptional response as a function of ARF, one may expect ultrasensitivity to occur for the full pathway for high values of n_A. To assess the validity of this claim, we calculated the sensitivity σ(x) for a sample of values x ∈ [0, 300], for different choices of α and θ (common for all binding reactions). This led us to empirically select the values (α,θ) = (10,0.001), for which the sensitivity was higher. For this choice (corresponding to very strong dimer associations, i.e. a very low K_d = 10⁻⁴), we computed the curve (x,σ(x)) for increasing values of n_A as shown in Fig. 6. The sensitivity did raise above 1 for higher values of n_A and at lower values of x, below 10.

Role of IAA:IAA dimers.

The role of IAA:IAA dimers is not intuitively obvious. This led us to consider a highly simplified model in which a single population of proteins is produced, degraded and can form homodimers, in order to derive first predictions which we then test in our complete model of the auxin pathway. Note also that the existence of higher orders IAA complexes [26, 27] has been suggested but has not been considered here for the sake of simplicity (but see discussion).

As is very often considered, we assume that dimers reach equilibrium at a faster rate than the production (transcription + translation) or degradation of proteins. This leads to the following equation (cf. Suppl. Info. for details): where the parameter γ is proportional to the dimer association rate and π_I is the steady state value in absence of homodimerization. A first calculation shows that the unique steady state I* of this equation is a decreasing function of γ, and in particular it is always lower than π_I for γ > 0.

Since increasing γ results in a lower steady state value, a given percentage of this steady state might be attained quicker even though the system evolves at the same speed. In order to assess the response time we rescaled the system to ensure the steady state is 1 regardless of γ:

Due to its simplicity, a closed form solution can be found for this equation, allowing to calculate a characteristic response time. From these calculations (cf. Suppl. Info.), one deduces that the characteristic time is a decreasing function of γ, which tends to 0 as γ → ∞. In other words, the stronger the homodimer formation rate, the faster the system reaches equilibrium (for instance after a perturbation), in principle arbitrarily fast if the dimers are strong enough.

The discussion above concerns a highly simplified description of the whole system, and is discussed only for it suggests a potential role of the IAA:IAA dimers improving the response time of the system. To confirm this suggestion, we return to the full model (5) and fixing all other parameters, assess the influence of α_II on the system’s response time.

More precisely, we varied all association and dissociation constants in the same coarse grid as in Fig. 3. For each of these values, we ran simulations of the model (5) with a constant input of auxin $\bar{x}=3000$ for a sample of values of the α_II association constant in the range [0, 100], recording the response times τ(x).

Since the analysis of the simplified model above considered a response time with a normalized steady state, we also calculated a response speed, defined as the ratio of the absolute response ρ_abs over the response time. When expressed as a ratio, this speed would for instance be reduced if the steady state is approached within the same time span, but has a lower value.

As anticipated, the response time did actually often increase with α_II, see Fig. 7(a), but this was generally correlated with an increase of the steady state value (not shown). The response speed, on the other hand, displayed a clearer pattern shown in Fig. 7(b): for most (α,θ) in the region θ > 0 and α > 0.1 the system responded in shorter time (at higher speed) as α_II was increased and overall it appeared that any response with α_II > 0 was faster than with α_II = 0. This confirmed that the observation made on the simplified model above is still valid for wide ranges of parameters in our complete model.

This section suggests that the overall auxin signalling dynamics is faster when IAAs are able to form dimers, with the caveat that this requires other parameters to belong to specific (wide) ranges. This type of speed up induced by a nonlinear feedback term has been previously observed [35], although in that case it was self-activation, whereas we deal here with a process whose effect is ultimately self-inhibitory.

Forms of auxin induced degradation.

As discussed in the first section it is assumed in our main model that dimers comprising an IAA unit can lose this unit upon auxin binding, releasing the other protein forming the dimer. However, other scenarii are also plausible. In particular, we have investigated the following alternative mechanisms:

• (m1) Dimers are not affected by auxin.
• (m2) Dimers bearing an IAA are entirely degraded by the proteasome as a consequence of the IAA binding TIR1 through auxin.

The mechanism (m1) may correspond to a situation where the structural conformation of IAA:X dimers prevents binding of auxin to the domain II of the IAA, whereas (m2) would account for ARF:IAA and IAA:IAA being obligate dimers.

In terms of the ODE model (5), (m1) can be modelled by suppressing all the terms δ_II xD_II and δ_AI xD_AI from the system. Note that in an extreme variant of this mechanism, dimers are not degraded at all (or at rate negligible compared to dissociation), in which case δ_II and δ_AI are simply set to zero.

For the mechanism (m2) on the other hand, the terms δ_II xD_II and δ_AI xD_AI are suppressed from the rates of I and A, but remain as decay terms in the rates of D_II and D_AI.

To begin exploring the properties of these two mechanisms, we computed the same landscapes as in Fig. 4 (see S2 Additional Information) for the two alternative models (m1) (for its extreme variant) and (m2). The overall effect for (m1) was not obvious, all the main characteristics discussed in the previous section remaining unchanged.

On the other hand, in the case (m2) a noticeable effect could be observed: in some parameter ranges, auxin could act as repressor, despite our default situation involving an activator ARF. These exploratory simulations seemed to indicate that the phenomenon was mostly seen near λ = 1, i.e. in a regime where the effect of auxin is mostly due to ARF being sequestered in ARF:IAA dimers. In agreement with the previous section, we sought to reproduce this phenomenon in a model without λ (i.e. with both core mechanisms).

A regime that would be comparable to high λ values is one in which the repressive action of ARF:IAA dimer is weak, which corresponds to h being close to or maybe higher than h_AI. We thus considered a variant of the default parameters shown in Table 1, fixing h = 0 and h_AI = 1 (reversing the default values) and leaving all other parameters unchanged. For this choice we found, as in the landscapes discussed above, that auxin could act as a repressor in the case of model (m2), but not in model (m1). Some typical curves are shown in Fig. 8. Similar curves were found for pairs (α,θ) such that the dissociation constant $K_d=\frac{\theta }{\alpha }\ge 1000$. For other values of the dissociation constants, the absolute response was systematically higher with mechanism (m1) than with (m2), suggesting that even in a regime where auxin activates the transcription of its targets, its efficacy is dampened in the case of (m2). Remarkably, the repressor effect was only observed for a limited range of auxin levels and in every model transcription is increased upon high inputs of auxin.

An intuitive explanation of the repressive action of auxin with degradation mechanism (m2) is that ARF:IAA degradation (induced by auxin) may result in an overall decrease of the ARF population. This can result in a reduction of the transcriptional activation by ARFs and their dimers. If the basal level of transcription is itself low, the end result of an input of auxin might be a repression of the target gene, even if the main population of ARF is composed of activators. Note that this phenomenon does not appear with other forms of dimer degradation and is furthermore limited to a peculiar set of parameter values (high K_d, i.e. weak dimers).

Role of feedback.

Our primary definition of feedback is due to IAA proteins having genes including ARF binding sites. Beyond this primary loop, there exists experimental evidence that some ARF proteins are also regulated by auxin [17]. In the case of activator ARFs, this results in a positive feedback loop, whereby an increase of ARF level ultimately increases further the production of ARF proteins (by inducing an increase on the transcriptional activation of auxin responsive genes, which include ARF coding genes).

Feedback loops have a crucial influence on the qualitative dynamics of a system, and it is established that the presence of a negative feedback loop (resp. positive feedback) is necessary for the occurrence of periodic solutions, i.e. oscillations (multiple steady states, i.e. multistability) [36–38]. These results are necessary conditions only, and the presence of bistability or oscillations, though it can be excluded in absence of any feedback loop, can only be confirmed by an analysis which is often non trivial.

In fact, the so-called mixed feedback loop is essentially a sequestration system with an additional feedback loop, and has been shown in [39] to have the potential for bistability when the loop is positive and oscillatory behaviour when the loop is negative. Although suggestive, this result can only be transposed to our situation with great care, as the auxin pathway presents several differences with the mixed feedback loop. Firstly, in the latter it is assumed that only the (sequestered) monomer is a regulator, whereas both ARF and ARF:IAA regulate transcription. Secondly, in the auxin pathway dimers are formed not only between ARF and IAA but also amongst each of these two families.

Concerning the negative feedback, a previous study has shown that the auxin pathway could in principle present stable oscillations [40]. Since this study involved a model that had a similar level of detail as our present system, we did not attempt a thorough exploration of oscillations, as it would have been redundant with [40].

We considered the role of a positive feedback on ARF+, which we implemented as follows in Equations (5):

• π_A is replaced by π_A R, becoming dependent on auxin regulated transcription,
• the production rate of I is a constant π_I.

The last assumption was not strictly speaking required to investigate positive feedback, but was used as a simplifying assumption which allowed to perform analytical calculations. Namely, in a specific parameter regime where all dimers are supposed to occur more rapidly than protein production and degradation (using a quasi-steady state assumption), we were actually able to characterize analytically some parameter values for which the system could admit multiple equilibria, see S1 Text and Fig. 9. These results were obtained under a quasi-steady state assumption, only valid in fairly specific regions of the parameter space of (5).