 Research Article
 Open Access
 Published:
How time delay and network design shape response patterns in biochemical negative feedback systems
BMC Systems Biology volume 10, Article number: 82 (2016)
Abstract
Background
Negative feedback in combination with time delay can bring about both sustained oscillations and adaptive behaviour in cellular networks. Here, we study which design features of systems with delayed negative feedback shape characteristic response patterns with special emphasis on the role of time delay. To this end, we analyse generic twodimensional delay differential equations describing the dynamics of biochemical signalresponse networks.
Results
We investigate the influence of several design features on the stability of the model equilibrium, i.e., presence of autoinhibition and/or mass conservation and the kind and/or strength of the delayed negative feedback. We show that autoinhibition and mass conservation have a stabilizing effect, whereas increasing abruptness and decreasing feedback threshold have a destabilizing effect on the model equilibrium. Moreover, applying our theoretical analysis to the mammalian p53 system we show that an autoinhibitory feedback can decouple period and amplitude of an oscillatory response, whereas the delayed feedback can not.
Conclusions
Our theoretical framework provides insight into how time delay and design features of biochemical networks act together to elicit specific characteristic response patterns. Such insight is useful for constructing synthetic networks and controlling their behaviour in response to external stimulation.
Background
Negative feedback is one of fundamental mechanisms in cellular networks [1–7], which fulfils a variety of functions such as mediating adaptation [8–10], stabilizing the abundance of biochemical components [1, 5, 11, 12], inducing oscillations [7, 13–15] and decoupling signal and response time [5]. Negative feedbacks are shown to be present in many biochemical systems including bacterial adaptation [9, 16], mammalian cell cycle [17, 18], stress response in mammalian cells [19] and yeast [20, 21].
Negative feedbacks may involve a time delay, which is needed for signal transduction and transcription, translation and formation of biochemical species [21–24]. The combination of negative feedback and time delay may result in oscillatory dynamics of components of the cellular network [25–27]. Oscillations induced by delayed negative feedbacks (DNFs) were experimentally observed in several biochemical systems as a response to external stimuli and stress, e.g., mammalian Hes1 [24, 28], p53 [29–31] and NF κB [23, 32, 33] systems.
It is conceivable that oscillatory behaviour might be inappropriate in biological systems mediating adaptive responses. For example, in the hyperthermia treatment of cancer, largeamplitude temperature oscillation could result in tissue damage or patient discomfort [34]. We wondered, if there exist design features and mechanisms of systems containing DNF, which may suppress oscillatory behaviour caused by external stimulation. Recent studies [22, 35] demonstrated that nested negative feedbacks may suppress oscillations of biochemical species involved in DNF. However, these studies provided no insight into how time delay influences the dynamics of DNF systems and interacts with nested negative feedbacks.
In our previous study [36] we derived explicit thresholds and boundaries showing how time delay determines characteristic response patterns of biochemical networks containing DNF. In this manuscript, we continued our research and investigated how the combination of time delay and certain design features influences the dynamics of biochemical DNF systems. To this end, we constructed a range of generic twodimensional DNF models using delay differential equations. Models differed in several properties:

(i)
presence of a nested negative (autoinhibitory) feedback,

(ii)
presence of mass conservation for biochemical compounds,

(iii)
mechanism of DNF, i.e., inputinhibition or outputactivation.
Further, we subjected these models to computational and analytical stability analyses with respect to time delay. Our computational analyses demonstrate that

the presence of autoinhibition and mass conservation have a stabilizing influence on the model equilibrium independent of the DNF strength.

increasing abruptness and decreasing DNF threshold have a destabilizing effect on the model equilibrium.
Our theoretical analyses show that

nested autoinhibitory feedbacks may increase the range of time delay, where the system is stable, through the abruptness of the feedback function.
Applying our theoretical framework to the oscillating p53 system in mammalian cells [37] indicates that

both period and amplitude of p53 oscillations increase with time delay, and

