Abstract
Negative feedback control is a ubiquitous feature of biochemical systems, as is time delay between a signal and its response. Negative feedback in conjunction with time delay can lead to oscillations. In a cellular context, it might be beneficial to mitigate oscillatory behaviour to avoid recurring stress situations. This can be achieved by increasing the distance between the parameters of the system and certain thresholds, beyond which oscillations occur. This distance has been termed resistance. Here, we prove that in a generic threedimensional negative feedback system the resistance of the system is modified by nested autoinhibitory feedbacks. Our system features negative feedbacks through both inputinhibition as well as outputactivation, a signalling component with mass conservation and perfect adaptation. We show that these features render the system applicable to biological data, exemplified by the high osmolarity glycerol system in yeast and the mammalian p53 system. Outputactivation is better supported by data than inputinhibition and also shows distinguished properties with respect to the system's stimulus. Our general approach might be useful in designing synthetic systems in which oscillations can be tuned by synthetic autoinhibitory feedbacks.
1. Introduction
Negative feedback control is a fundamental and a ubiquitous feature of biochemical systems [1–6] and can mediate adaptation [7–10], stabilize the abundance of biochemical species [4,11,12], induce oscillations [3,5,13–16] and accelerate response times [11,17]. In fact, negative feedbacks have been observed in a wealth of biological systems ranging from mammalian cell cycle [13,18] to bacterial adaptation [8,19] and stress response in mammals [20] and yeast [21,22]. Another ubiquitous principle of biochemical systems is time delay between a signal and its response, which can, for example, be caused by the time needed to transcribe and translate biochemical information into cellular compounds. It is a longstanding theoretical result that negative feedbacks in conjunction with time delay can lead to stable oscillations [15,16,23]. Oscillatory behaviour brought about by delayed negative feedbacks has been observed and analysed in a range of biological systems, for example, the mammalian p53 system [24–26], and the NFκB system [27–29]. Both the p53 as well as the NFκB system mediate adaptation to external stimuli and stress, such as, for example, DNA damage. It is comprehensible that it might be beneficial to mitigate oscillatory behaviour during adaptation in order to avoid recurring stress situations. This can be achieved by moving the steady state of the system far away from certain thresholds, beyond which oscillations occur. The distance between such thresholds and the parameters of the system has been termed resilience and/or resistance. The larger the resistance of a system, the better perturbations in external or internal conditions, i.e. the systems parameters, can be absorbed, which otherwise would trigger a change in stability properties of the system [30,31]. There exist several definitions of resistance in the literature [31]. Here, based on [31], we will refer to ‘resistance’ as the system's response to perturbations of parameter values. As a quantitative measure of resistance, we use the Euclidian distance of the parameter vector to a critical threshold, beyond which stability properties of the system change. The larger this distance, the more resistant is the system.
In a recent study, using a parametrized mathematical model, evidence was presented that during osmoadaptation in yeast, which is largely mediated by a delayed negative integral feedback, the potential of oscillatory behaviour is reduced by introducing nested direct negative feedbacks [21]. Thus, there is evidence that in a concrete biological system nested negative feedbacks can increase the resistance of a biochemical network.
Here, using a generic threedimensional model for integral negative feedback control of a biochemical network, we explore whether coupling autoinhibitory and delayed negative feedbacks might be a general cellular mechanism to increase resistance of a system. Our system has several distinguished features which generalize and extend former studies [15,16,32]:

— our models include components, which resemble posttranslational modification of proteins conserving the total protein abundance (mass conservation), rendering the system more general and, at the same time, more realistic, especially with respect to signalling cascades;

— our models use integral feedback properties, i.e. some state variables of the system robustly track their desired values independent from the input [8], which also renders the systems more applicable to realistic situations, as shown below; and