a nested autoinhibitory feedback can decouple period and amplitude of oscillations, whereas the delayed feedback can not.
Our analysis provides insight into how time delay and specific design features act in concert to shape the systems doseresponse relationship. This knowledge can be used for constructing synthetic networks with the finetuned behaviour.
Methods
Data
The dataset used for the parametrization of the p53 model was digitized from the supplementary material of [30] from Additional file 1: Figure S6 as described in [22]. It represents an averaged oscillation pattern that was meant to resemble an idealized undamped oscillation with peak characteristics that correspond to the average peak characteristics of oscillating cells.
Model simulation and analysis
All simulations of the delay differential equations were carried out in Mathematica 9 (Wolfram Research, Champaign, Illinois) using the function NDSolve based on the method of steps.
We used DDEBIFTOOL v. 2.00 [38] and MATLAB R2008b (The MathWorks, Natick, MA) to calculate dependencies between the value of time delay τ and amplitude and period of oscillations of the p53 model.
MonteCarlo analysis was performed in Mathematica 9.
For the parameter estimation we used the leastsquares method minimizing the sum of squared residuals (SSR):
where p=(p _{1},p _{2},…,p _{ m }) denotes a set of parameters to be estimated, x(t,p) is a numerical solution depending on parameters p, x _{ i } is a measured data point at the time t _{ i }, n is a number of the data points.
For minimizing S S R(p) with respect to parameter values we utilized the numerical function NMinimize in Mathematica 9, which, by default, uses the “NelderMead” method. In case “NelderMead" performs poorly, it automatically switches to the “Differential Evolution” method. The parameter optimization process is assumed to have converged to a local minimum if the difference between the new best and the old best function value S S R(p), as well as the distance between the new best and the old best parameter values, are less than a tolerance of 10^{−8}.
Robustness with respect to model parameters
We analysed the robustness of the parameter fit for the p53 model with respect to noise. Parameter values were randomly sampled within ±10 % of their respective fitted values using a uniform distribution for 100 times. Then the p53 model was simulated using perturbed parameter sets. The relative variation of the integral of the first transient response after the initial stimulus was calculated (Additional file 1: Figure S6). Namely, we calculated the integral of the simulated p53 model from the time point 0 until the time of the first minimum after stimulation. This way, the robustness of both initial activation amplitude and timing of the first transient response, two characteristic measures of system dynamics, have been estimated concomitantly.
Results and discussion
Model formulation
We investigated four different twodimensional models containing DNF that describe generic signalresponse networks (Fig. 1). Models differ in design of the DNF and presence of mass conservation for a biochemical compound. All models can have a nested negative (autoinhibitory) feedback.
We mathematically formulated models from Fig. 1 using deterministic delay differential equations: Model 1: DNF with inputinhibition and without mass conservation (Fig. 1 a):
Model 2: DNF with inputinhibition and with mass conservation (Fig. 1 b):
Model 3: DNF with outputactivation and without mass conservation (Fig. 1 c):
Model 4: DNF with outputactivation and with mass conservation (Fig. 1 d):
The function \(S_{1}\colon [0,\infty)\to \mathbb {R}_{+}\) is twice continuously differentiable and monotonically decreasing with R. The function \(S_{2}\colon [0,\infty)\to \mathbb {R}_{+}\) is twice continuously differentiable and monotonically increasing with R. The twice continuously differentiable monotonically decreasing function F:[0,∞)→(0,1] mimics an optional autoinhibitory feedback.
Parameters of Models 14 have real positive constant values. For convenience, we combined them into the vector p:
Note that all model parameters represent lumped biological processes and therefore have only limited physical or biological meaning. For Models 1, 2 the parameter δ equals to 0.
In all models the input I mimics some constant stimulus (e.g., radiation, see below) that activates a gene transcription network (Fig. 1 a, c) or a signalling cascade (Fig. 1 b, d) starting with the component C. Further, the component C activates the response variable R with a certain time delay τ caused by processes like transport, transcription, translation, etc. Further, the response R negatively feeds back through S _{1} or S _{2} performing the DNF. Depending on the model design the response R mediates DNF using different mechanisms.
We refer to models from Fig. 1 a, c as a transcription network, because C is not reversibly converted into different states, but rather produced and degraded. In the model from Fig. 1 a the response variable R has the inhibiting influence on C by means of the function S _{1}(R). Therefore, we refer to this DNF mechanism as inputinhibition. In the model from Fig. 1 c the DNF mediated by the response R acts through activating the degradation of C by means of the function S _{2}(R). We refer to this DNF mechanism as outputactivation. We used this model to simulate p53 protein dynamics [30, 37].
We considered both inputinhibition and outputactivation architectures together with a socalled signalling component (see Fig. 1 b, d, respectively). The component C activating the response R originates from another component \(\breve {C}\) to which it is constitutively converted back. The sum of both components \(C_{t}=C+\breve {C}\) is a constant and assumed to be unity (C _{ t }=1). We refer to this model feature as mass conservation. This modelling technique mimics a fast and reversible posttranslational protein modification, e.g., phosphorylation, leaving the total protein content unchanged, as it is often described in signalling cascades.
In the following sections we presented theoretical and computational analyses of Models 14 with application to p53 system in mammalian cells.
Stability analysis of systems with delayed negative feedback
Presented in this section stability analysis can be applied to all Models 14. Therefore, we do not make an explicit distinction between models, unless necessary.
As the equilibria of Models 14 we considered the vector E=(C _{ s },R _{ s })^{T}. The equilibria E always exist (see Additional file 1) and implicitly depend on the parameter vector p (5). Model equilibria E also depend on the input I and, therefore, Models 14 are not able to show a perfect adaptation.
We linearised Models 14 about their respective equilibria E:
where for Model 1, we have
for Model 2, we have
for Model 3, we have
for Model 4, we have
The analysis of the model (6) revealed the following stability properties of the equilibrium E with respect to x, y, β, τ (for details refer to Additional file 1):

If x β≥y holds, then the equilibrium E is absolutely stable, i.e., stable for any τ≥0.

If x β<y holds, then there exists a marginal time delay τ _{ m } (the Hopf bifurcation point) such that the equilibrium of the model (6) is stable for any 0<τ<τ _{ m } and unstable for any τ≥τ _{ m }. The marginal time delay τ _{ m } is calculated as a product of functions f and g that depend on β and x, y from (7)–(10), respectively:
$$ \tau_{m}(x,y,\beta)=f(x,y,\beta)\times g(x,y,\beta), $$(11)where
$$ \begin{aligned} &f(x,y,\beta)\\&= \frac{\sqrt{2}}{\sqrt{x^{2}\beta^{2}+\sqrt{(x^{2}+\beta^{2})^{2}+4(y^{2}x^{2} \beta^{2})}}}>0,\\ &g(x,y,\beta)\\&=\arccos{\frac{(x+\beta)^{2}+\sqrt{(x^{2}+\beta^{2})^{2}+4(y^{2}x^{2} \beta^{2})}}{2 y}}>0. \end{aligned} $$(12)The derivation of functions f and g is represented in Additional file 1.
In the next section, we considered features and mechanisms of systems with DNF that may stabilize the system equilibrium after stimulation.
Design features stabilizing biochemical delayed negative feedback systems
Recently, two research articles indicated that nested autoinhibitory feedbacks repress oscillatory dynamics of simple biochemical networks involving a nonlinear DNF [21, 35]. We wondered how other design features of systems with DNF influence the stability of the model equilibrium. Thus, in addition to autoinhibition we investigated the influence of the following designs on the stability of the model equilibrium:

Mechanism of DNF: inputinhibition or outputactivation,

Strength of DNF: strong or weak,

Presence of mass conservation for a chemical compound.
We also considered how different combinations of delayed and autoinhibitory negative feedbacks affect the stability of the equilibrium:

Weak DNF with and without autoinhibition,

Strong DNF with and without autoinhibition.
For analysing the influence of these design features on the stability of the equilibrium we performed a MonteCarlo analysis of Models 14. First, we defined concrete DNF functions S _{1} and S _{2} and an autoinhibitory feedback function F.
We defined a reverse Hill function as the inputinhibition DNF function S _{1}(R). As the outputactivation function S _{2}(R) we defined a Hill function. Thus, functions S _{1} and S _{2} have the following form:
with K _{ m }>0 being the halfsaturation constant, characterizing the activation threshold beyond which the feedback takes effect, and n≥1 being the Hill coefficient, characterizing how abrupt the feedback takes effect after having passed the activation threshold. Thus, parameters K _{ m } and n specify the strength of the DNF: the lower the activation threshold K _{ m } and the higher the abruptness n, the stronger the DNF is. Note that applying the same parameters make functions S _{1} and S _{2} symmetric, allowing a fair comparison of the influence of inputinhibition and outputactivation on the model stability.
As the autoinhibitory feedback we employed a reverse Hill function F(C):
Then, we randomly generated parameter values I=0.87, α=0.11, β=0.17, δ=58.2, n=12.77, K _{ m }=0.23 in the way that Models 14 without autoinhibition, i.e., κ=0, have similar values of τ _{ m }. They correspond to τ _{ m }=1.23, τ _{ m }=1.27, τ _{ m }=1.22, τ _{ m }=1.24 for Models 14, respectively. Thus, applying this parameter set we guarantee that any value of the time delay τ has a similar distance to the Hopf bifurcation point τ _{ m } for all considered models. Simulations of Models 14 with these parameters and τ=2.5 are depicted in Fig. 2 a.
Further, using these parameter values we performed a MonteCarlo analysis for Models 14 with the constant input I=0.87. Namely, the rate constants α, β, δ have been sampled in the range from 0.1 to 10 times of their respective values for 10000 times. The parameter values defining the system design, i.e., n, K _{ m }, ν and κ, were sampled in the following way:

(i)
in the case of weak DNF without autoinhibition (κ=0) we sampled n in the range from 0.1 to 1 time of its value and K _{ m } in the range from 10 to 20 times of its value 10000 times.

(ii)
in the case of strong DNF without autoinhibition (κ=0) we sampled n in the range from 1 to 2 times of its value and K _{ m } in the range from 0.1 to 10 times of its value 10000 times.

(iii)
in the case of weak DNF with autoinhibition we sampled n and K _{ m } as in (i), κ in the interval [0.1,10], ν in the interval [1,20] 10000 times.