— in our models, the delayed negative feedback may operate through both by inhibiting sensor inflow, like in [15], and by activating sensor outflow, like in [9]. This has been termed input and output control [10]. Here, we refer to these two control types as inputinhibition and outputactivation, respectively. The type of used delayed negative feedback has important implications with respect to what actually stimulates the system. We show that an outputactivation feedback mechanism is better supported by a range of data.
We prove that, in our systems, stable limit cycle oscillations can occur owing to a Hopf bifurcation. Further, we prove that the parameter region, where oscillations occur, can be reduced by introducing autoinhibitory feedbacks. Thus, by nesting autoinhibitory negative feedbacks into delayed negative feedbacks the structural stability of the system can be altered. This is true, in general, for inputinhibition systems. However, there exist certain limitations for this effect in outputactivation systems. We also provide computational evidence that the sensitivity of the steady state with respect to parameter perturbations is decreased in system with nested autoinhibitory feedback. We apply our generic model to the high osmolarity glycerol (HOG) system in yeast, and the mammalian p53 system demonstrating the applicability of our general framework to concrete situations. We propose that this simple framework can be used to design synthetic systems in which oscillatory behaviour can be tuned by nesting direct, possibly autoinhibitory, and delayed negative feedbacks.
2. Results
2.1. The model
In the following, we consider the threedimensional system 2.1with positive parameters k_{0}, k_{1}, k_{2}, k_{3}, k_{4}, k_{5}, k_{6}, y_{T}, κ_{y}, κ_{z} and nonnegative initial conditions x_{0}, y_{0} and z_{0}. External perturbations are simulated by modifying the value of k_{0} which otherwise mimics a basal stimulation of the system. This system can be even more generalized by replacing all linear functions x and y by smooth strictly increasing functions (see the electronic supplementary material). Figure 1 displays a wiring diagram of the system (2.1) also indicating different alternative model formulations that were tested in the application described below.
In general, x(t) can be considered as a cellular sensor, which reacts to an external stimulus k_{0}. The component y(t) mimics a signal transduction, which relays the input signal coming from the sensor x(t) to the response component z(t), which in turn negatively feeds back into the sensor x(t) via f_{3}(z) or via q(z). We assume f_{1},f_{3}, and h to be strictly increasing functions on vanishing at 0. The function q(z) is a smooth, positive and strictly decreasing function with q(0) = 1. The functions f_{3}(z) and q(z) represent the overall feedback in the system, respectively. In the examples below, we consider them to be mutually exclusive and refer to a feedback through q(z) as inputinhibition, and a feedback through f_{3}(z) as outputactivation. The functions g(y, κ_{y}) and H(z, κ_{z}) are smooth, positive and decreasing in both arguments with . The functions g(y, κ_{y}) and H(z, κ_{z}) may represent autoinhibitory feedbacks, where we later consider the additional parameters κ_{y} and κ_{z}, which, among other parameters, shape the form of these functions. We call them autoinhibitory, because the inhibition of a component depends only on itself and not on other system variables.
The term (y_{T} − y) leads to a model with an a priori bound y for the second component in case its initial value y(0) is in [0, y_{T}]. This term is obtained by reducing a fourdimensional system with mimicking a reversible posttranslational modification, such as, for example, phosphorylation, of a protein y that does not affect the total protein abundance y_{T} (mass conservation).
Note that the nonnegative orthant is positive invariant for all models and, therefore, all models are biologically sound in the sense that no negative values for the components can occur. For further details, please refer to the electronic supplementary material.
Taken together, we analyse a generic model that comprises a range of special cases that have been addressed in the literature, for example, the Goodwintype models [15,16,32], but also addresses models that have not been thoroughly analysed yet, especially the outputactivation models, which will be shown to be especially relevant in concrete situations.
2.2. Integral feedback property
Integral feedback control is an engineering strategy that is supposed to ensure that the output of a system always adapts to its desired value independent of noise and of perturbations of the system parameters [8]. For twodimensional systems, it has been reported that the kinetic nature of h(z) is important in this respect; for example, mass action kinetics for h(z) is not sufficient to obtain perfect adaptation for y(t) [9,10]. For convenience of the mathematical analysis, we approximate a zeroorder h(z) by a smooth function with h(0) = 0, and h(z) = 1 for z ≥ a > 0 and that is strictly increasing on (0, a) and require that an equilibrium exists with . Thus, and, for , the equilibrium of the second component is constant and independent of the input signal in this limiting case.
Taken together, the system (2.1) approximates a perfect adaptor for zeroorder h(z) with respect to y(t). Note that the above approximation is only a theoretical one ensuring that the solution stays in the positive orthant , unlike the approximation K_{z} = 0 for Michaelis–Menten type h(z) = z/(K_{z} + z) for which negative solutions can occur. For concrete situations, it suffices to assume a sufficiently small K_{z} for the Michaelis–Menten type h(z) implying .
2.3. The outputactivation system is stimulated by the difference between the internal and external state
Zeroorder kinetics for f_{1}(x) has important implications with respect to what is actually sensed and integrated as error function by the system. Let us assume , a linear feedback function f_{3}(z) = z, and that the sensor x(t) is in quasiequilibrium with respect to the response variable z. Further, we approximate a zeroorder f_{1}(x), for example, , by a smooth function with f_{1}(0) = 0, and f_{1}(x) = 1 for and that is strictly increasing on (0, a) and require that at the equilibrium . Thus, and 2.2is a linear function of z. Thus, for zeroorder and linear outputactivation, the system (2.1) is stimulated by the positive difference between the external stimulus k_{0} and the scaled response variable z, i.e. the internal state. As above, this approximation is introduced for convenience of the theoretical analysis only. In real situations, it suffices to assume a sufficiently small K_{x} as in figure 2.
In the case of mass action kinetics for f_{1}(x), 2.3
is related to the ratio between the external stimulus k_{0} and the response variable z. A similar form for as in (2.3) is obtained in a model with inputinhibition, i.e. k_{2} = 0 and which was assumed in the classical model of Goodwin [15,16,32]. Here, we obtain 2.4In figure 2, we display quasiequilibrium sensor values as a function of the response variable z for the system (2.1) with either a Hilltype inputinhibition term as in [15] and k_{2} = 0, or a zeroorder outputactivation term with linear feedback in z, i.e. with , respectively.
Note that in a system where the feedback is mediated through outputactivation by a linear feedback with zeroorder degradation, i.e. and , the response threshold k_{0}/k_{2} is well defined for small K_{x} and can directly be tuned through k_{2} (figure 2). Similarly, for systems with Hilltype inputinhibition, i.e. k_{2} = 0 and , the response threshold is well defined for large Hill factors h (grey dashed line in figure 2), but can only indirectly be tuned through the halfsaturation constant K. In addition, the better the threshold is defined by large Hill factors h, the sooner the sensor activation saturates for decreasing z, whereas in a system with outputactivation, sensor activation is linear below the activation threshold k_{0}/k_{2}. Conversely, when in a system with Hilltype inputinhibition, the sensor approximates a linear function in z (grey solid line in figure 2), the response threshold is poorly defined.
Taken together, the type (inputinhibition or outputactivation) and kinetic nature (mass action or saturating) of the overall negative feedback determines the signal, which stimulates the system. In the case of zeroorder outputactivation, the system is linearly stimulated by the difference between internal and external conditions; otherwise, the stimulus is nonlinear and may even be steplike. In the examples below, we show that the data clearly support outputactivation models rather than support nonlinear Goodwintype inputinhibition models.
2.4. The Hopf bifurcation
In the system (2.1), stable limit cycles can occur due to a Hopf bifurcation. Shortly, any steady state of the system (2.1) is given by 2.5
We define G(y, κ_{y}) = (y_{t} − y)g(y, κ_{y}) and the parameter value 2.6
The Jacobian at the equilibrium is of the form
and has the characteristic Hurwitzpolynomial , with positive
So, any real eigenvalue of J is negative. The necessary and sufficient condition for a single pair of pure imaginary eigenvalues is , i.e. 2.7evaluated at . With
we consider (2.7) as an equation of the parameters k_{1} and that is to be solved in the form . The curve , given by the unique positive solution of the quadratic equation, derived from (2.7), 2.8evaluated at , indicates possible Hopfbifurcation points in the plane. The necessary and sufficient condition of the existence for a positive solution of (2.8) is , i.e. 2.9
Taken together, having chosen the equilibrium component we define and according to (2.5), solve (2.8) for , provided (2.9), and set according to (2.6). Then, the possible Hopfbifurcation point is given by for the critical parameters k_{1} and .
In this case, the transversality condition for a Hopf bifurcation can be generically fulfilled (see the electronic supplementary material, where also further details of the proof are supplied). Thus, the system (2.1) can show stable oscillations owing to a Hopf bifurcation.
A Hopf bifurcation can also occur for a more general class of systems, where all linear functions x and y in the system (2.1) are replaced by smooth strictly increasing functions and with or without the term (y_{T} − y) (see the electronic supplementary material).
2.5. Autoinhibition decreases the oscillatory
In the electronic supplementary material, we prove that the plane permissive for oscillations decreases with increasing autoinhibition either through g(y, κ_{y}) or H(z, κ_{z}), i.e. for the curve that divides the plane into regions with and without stable oscillations it holds 2.10
Again, there is a notable difference between models with inputinhibition, i.e. k_{2} = 0, and outputactivation, i.e. .
For models with inputinhibition relation, (2.10) is always true. This has been shown before for the classical Goodwintype models with Hilltype q, i.e. [32]. However, opposed to these classical models, where a high Hill coefficient (cooperativity) of h ≥ 8 was necessary to obtain oscillations, this is not necessary in our framework (2.1) (see the electronic supplementary material, figures S4 and S7).
For models with outputactivation condition (2.10) only applies, if , i.e. if at the equilibrium f_{1} is of zeroorder such that . This situation can be approximated, e.g. by low K_{x} values for . Note that this was also a prerequisite for the quasisteady state to be a linear function of z as in (2.2). However, to reduce the parameter region for oscillations for the outputactivation system, is sufficient, but not necessary (for examples, refer to the electronic supplementary material).
Taken together, we provide formal proof that the structural stability of the system (2.1) can be altered by introducing autoinhibitory feedbacks. Moreover, the distance between the bifurcation threshold and a given pair can be modulated by introducing autoinhibitory negative feedbacks. Thus, nested autoinhibitory feedbacks can modulate the resistance of the system (2.1). Figure 3 illustrates these theoretical results for an outputactivation system. For other examples, please refer to the electronic supplementary material, figures S2–S7.
In figure 3a, we show bifurcation curves in the plane, for different values of κ_{y}, i.e. with (κ_{y} = 40) or without (κ_{y} = 0) autoinhibitory feedback for a concrete system. The area below the curve, where oscillations occur, is reduced with κ_{y} increasing. The larger κ_{y}, i.e. the stronger the autoinhibitory feedback, the smaller the area below the curve. The dot in figure 3a indicates a concrete pair of , to which the computed eigenvalues in figure 3b and dynamics in figure 2c,d correspond. The distance between a point in the and the bifurcation curve can be interpreted as a measure for resistance. Note that a change in can also be interpreted as a change in parameter , because there is a 1 : 1 relationship between and (see (2.6)). Without autoinhibitory feedback (κ_{y} = 0), the bifurcation parameters are below the bifurcation curve and, consequently, we have a single pair of complex eigenvalues with positive real parts (black dots in figure 3b) corresponding to stable oscillations (figure 3c). With autoinhibitory feedback (κ_{y} = 40), the bifurcation parameters are above the bifurcation curve and, consequently, all real parts of the eigenvalues are negative (circles in figure 4b), and the system tends to a stable equilibrium (figure 3d).
For convenience, we conducted the theoretical analysis by parametrizing the system with respect to the steady state in and considered k_{1} as the bifurcation parameter. However, as illustrated by the computational analysis of the HOG system and the p53 system below, all tested parameters in the system may be taken as bifurcation parameters (figures 4 and 7; electronic supplementary material, S9 and S10). Accordingly, the introduction of autoinhibitory feedbacks reduces the region for oscillations for those parameters as well.
2.6. Application to the high osmolarity glycerol system
The HOG system in yeast mediates adaptation to a hyperosmotic shock and is one of the beststudied eukaryotic signalling pathways [33]. Several mathematical models of different complexity have been developed for this system [21,22,34,35]. In short, the signal that triggers response, and adaptation is supposedly related to volume [36,37], which, in turn, is proportional to the difference between internal and external osmotic pressure [38]. The signal coming from the membrane is transduced via a stressactivated protein kinase (SAPK) cascade, which culminates in the activation the SAPK Hog1. Hog1 translocates to the nucleusactivating transcription factors that lead to the upregulation of glycerol production, which, in turn, increases the intracellular osmolarity and turgor, thereby mediating adaptation. There is evidence that in this system oscillatory behaviour might indeed be avoided by nested negative feedbacks [21]. Therefore, we tested whether our general framework is supported by available data, and, whether the data support autoinhibitory feedbacks. We fitted different candidate models representing different hypothesis about the underlying biochemical mechanisms and ranked them according to the Akaike information criterion (AIC). The candidate models were specified versions of our general framework (2.1) (figure 4a): 2.11where and indicate model alternatives. Here, x represents a putative sensor of volume change or the difference between internal and external water potentials, i.e. z and k_{0} + NaCl, respectively. The component y represents the adaptive phosphorylated Hog1 and z represents the integrator glycerol (figure 4a). Specifically, we tested the kind of delayed feedback, i.e. inputinhibition (δ_{1} = 1) or outputactivation (δ_{1} = 0), the existence of autoinhibition in the signalling component y (δ_{3}) and two different kinetics for f_{1}(x) and h(z), respectively, i.e. mass action (δ_{2}, δ_{4} = 0) or Michaelis–Menten kinetics (δ_{2}, δ_{4} = 1). The combination of all these model alternatives yielded 12 different models. Here, we assumed f_{3}(z) = z and H(z,κ_{z}) ≡ 1. For parameters of the bestranked model, please refer to the electronic supplementary material, table S2. For more details on the model and parameter estimation, please refer to the Methods section and the electronic supplementary material. A COPASI implementation of the bestranked model together with the data and an systems biology markup language (SBML) [39] version is also provided in the electronic supplementary material. The results of fitting and ranking are displayed in table 1.
The model with zeroorder kinetics in f_{1}(x), autoinhibition in y and mass action kinetics in h(z) is ranked best (model no. 6, δ = (0, 1, 1, 0)). This corresponds to a model, which senses the difference between external and internal osmolarity and shows no perfect adaptation. This model is closely followed by the same model, but with perfect adaptation (model no. 8, δ = (0, 1, 1, 1)). The best model without autoinhibition in y is on the third place (model no. 2, δ = (0, 1, 0, 0)). The best two models are able to recapitulate Hog1 phosphorylation and intracellular glycerol data for a range of different conditions (figure 5). Model no. 2 can also well recapitulate Hog1 osmotic stress and glycerol data, but cannot recapitulate well the Hog1 inhibition experiment (see electronic supplementary material, figure S8). The models with an inputinhibition did not give a good approximation to the data compared with the bestranking models (table 1).
Apparently, both the system with and without autoinhibitory feedback can show adaptive behaviour. Analysing the stability of the steady state as a function of the parameters, it becomes obvious that the parameter region, where oscillations occur, is much more distant from the actually fitted parameters for the model with (model no. 6) than for the model without autoinhibition (model no. 2). Thus, the resistance of the adaptation is increased in the system including the autoinhibitory feedback. Perturbing the initial steady state, which was also set in this case, had no influence on the stability, i.e. the system is resistant with respect to a change in initial steadystate concentrations (see the electronic supplementary material, figure S9). In figure 6, we plot the stability regions of the steady state in the twodimensional and plane. For other parameter combinations, please refer to the electronic supplementary material, figure S9. Notably, the parameters of both the system with autoinhibition (model no. 6) and the system without autoinhibition (model no. 2) are rather similar (black and grey dots in figure 6 and the electronic supplementary material, S9). This indicates that the stability of a system can be modified by changing the system's structure by autoinhibition without significantly affecting other system parameters and, therefore, its dynamics (figure 5 and the electronic supplementary material, S8).
It can be anticipated that in homeostatic adaptive systems the steady state should be robust against parameter perturbations. It has been shown that negative feedbacks can increase the robustness of the steady states with respect to input noise and parameter perturbations [4,11,12]. We hypothesized that nested autoinhibitory feedbacks can increase the robustness of the steady state. Therefore, we compared the steadystate variability with respect to parameter perturbation of the best model with feedback (model no. 6) and the best model without feedback (model no. 2) after an osmotic shock of 0.2 M NaCl in a Monte Carlo analysis. Specifically, we perturbed all free parameters of the system simultaneously by sampling 1000 times from a uniform distribution ranging from half to double of its original value. Subsequently, we calculated the distance between the original fitted steady state and the perturbed steady states (figure 7).
The variance of the distance between the original and the perturbed steady states of the sensor x is significantly smaller (p < 0.01) for the model with autoinhibitory feedback (model no. 6) compared with the model without autoinhibitory feedback (model no. 2; figure 7). In addition, the respective variance for components y is smaller (p < 0.01; the electronic supplementary material, table S1). For the response z, no significant difference was detected. Thus, in this concrete case, the autoinhibitory feedback increases robustness of the steady states after osmotic shock for the sensor and activated Hog1.
Taken together, our threedimensional framework is able to recapitulate well a range of data for different conditions for the HOG system. Model discrimination suggests that in the HOG system there are autoinhibitory feedbacks nested within the glycerolmediated feedback, and the system is stimulated by the difference between internal and external conditions. This is well supported by other studies [17,21]. Moreover, the autoinhibitory feedback renders the steady state of the system more resistant in the sense that parameter perturbations and external stress conditions are unlikely to drive the system beyond the threshold where oscillations occur.
2.7. Application to the p53 system
The p53 system is one of the beststudied human signalling pathways, which is activated by various stress signals, including DNA damage [40,41]. Interestingly, p53 phosphorylation can exhibit both oscillatory behaviour and sustained activation, depending on the stimulus, which imply different cell fates [26,42]. A range of models have been developed for this pathway to understand dynamics and variability of the protein circuitry [24,25,43]. Here, we tested whether our modelling framework can also explain p53 and Mdm2 dynamics, possibly giving new insights into the feedback regulation circuitry of the system. To this end, we fitted again different model alternatives based on our general framework (2.1) to a published average oscillation pattern of p53 and Mdm2 dynamics after DNA damage [25] (figure 4b): 2.12where δ = (δ_{1}, δ_{2}, δ_{3}, δ_{4}, δ_{5}), , indicate model alternatives. Now, our model components are interpreted such that the signal x is p53 activation (sensing e.g. DNA damage), and the transducer y is an intermediate component, e.g. Mdm2 RNA. Consequently, for the latter component, no mass conservation is assumed, i.e. without the term (y_{t} − y) in (2.1). The response z represents Mdm2 protein concentration which, in turn, mediates p53 degradation. Like for the HOG model, we tested two kinetic alternatives for reactions f_{1}(x), and h(z), i.e. mass action (δ_{2}, δ_{5} = 0) or Michaelis–Menten kinetics (δ_{2}, δ_{5} = 1). In the p53 models, we assumed , because, assuming the transducer to be RNA, a fast autoinhibitory feedback seemed unlikely. Therefore, we tested autoinhibition in component z by alternatively introducing (δ_{4}), assuming the fast autoinhibitory feedback to act at the protein level by, e.g. posttranslational modifications. The kinetic nature of the negative feedback of Mdm2 on p53 remains elusive. Therefore, we also tested here two different alternatives for f_{3}(z), i.e. mass action (δ_{3} = 0) and Hilltype kinetics (δ_{3} = 1). Additionally, we tested, as for the HOG system, the possibility that the negative feedback acts by inputinhibition (δ_{1} = 1) or by outputactivation (δ_{1} = 0). Combination of the different possibilities results in 20 different models. The result of the fitting and ranking is displayed in table 2. For parameters of the bestranked model, please refer to the electronic supplementary material, table S3. A COPASI implementation of the bestranked model together with the data and an SBML version is also provided in the electronic supplementary material.
The two bestranked models (no. 11, δ = (0, 0, 1, 0, 1) and no. 15, δ = (0, 0, 1, 1, 1)) feature mass action kinetics in f_{1}(x), Hilltype kinetics in f_{3}(z), and zeroorder kinetics in h(z) and their fit is significantly better than for the other model candidates. Whether or not these two models have an autoinhibitory feedback does not influence the goodness fit itself (sum of squared residual (SSR) in table 2), but as the model without autoinhibitory feedback has two parameters less, it is clearly ranked first. The fit of the best approximating model no. 11 is shown in figure 8.
The p53 system can show both oscillatory as well as sustained behaviour, depending on the stimulus [26,42]. Therefore, we asked the question whether the oscillations of the best approximating model can be stabilized by introducing a fast autoinhibitory feedback. This is only true in general, i.e. irrespective of the other parameters, when f_{1}(x) has zeroorder kinetics at the equilibrium, which is not the case for the best approximating p53 model no. 11. However, with the given set of parameters, our fitted model can indeed be stabilized by the introduction of an autoinhibitory feedback of the form , with κ_{z} = 1.95 and m = 3 (figure 9).
Figure 9 depicts the stability of the steady state of the best approximating model with and without autoinhibitory feedback as a function of selected parameters (for additional pairs of parameters, refer to the electronic supplementary material, figure S10). Like for the HOG model, the unstable region diminishes through the introduction of a nested autoinhibitory feedback. Moreover, the stability of the system changes upon addition of the nested autoinhibitory feedback, rendering the system in the stable zone after introduction of an autoinhibitory feedback.
For the p53 system, a Monte Carlo analysis of the steady state with respect to parameter perturbations also indicated that the system with autoinhibitory feedback is less sensitive (figure 10).
The variance of the distance between the original and the perturbed steady states for all steady states is significantly smaller (p < 0.001) for the model with autoinhibitory feedback (model no. 15) compared with the model without autoinhibitory feedback (model no. 11; figure 10 and the electronic supplementary material, table S4).
Taken together, our simple framework suggests a mechanism how the p53 signalling system can change its dynamic behaviour upon different stimuli. Certain stimuli might activate components which introduce a nested autoinhibitory feedback. This changes the stability landscape of the system, shifting it from an oscillatory regime into a stable one. Thus, the p53 system depicts low resistance to parameter perturbations in order to be able to change its stability properties depending on environmental conditions.
3. Discussion
The ability to adapt to perturbations in external or internal conditions without losing structural stability is a fundamental feature of biological systems, including ecological, climate or biochemical systems. Adaptation is often mediated by negative feedbacks [1,7]. In biochemical systems, negative feedbacks inevitably come with time delays, which may lead to oscillatory behaviour both damped and sustained [15,23]. In some instances, oscillatory behaviour might oppose efficient adaptation owing to recurring stress. In such cases, the distance between the state of a system and the threshold beyond which oscillations occur, i.e. the systems resistance, should be large. This way, perturbations can be absorbed without affecting the structural stability of the system. In other instances, however, it might be beneficial to be able to switch between different dynamic regimes. It has been shown that the difference between a sustained or oscillatory signal can control cell fate [26,44], and that oscillation frequency can encode biochemical information [45]. In that case, resistance of a system should be low to be able to easily shift between different stability regimes. It might also be desirable to synthetically engineer cellular systems in a way such that oscillatory behaviour can be tuned by an independent artificial component.
For a threedimensional system, it has been observed that coupling autoinhibitory and delayed negative feedbacks reduces the probability of occurrence of stable limit cycles [5]. For the simple gene transcription network model with inputinhibition proposed by Goodwin [15], it has been shown that nested selfrepressing feedback loops have the potential to suppress oscillations [32]. Here, we propose a generic mechanism, how adaptive homeostatic biochemical systems can control both its dynamic response and its distance to a threshold beyond which these dynamics are drastically altered. We provide complete proof that in generic threedimensional homeostatic adaptive biochemical networks both with inputinhibition as well as with outputactivation oscillations may arise due to a Hopf bifurcation. We further prove that nested autoinhibitory feedbacks diminish the parameter space in which the steady state becomes unstable and oscillations occur. For systems with inputinhibition, the region for oscillations is generally reduced by autoinhibitory feedbacks. This is also true for models with a signalling module (mass conservation) and perfect adaptation (zeroorder h(z)). The latter renders inputinhibition systems susceptible to oscillations also for low cooperativity in the inputinhibition which extends former studies [16,32,46]. For our system with outputactivation, this is only true irrespective of the parameters, when the feedbackactivated output is of zeroorder. Thus, this is a sufficient, but not necessary condition. We show that this condition also has as the consequence that the system is stimulated by the positive difference between external and internal conditions. Our applications to the HOG and the p53 system suggest that zeroorder outputactivation might be a biologically more relevant feedback mechanism than inputinhibition. For the adaptive HOG system, we also demonstrate that a nested autoinhibitory feedback can alter the structural stability of the system without significantly affecting parameters of the system that are not involved in this feedback. Therefore, autoinhibition can alter stability properties of a system without affecting dynamic properties within a certain range of conditions. Owing to the generality of our model, we may hypothesize that other kinds of nested feedbacks also have the potential to suppress oscillatory behaviour. Here, we analysed only autoinhibitory feedback for mathematical convenience. However, in the special case of a Goodwintype model, it has been shown that a nested feedback from z to y also diminishes the parameter region in which oscillations occur [32]. In addition, for a fourdimensional model, it has even been demonstrated that a feedforward loop within an integral negative feedback also diminishes the parameter regions in which oscillations occur [21].
The application of our system to the HOG and the p53 system also provides evidence that nested autoinhibitory feedbacks increase the robustness of the steady state with respect to parameter perturbations. For the ERK pathway, it has been observed that a fast posttranslational feedback mechanism confers robustness to steadystate phosphorylation of ERK [47], which supports our analysis.
Our results may have implications to understand the complex dynamics of a range of signalling pathways. Not only has the p53 system been shown to exhibit different dynamics depending on the stimulus. The ERK pathway can show both oscillatory and adaptive dynamics, which are likely due to different feedback mechanisms that act on different timescales and that are activated depending on the stimulus [47–49]. The NFκB system can show damped oscillations, which are likely due to different feedback mechanisms acting on different timescales [27]. It seems that the coupling of fast posttranslational and delayed transcriptional feedbacks is a general feature of signalling pathways that allows finetuning of dynamics and steadystate features. The role of fast posttranslational negative feedbacks in this respect is apparently either to suppress oscillatory behaviour or stabilize steadystate protein levels or both.
The presented theoretical results on suppressing oscillatory behaviour induced by Hopf bifurcations may be useful in designing synthetic systems in which oscillations can be tuned by synthetic autoinhibitory feedbacks. This may be useful for studying cell fate decisions, as, for example, in the p53 or the ERK system. For the HOG system, the parametrized models show that even without autoinhibitory feedback osmoadaptation is extremely stable. For this system, it seems unlikely that oscillations can be induced artificially by weakening the reported autoinhibitory feedbacks.
4. Material and methods
4.1. Data
We made extensive use of published data to parametrize dynamic models of the HOG pathway and the p53 pathway. The dataset used for model parametrization and discrimination of the HOG model was taken from [17]. This dataset consists of time series of phosphorylated Hog1 under several hyperosmotic shock conditions, for wildtype and different mutants yeast, for up to 2 h after hyperosmotic shock (figure 5a,b). Additionally, we used a time series of glycerol published in [22] (figure 5c). These datasets, although coming from different sources, are comparable because they were produced using the same genetic background and under the same culture conditions. The dataset used for model parametrization and discrimination of the p53 model was digitized from the electronic supplementary material, figure S6 of the supplementary material of [25]. These data are meant to resemble an idealized undamped oscillation with peak characteristic that correspond to the average peak characteristic of oscillating cells. For the ranking procedure, we considered only 91 data points, because the last three periods in figure 8 are repetitions of the former oscillations.
4.2. Model fitting, ranking and selecting
The models were implemented and fitted with the free software COPASI (v. 4.7, build 34) [50]. We used the evolutionary programming algorithm to fit the models, where the population size was set to 10 times the number of parameters and the number of generations was limited to 10 times the number of parameters. When estimated parameters hit parameter boundaries, the boundaries were relaxed and the model refitted until the fit converged within defined parameter boundaries. Model ranking was performed using modelMaGe [51,52]. For model ranking, we calculated the Akaike information criterion corrected (AICc) for small sample sizes [53] for each candidate model: where SSR is the sum of squared residuals of the fit, k is the number of parameters and n is the number of data points. The AICc is an informationtheorybased measure of parsimonious data representation that incorporates the goodness of the fit (SSR) as well as the complexity of the model (k), thereby giving an objective measure for model selection and discrimination.
In order to select and compare the best approximating model(s), we calculated the Akaike weights (AICw) [53] where Δ_{i} = AIC_{i} − AIC_{min}, with AIC_{i} being the AICc for model i, i = 1,…, R according to ranking and AIC_{min} the minimal AICc. The AICws can be considered as the weight of evidence in favour of a model given as a number between 0 and 1, i.e. the higher the weight, the closer the model is to the hypothetical true model [53]. We considered those models as best approximating that had an .
Funding statement
This study was supported by the German Ministry of Science and Education (BMBF project no. 0135779 and 0316188E to J.S.) and the International Max Planck Research School Magdeburg for Advanced Methods in Process and Systems Engineering. The authors declare no conflicting interests.
 Received October 21, 2013.
 Accepted November 12, 2013.
© 2013 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0/, which permits unrestricted use, provided the original author and source are credited.