(iv)
in the case of strong DNF with autoinhibition we sampled n and K _{ m } as in (ii), κ in the interval [0.1,10], ν in the interval [1,20] 10000 times.
Thus, for each Model 14 and each combination of DNF and autoinhibition we obtained 10000 parameter sets. For each parameter set we calculated x and y according to (7)–(10) for Models 14, respectively. Then, we defined the percentage of parameter sets for which x β≥y holds indicating absolute stability of the model equilibrium (see Fig. 2 b). For all other parameter sets, we calculated the mean value of the marginal time delay τ _{ m } (11), i.e., the Hopfbifurcation point (see Fig. 2 c).
Figure 2 b, c shows that considered models with autoinhibition have a higher percentage of parameter sets leading to absolute stability and higher mean value of τ _{ m } than the same models without autoinhibition. This observation confirms previous results [21, 35] showing that nested autoinhibitory feedbacks repress oscillatory dynamics in networks containing DNF.
Additionally, for models with weak DNF there are more parameter sets, which induce absolute stability of the model equilibrium, than for models with the strong DNF. Accordingly, models with weak DNF have higher mean value of τ _{ m }, i.e., are less prone to oscillations, than models with the strong one. Thus, the nested autoinhibitory feedback and the DNF have opposing effects on the system’s stability. The stronger the autoinhibitory feedback and the weaker the DNF, the less probable an oscillatory response of the system is.
Moreover, Models 2 and 4 with mass conservation have a higher percentage of parameter sets leading to absolute stability of the model equilibrium and higher mean value of τ _{ m } than respective Models 1 and 3 without mass conservation. In order to explain this effect, for each model we quantified the dependence between τ _{ m } and each parameter value in the range from 0.1 to 10 times of its default value leaving the other parameters fixed (see Additional file 1: Figure S11a–d). This sensitivity analysis shows that the presence of mass conservation influences the sensitivity of τ _{ m } only with respect to the halfsaturation rate K _{ m } leaving all other parameter sensitivities unchanged (compare Additional file 1: Figure S11a to b and S11c to d for models without and with mass activation, respectively, and Additional file 1: Figure S11e). In the presence of mass conservation τ _{ m } increases much stronger with increasing feedback activation threshold K _{ m } and, therefore, stabilizes the equilibrium. Moreover, in the presence of mass conservation the value of K _{ m } beyond which τ _{ m } does not exist any more, i.e., the system’s equilibrium becomes absolutely stable, also decreases (Additional file 1: Figure S11e).
Concerning the design of the DNF, MonteCarlo analysis shows that Models 1, 2 with inputinhibition have a higher percentage of parameter sets leading to absolute stability and a higher mean value of τ _{ m } than Models 3, 4 with outputactivation (see Fig. 2 b, c). However, we were not able to support these results with an alternative parameter set I=0.48, α=0.14, β=0.44, δ=83.71, n=10, K _{ m }=0.9. Refer to Additional file 1: Figure S12a for simulations of Models 14 using the alternative parameter set and τ=10. In comparison, for the alternative parameter set Models 1, 2 with inputinhibition have higher values of τ _{ m }, i.e., τ _{ m }=1.27 and τ _{ m }=1.86, than Models 3, 4 with outputactivation, i.e., τ _{ m }=0.36 and τ _{ m }=0.52. Thus, the distance of τ to τ _{ m } is smaller for Models 1, 2, and one may expect a higher percentage of parameter sets inducing absolute stability. Nevertheless, for the alternative parameter set the MonteCarlo analysis showed that models with inputinhibition have approximately the same percentage of parameter sets leading to absolute stability as corresponding models with outputactivation. In comparison, all other conclusions presented above were confirmed for model simulations with the alternative parameter set (Additional file 1: Figure S12b).
Taken together, we conclude that autoinhibition as well as mass conservation have a stabilizing influence on the model equilibrium independent of the strength of DNF and allow systems with DNF to adapt to an external stimulus without producing sustained oscillations. Moreover, the higher the activation threshold and the less abrupt the DNF, the less prone the system is to an oscillatory behaviour.
Autoinhibition increases τ _{ m }
Computational analysis presented in the above section demonstrated the opposing behaviour of autoinhibitory and delayed negative feedbacks with respect to stability (Fig. 2). In this section, we analytically investigated how the autoinhibitory feedback stabilizes the equilibrium of the system.
We proved that τ _{ m } (12) increases with x and decreases with y (see Additional file 1):
Further, we derived upper and lower bounds for x and y from (7)–(10) for Models 14, respectively (see Additional file 1):
Both the lower and upper bound of x, i.e., ε _{ lb } and ε _{ ub }, are increasing with F ^{′}(C _{ s }) for Models 14. Therefore, we can always increase a given x by choosing an appropriate value of F ^{′}(C _{ s }). The lower and upper bound of y, i.e., σ _{ lb } and σ _{ ub }, have nonnegative constant values. Consequently, according to (16), we can always increase τ _{ m } by increasing F ^{′}(C _{ s }).
Taken together, we showed that autoinhibitory feedback can increase the range of the interval [0,τ _{ m }) and stabilise the model equilibrium.
Application to the p53 system
In this section we applied Model 3 to the p53 system to gain novel insights into the functioning of this system.
The tumour suppressor protein p53 is activated in response to many stress signals and activates various stressresponse programs including cellcycle arrest, senescence and apoptosis [39]. It is also well established that p53 acts within a negative feedback loop, including Mdm2 as the negative regulator of p53: p53 transcriptionally activates Mdm2, which in turn targets p53 for degradation [29, 40].
Several mathematical models of p53Mdm2 feedback loop have been published [22, 28, 30, 37]. One of these models (model III from Table 1 in [30]) is a particular case of the model from Fig. 1 c corresponding to the Model 3 with F(C)≡1. Therefore, we wondered, whether our framework would also be able to explain measured p53 dynamics upon DNA damage. In our designations C and R correspond to p53 and Mdm2, respectively. The input I is defined here as a scaled DNA damage signal and is measured in arbitrary units. The negative feedback by outputactivation is modelled by a nonlinear Hill function S _{2}(R) (14).
We fitted parameters of the p53 model (3) to the experimental data of an averaged oscillation pattern of the p53Mdm2 system after DNA damage from [30] (Additional file 1: Fig. S6 therein). The results of the fit are presented in Additional file 1: Table S1. Figure 3a shows the simulation of the p53 model (3) with fitted parameters. The model well recapitulates measured p53 dynamics after DNA damage. Moreover, the fitted optimal solution is also robust with respect to noise in the fitted parameters (see Additional file 1). Indeed, the integrated first transient response varies only by 8.7 % assuming a parameter noise of ±10 % (see Methods Section and Additional file 1).
The model analysis shows that the fitted time delay τ=1.37 h is almost two times larger than the corresponding τ _{ m }=0.76 h that was calculated for the fitted DNA damage signal I=0.23 (Fig. 3 b). Therefore, the p53 model (3) with fitted parameters from Table S1 (see Additional file 1) produces sustained oscillatory response.
It was earlier reported that distinct p53 dynamics such as oscillations or sustained activation may lead to different cell fate decisions [31, 39]. Recent study [41] showed that the system’s response is modulated by DNA damage strength. Namely, after high DNA damage p53 level was monotonically increased and cells activated apoptosis, whereas after low DNA damage p53 level underwent periodic pulsing resulting in a cellcycle arrest. We checked if our generic p53 model (3) is able to reproduce this transition with respect to the DNA damage level. Figure 3 b shows that τ _{ m } is decreasing with respect to the DNA damage signal I. Hence the fitted value of time delay τ is greater than τ _{ m } for any I>0.23 (fitted value). Therefore, the p53 model (3) produces sustained oscillations for any I greater than the fitted value and is not able to perform the transition from oscillatory to adaptive behaviour with respect to increased DNA damage signal I. This conclusion is also applicable to other model alternatives (1)–(4) used in our previous study presenting models of HOG pathway in yeast and NF κB signalling in mammalian cells (see Figs. 3b and 6b in [36]). However, according to [41] the switch from sustained oscillations to monotonic increase of p53 level is regulated by a mechanism attenuating Mdm2 expression that is not present in the current p53 model. Studies [39, 41] considered DNA damage kinases ATM and ATR as negative regulators of Mdm2 expression. Using this knowledge, we extended the p53 model (3) by including the additional component ATM activating p53 and attenuating Mdm2 (see Additional file 1: Figure S7). As a result, the extended p53 model was able to qualitatively reproduce the switch from oscillations to monotonic increase of p53 level (see Additional file 1: Figure S8). Simulations of the extended p53 model suggest that this transition between response types originates from the competition between ATM and p53 for the inhibition and activation of Mdm2, respectively. In case of high DNA damage, ATM level is high and suppresses Mdm2 giving a monotone increase of p53 level. In case of low DNA damage, Mdm2 activity is not effectively suppressed by ATM resulting in sustained oscillations of both p53 and Mdm2.
Further, we applied our theoretical analysis to explore under what conditions sustained oscillations of p53 model (3) can be suppressed by the activation of a nested autoinhibitory feedback to the model component C preserving all values of fitted parameters. Our theoretical analysis suggested that the marginal time delay τ _{ m }, beyond which any time delay leads to sustained oscillations, can be increased by increasing the slope of the autoinhibitory feedback function at the equilibrium F ^{′}(C _{ s }). As in the previous section, as the autoinhibitory feedback function we utilized a reverse Hillfunction F(C) (15). Further, we adjusted parameters ν and κ of F(C) and calculated τ _{ m } (see Additional file 1). For ν=3 the marginal value of time delay τ _{ m } is larger than the fitted time delay τ=1.37 h (Fig. 3 b). As a result, p53 model (3) with parameters from Table S1 in Additional file 1 produced damped oscillatory response (Additional file 1: Figure S4).
In a similar DNF system it was shown that the period of oscillations increases with the Hill coefficient n of the DNF function for a given time delay [24]. We wondered how parameters of the delayed negative and autoinhibitory feedbacks influence the amplitude and period of oscillations in our system. Figure 4 demonstrates that the autoinhibitory feedback (with parameters ν=3, κ=1.73) decreases and stabilizes the amplitude of oscillations, whereas the amplitude of oscillations increases with respect to the Hill coefficient n of the DNF function S _{2}. Moreover, increasing the abruptness of the DNF has no substantial influence on the increase of period with respect to time delay τ. The period of oscillations is a linear function of the time delay τ irrespective of values of ν, κ and n. Thus, opposed to the delayed feedback, the autoinhibitory feedback has the potential to decouple the increase of amplitude and period of oscillations with respect to τ. Moreover, autoinhibitory and delayed negative feedbacks have an opposing influence on the amplitude of oscillations.
Thus, our analysis showed that for the p53 model (3) an autoinhibitory feedback can be a potential mechanism increasing the marginal time delay τ _{ m }, decreasing the amplitude of oscillations and turning sustained oscillations into damped oscillations.
The experimental study of p53 oscillations [29] concluded that the mean number of p53 pulses in individual cells increased with DNA damage. Moreover, the authors suggest that the p53Mdm2 feedback loop generates a “digital” clock making the number of p53 pulses relevant for the cell fate, and not their amplitude and duration. However, this hypothesis has not been proven yet. Therefore, we wondered which parameters of p53 model (3) play a prominent role in controlling the length of p53 pulses. Using p53 model (3) we split the p53 simulation curve on “On” and “Off” states (see Additional file 1: Figure S9). Then we checked how different model parameters control the duration of p53 pulses (see Additional file 1: Figure S10). The analysis showed that time delay τ is the only parameter that significantly changed the duration of “On” and “Off” states of the model response (see Additional file 1: Figure S10b). Namely, time delay τ increases the duration of “Off” states and decreases the duration of “On" states. In addition, τ increases the amplitude and period of pulses (see Fig. 4). Note that the same conclusion can be applied to the relation between time delay τ and marginal time delay τ _{ m }: the higher τ/τ _{ m }, the higher the amplitude and period of pulses are. It would be interesting to validate these predictions experimentally and check the physiological effect of changing time delay between p53 and Mdm2 activation after DNA damage.
Conclusions
Negative feedback in combination with time delay can induce oscillations in cellular networks [25–27]. However, oscillations might be inappropriate in biological systems with adaptive behaviour [34].
Here, we systematically study how design features in combination with time delay tune the response patterns of biochemical networks. To this end, we create a range of models containing an explicit time delay and a DNF differing in several aspects: presence of a nested negative (autoinhibitory) feedback, presence of mass conservation for a system component and mechanism of DNF, i.e., inputinhibition or outputactivation (Fig. 1). The obtained models (Models 14) are mathematically described by twodimensional delay differential Eqs. 1–(4) and subjected to computational and theoretical stability analyses with respect to time delay. The general idea to specifically address the interaction of time delay and design features was that all these design features act on the system stability and response pattern by modifying the time delay threshold, i.e., the bifurcation point, beyond which the system’s stability properties change.
We show that

nested autoinhibitory feedbacks and overall delayed negative feedbacks have opposing roles with respect to the characteristic response pattern. Indeed, nested autoinhibitory feedbacks have the potential to suppress oscillatory behaviour, whereas the increasing strength of the DNF promotes oscillations. Moreover, in oscillatory systems autoinhibitory feedbacks decouple amplitude and period of oscillations.

mass conservation has a stabilizing effect on the system’s equilibrium.

depending on the parameter set, the type of DNF can also influence the response pattern. We found that inputinhibition can be more stabilizing compared to outputactivation.
Thus, biochemical networks have a range of design possibilities shaping both their dynamic as well as their equilibrium properties. Our systematic analysis of different design features allows predicting what kind of biochemical network underlies a certain characteristic response. For example, in oscillatory systems with a long time delay, it is reasonable to assume a limited number of posttranslational modifications (mass conservation), no nested feedbacks and a strong overall negative feedback. Whereas adaptive systems with long time delay are likely to harbour nested negative feedbacks and posttranslational modifications. Systems with low number of components and short time delay that are meant to oscillate, will need an abrupt negative feedback with low activation threshold, whereas short time delay and a weak negative feedback are good designs principles for adaptive systems.
Our framework of delayed and nondelayed feedbacks can serve to support a design process for novel synthetic generegulatory networks. Indeed, our study allows to approximate the value of time delay and the structure of the DNF system for obtaining a certain type of the system dynamics. As an example, we considered p53 system in mammalian cells that contains DNF and is able to produce both oscillatory and adaptive responses depending on the source and strength of DNA damage [37, 41]. Although many studies are dedicated to studying DNA damage response in cells, the purpose of oscillations in p53 system remains unclear [29, 37, 39, 41]. The earlier study [37] suggested that oscillatory behaviour can be advantageous for the p53 system to achieve a tradeoff between irreversible biological outcome, e.g., irreversible cell cycle arrest or apoptosis, and insufficiently low levels of p53. Thus, oscillations have been viewed as repetitive repair efforts allowing the system to check after every p53 pulse whether the damage has been properly repaired. Our analysis showed that time delay increases the duration of “Off” states and decreases the duration of “On” states. Additionally, time delay may increase the amplitude and period of oscillations. According to our analysis the autoinhibitory feedback is able to decouple the amplitude and period of oscillations with respect to time delay. Thus, our study suggests that autoinhibition and time delay may control oscillations in p53 system. The experimental validation of these predictions may help to better understand the role of p53 oscillations and indicate more efficient treatment of diseases caused by violation of p53 regulation.
Abbreviations
 DNF:

delayed negative feedback
 ODE:

ordinary differential equation
 SSR:

sum of squared residuals
References
Alon U, Vol. 10. An Introduction to Systems Biology: Design Principles of Biological Circuits. Boca Raton: Chapman & Hall/CRC; 2006.
Alon U. Network motifs: theory and experimental approaches. Nat Rev Genet. 2007; 8(6):450–61.
Legewie S, Herzel H, Westerhoff HV, Blüthgen N. Recurrent design patterns in the feedback regulation of the mammalian signalling network.Mol Syst Biol. 2008; 4:190.
Kholodenko BN, Hancock JF, Kolch W. Signalling ballet in space and time. Nat Rev Mol Cell Biol. 2010; 11(6):414–26.
Tyson JJ, Chen KC, Novak B. Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell. Curr Opin Cell Biol. 2003; 15(2):221–31.
Sauro HHM, Kholodenko BBN. Quantitative analysis of signaling networks. Prog Biophys Mol Biol. 2004; 86(1):5–43.
Tsai TYC, Choi YS, Ma W, Pomerening JR, Tang C, Ferrell JE. Robust, tunable biological oscillations from interlinked positive and negative feedback loops. Science. 2008; 321(5885):126–9.
Ma W, Trusina A, ElSamad H, Lim WA, Tang C. Defining network topologies that can achieve biochemical adaptation. Cell. 2009; 138(4):760–73.
Yi TM, Huang Y, Simon MI, Doyle J. Robust perfect adaptation in bacterial chemotaxis through integral feedback control. Proc Natl Acad Sci. 2000; 97(9):4649–653.
Ni XY, Drengstig T, Ruoff P. The control of the controller: molecular mechanisms for robust perfect adaptation and temperature compensation. Biophys J. 2009; 97(5):1244–53.
Hasty J, Dolnik M, Rottschäfer V, Collins JJ. Synthetic gene network for entraining and amplifying cellular oscillations. Phys Rev Lett. 2002; 88(14):148101.
Sturm OE, Orton R, Grindlay J, Birtwistle M, Vyshemirsky V, Gilbert D, Calder M, Pitt A, Kholodenko B, Kolch W. The mammalian MAPK/ERK pathway exhibits properties of a negative feedback amplifier. Sci Signal. 2010; 3(153):90.
Novak B, Tyson JJ, Gyorffy B, CsikaszNagy A. Irreversible cellcycle transitions are due to systemslevel feedback. Nat Cell Biol. 2007; 9(7):724–8.
Elowitz MB, Leibler S. A synthetic oscillatory network of transcriptional regulators. Nature. 2000; 403(6767):335–8.
Kholodenko BN. Negative feedback and ultrasensitivity can bring about oscillations in the mitogenactivated protein kinase cascades. Eur J Biochem. 2000; 267(6):1583–1588.
Kollmann M, Løvdok L, Bartholomé K, Timmer J, Sourjik V. Design principles of a bacterial signalling network. Nature. 2005; 438(7067):504–7.
Novak B, Kapuy O, DomingoSananes MR, Tyson JJ. Regulated protein kinases and phosphatases in cell cycle decisions. Curr Opin Cell Biol. 2010; 22(6):801–8.
Ferrell JE, Tsai TYC, Yang Q. Modeling the cell cycle: why do certain circuits oscillate?. Cell. 2011; 144(6):874–85.
Blüthgen N. Transcriptional feedbacks in mammalian signal transduction pathways facilitate rapid and reliable protein induction. Mol Biosyst. 2010; 6(7):1277–1284.
Klipp E, Nordlander B, Krüger R, Gennemark P, Hohmann S. Integrative model of the response of yeast to osmotic shock. Nat Biotechnol. 2005; 23(8):975–82.
Schaber J, Baltanas R, Bush A, Klipp E, ColmanLerner A. Modelling reveals novel roles of two parallel signalling pathways and homeostatic feedbacks in yeast. Mol Syst Biol. 2012; 8:622.
Schaber J, Lapytsko A, Flockerzi D. Nested autoinhibitory feedbacks alter the resistance of homeostatic adaptive biochemical networks. J R Soc Interface. 2014; 11(91):20130971.
Hoffmann A, Levchenko A, Scott ML, Baltimore D. The I κBNF κB signaling module: Temporal control and selective gene activation. Science. 2002; 298(5596):1241–1245.
Bernard S, Cajavec B, PujoMenjouet L, Mackey MC, Herzel H. Modelling transcriptional feedback loops: the role of Gro/TLE1 in Hes1 oscillations. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences. 2006; 364(1842):1155–1170.
Goodwin BC. Temporal Organization in Cells: a Dynamic Theory of Cellular Control Processes. London: Academic Press; 1963.
Griffith JS. Mathematics of cellular control processes i. negative feedback to one gene. J Theor Biol. 1968; 20(2):202–8.
Mahaffy J, Pao C. Models of genetic control by repression with time delays and spatial effects. J Math Biol. 1984; 20(1):39–57.
Monk NAM. Oscillatory expression of hes1, p53, and NF κB driven by transcriptional time delays. Curr Biol. 2003; 13(16):1409–1413.
Lahav G, Rosenfeld N, Sigal A, GevaZatorsky N, Levine AJ, Elowitz MB, Alon U. Dynamics of the p53Mdm2 feedback loop in individual cells. Nat Genet. 2004; 36(2):147–50.
GevaZatorsky N, Rosenfeld N, Itzkovitz S, Milo R, Sigal A, Dekel E, Yarnitzky T, Liron Y, Polak P, Lahav G. Oscillations and variability in the p53 system. Mol Syst Biol. 2006; 2(1):2006–0033.
Purvis JE, Karhohs KW, Mock C, Batchelor E, Loewer A, Lahav G. p53 dynamics control cell fate. Sci Signal. 2012; 336(6087):1440.
Nelson D, Ihekwaba A, Elliott M, Johnson J, Gibney C, Foreman B, Nelson G, See V, Horton C, Spiller D, et al. Oscillations in NF κB signaling control the dynamics of gene expression. Sci Signal. 2004; 306(5696):704.
Ashall L, Horton CA, Nelson DE, Paszek P, Harper CV, Sillitoe K, Ryan S, Spiller DG, Unitt JF, Broomhead DS, et al. Pulsatile stimulation determines timing and specificity of NF κBdependent transcription. Sci Signal. 2009; 324(5924):242.
Tharp HS, Zhang WW. Robust control of a nonlinear timedelay system. Dyn Control. 1994; 4(1):21–38.
Nguyen LK. Regulation of oscillation dynamics in biochemical systems with dual negative feedback loops. J R Soc Interface. 2012; 9(73):1998–2010.
Lapytsko A, Schaber J. The role of time delay in adaptive cellular negative feedback systems. J Theor Biol. 2016; 398:64–73.
Lev BarOr R, Maya R, Segel LA, Alon U, Levine AJ, Oren M. Generation of oscillations by the p53Mdm2 feedback loop: A theoretical and experimental study. Proc Natl Acad Sci. 2000; 97(21):11250–11255.
Engelborghs K, Luzyanina T, Samaey G. Numerical bifurcation analysis of delay differential equations using DDEBIFTOOL. ACM Transactions on Mathematical Software, 28. 2002;:1–21.
Batchelor E, Loewer A, Mock C, Lahav G. Stimulusdependent dynamics of p53 in single cells. Mol Syst Biol. 2011; 7:488.
Kubbutat MH, Ludwig RL, Ashcroft M, Vousden KH. Regulation of Mdm2directed degradation by the C terminus of p53. Mol Cell Biol. 1998; 18(10):5690–698.
Chen X, Chen J, Gan S, Guan H, Zhou Y, Ouyang Q, Shi J. DNA damage strength modulates a bimodal switch of p53 dynamics for cellfate control. BMC Biol. 2013; 11(1):73.
Acknowledgements
The authors thank Prof. Dietrich Flockerzi for fruitful discussions and Prof. Steffen Waldherr for critical comments on the manuscript. The authors are also affiliated to the “International Max Planck Research School (IMPRS) for Advanced Methods in Process and System Engineering (Magdeburg)”.
Funding
This work was supported by the German Ministry of Science and Education [Federal Ministry of Education and Research project 0135779 to JS].
Availability of data and materials
The data sets supporting the result of the article are included within the article and its additional files.
Authors’ contributions
AB and JS designed research, AB did the mathematical and data analysis, JS and AB wrote the paper. Both authors read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Author information
Affiliations
Corresponding author
Additional file
Additional file 1
Supporting material. The file contains detailed stability analysis of Models 14; theoretical analysis showing how autoinhibition increases τ _{ m }; demonstration how τ _{ m } can be increased by autoinhibition for the p53 model; details of robustness analysis of the optimal solution for the p53 model; modelling the switch from oscillatory to adaptive response of the p53 system; calculating the duration of “On” and “Off” states of p53 pulses; figure demonstrating the dependence between parameter values of Models 14 and τ _{ m }; figure with results of MonteCarlo analysis of Models 14 applied to the alternative parameter set; table with the bestfit parameters for the p53 model. (PDF 1044 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Börsch, A., Schaber, J. How time delay and network design shape response patterns in biochemical negative feedback systems. BMC Syst Biol 10, 82 (2016). https://doi.org/10.1186/s1291801603259
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1291801603259
Keywords
 Stability
 Bifurcation
 Autoinhibiton
 Mass conservation
 p53