How time delay and network design shape response patterns in biochemical negative feedback systems

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 two-dimensional delay differential equations describing the dynamics of biochemical signal-response networks. Results We investigate the influence of several design features on the stability of the model equilibrium, i.e., presence of auto-inhibition and/or mass conservation and the kind and/or strength of the delayed negative feedback. We show that auto-inhibition and mass conservation have a stabilizing effect, whereas increasing abruptness and decreasing feedback threshold have a de-stabilizing effect on the model equilibrium. Moreover, applying our theoretical analysis to the mammalian p53 system we show that an auto-inhibitory 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. Electronic supplementary material The online version of this article (doi:10.1186/s12918-016-0325-9) contains supplementary material, which is available to authorized users.

According to the lemma from [2], as τ varies, the sum of the multiplicities of zeros of (2) in the open right half-plane M (τ ) can change only if a zero appears on or crosses the imaginary axis. Thus, the only way that M (τ ) = M (τ ) for τ < τ is, if there exists a marginal value τ m between τ and τ , such that P (λ m , τ m ) = 0 and Re(λ m ) = 0 [3].
We summed up above equations and applied Pythagorean trigonometric identity: Equation (6) has real non-zero roots if and only if x 2 β 2 − y 2 < 0, which is equivalent to x β < y. Otherwise, τ m does not exist and the equilibrium (C s , R s ) of Models 1-4 is absolutely stable. We assumed that x β < y holds and defined the discriminant D of (6) by Then, we defined roots of (6): Using the fact that ω 2 m should be positive, we get: Further, we substituted the obtained expression for ω m into (4) and (5): Using (8) and (9) we expressed the value of τ m in the following form: . In both cases, the smallest time delay, which causes a purely imaginary pair of roots λ 1,2 = ±iω m , is Further, we proved the stability of the equilibrium (C s , R s ) for any τ ∈ [0, τ m ). Indeed, roots of P (λ, 0) are either real negative or have negative real part (see above).
Then, according to Lemma from [2], M (τ ) = 0 for any τ ∈ [0, τ m ). Hence, for any τ ∈ [0, τ m ) the equilibrium (C s , R s ) is asymptotically stable. Now, we prove that a Hopf bifurcation occurs at τ = τ m . For this we differentiated the characteristic polynomial (2) with respect to τ and equated it to zero: Then, we obtained the expression for dλ dτ : We substituted λ = iω m into the expression for dλ dτ and applied Euler's formula: .
We substituted the expression for τ m (10), used (8) and (9), multiplied the top and the bottom by the conjugate of the denominator, took real part and obtained: The positivity of (11) guarantees that the hypotheses of the implicit function theorem hold. Hence, we may conclude that for τ ≈ τ m the root of the characteristic polynomial (2) λ m = iω m crosses the imaginary axis from left to the right.
According to Proposition 6.5 from [3], for τ ≥ τ m the equilibrium (C s , R s ) is unstable.

Auto-inhibition increases τ m
We represented τ m (10) as a function of x, y and β: where Assuming x β < y and using D from (7), we obtained derivatives of functions f and g with respect to x and y: Thus, we conclude that f (x, y, β) and g(x, y, β) are both positive functions, which increase with x and decrease with y. Consequently, τ m (x, y, β) follows the same pattern. Using the relation between equilibria of models with and without auto-inhibitory feedback and properties of feedback functions F , S 1 , S 2 described in the main manuscript, we calculated lower and upper bounds of x and y. Further, we showed how these bounds can be increased or decreased leading to changing in the value of τ m .
Firstly, we defined the relation between equilibrium components R s and C s for Models 1-4.
For this we equated the right hand side of the differential equation dR dt to 0: Further, for each Model 1-4 we equated the right hand side of the differential equation dC dt to 0 and substituted the expression for R s (14): Thus, for each model we got the equation defining the equilibrium component C s in the form: According to (19), we find the equilibrium component C s as the intersection of functions θ 1 • F and θ 2 (see Fig. S1). Since θ 1 is either constant or decreasing function with respect to C and θ 2 is increasing function with respect to C, the equilibrium component C s always exists and is unique.
In case the auto-inhibitory feedback is not present in the system (F (C) ≡ 1), the equilibrium component C s of the model without auto-inhibition is always greater than the equilibrium component C s of the model with auto-inhibition (see Fig. S1). According to (19), the same conclusion holds for equilibrium components R s and R s of models with and without auto-inhibition, respectively.
Figure S1 Schematic intersection of functions θ 1 , θ 1 • F and θ 2 . Functions θ 1 and θ 2 were introduced to find equilibrium components of Models 1-4. a How to find equilibrium components C s , C s for Model 3. b How to find equilibrium components C s , C s for Models 1, 2, 4.
Taking into account obtained relations C s ≤ C s and R s ≤ R s and properties of feedback functions F , S 1 , S 2 described in the main manuscript we obtained following estimations of x and y: Refer to Table S2 for values of ε lb , ε ub , σ lb , σ ub for each Model 1-4. Both the lower and upper bound of x, i.e., ε lb and ε ub , are increasing with |F (C s )| for Models 1-4. Therefore, we can always increase a given x by choosing an appropriate value for |F (C s )|. The lower and upper bound of y, i.e., σ lb and σ ub , have non-negative constant values. Consequently, according to (13), we can always increase τ m by increasing |F (C s )|.
Taken together, we show how auto-inhibitory feedback allows to modify the range of the interval [0, τ m ).

Numerical evidence of increasing τ m by auto-inhibition
The theoretical analysis presented in the section above shows that increasing the slope of the autoinhibitory function at the equilibrium |F (C s )| leads to increasing the marginal time delay τ m . Here, we suggest an algorithm how to increase |F (C s )| and, consequently, τ m for the p53 model.
As auto-inhibitory feedback function we used a reverse Hill function F H (C) (Fig. S2a): For the reverse Hill-function, an obvious choice to increase the slope is increasing the Hillcoefficient ν. Therefore, we made ν a free positive parameter, which value we choose. For simplicity we considered only integer values of ν, however the algorithm can be extended to real values of ν. Then we adjusted values of κ and equilibrium (C s , R s ) and maximized |F H (C s )|, which is equivalent to solving F H (C s ) = 0. Thus, for the fixed value of ν we solved the following system of equations with respect to (C s , R s , κ) using fitted parameter values (α, β, δ, K m , n) from Table S1: Then, we quantified τ m using (12). If it is necessary, we increase ν, solve the system (22) with respect to (C s , R s , κ) and quantify τ m again. Fig. S3 demonstrates results of application of the algorithm to the p53 model. In Fig. S3a we illustrated the graph of the slope |F H (C)| for several values of ν and adjusted values of κ. One can see that |F H (C s )| is increasing, when we increase ν, whereas the value of C s is only slightly decreasing (indicated by black dots in Fig. S3a). Consequently, τ m is also increasing with respect to the Hill coefficient ν (Fig. S3b). For ν ≥ 3 the equilibrium is asymptotically stable. For ν > 8 the equilibrium is absolutely stable. Thus, the relation between ν, |F H (C)| and τ m demonstrated in Fig. S3 corresponds well to the theoretical analysis made above.
In Fig. S4 we depicted simulations of the p53 model with fitted parameters from Table S1 and synthetically activated auto-inhibitory feedback F H (C) (21) with ν = 3, κ = 1.73. One can see that the model with auto-inhibitory feedback produces damped instead of sustained oscillations.  Simulation of the p53 model with fitted parameters from Table S1 and synthetically activated auto-inhibitory feedback F H (C) (21) with ν = 3, κ = 1.73.
Further, we investigated the stability of the p53 model with and without synthetically activated feedback with respect to several parameters. In Fig. S5a-c we show that increasing the Hill coefficient ν and adjusting κ leads to a decreasing parameter region, where oscillations occur. Fig. S5d demonstrates that τ m increases with respect to κ for both ν = 2 (dark gray) and ν = 3 (light grey). Thus, these results correspond well to our previous research [4], where we proved for the model with delayed negative feedback that nested negative feedbacks lead to increasing resistance of the system. In addition, in Fig. S5a we illustrate that τ m increases with respect to α that indicates the stabilizing property of α and may explain such a small fitted parameter value. On the other hand, τ m decreases with respect to I (Fig. S5c).

Figure S5
Dependencies between parameters of the p53 model. The p53 model is considered without (black) and with synthetically activated (ν = 2, κ = 1.23 dark grey, ν = 3, κ = 1.73 light grey) auto-inhibitory feedback. a τ m and α. b τ m and β. c τ m and I. d τ m and κ. Designations: dots -the fitted value of the parameter used in the p53 model, dashed line -the fitted value of τ (Table S1).

Robustness of optimal solution of p53 model
We analysed the robustness of the optimal solutions for the p53 model with respect to noise. To this end, we randomly sampled parameter values within ±10% of their respective fitted values using a uniform distribution for 100 times. Then we simulated p53 model with perturbed parameters and calculated 0.05 and 0.95 quantiles of obtained simulations (grey regions in Fig. S6). Fig. S6 shows that the model solution with fitted parameters is located between 0.05 and 0.95 quantiles for all considered time points.
Further, for each perturbed parameter set we calculated the relative variation of the integral of the first transient response after the stimulation: where C(t) corresponds to the model solution with fitted parameters, C i p (t) corresponds to the model solution with i th perturbed parameter set (i = 1, 2, 3, ..., 100), t Int corresponds to the time of the first minimum after initial stimulation (Fig. S6).
This way, the robustness of both initial activation amplitude and timing of the first transient response, two characteristic measures of system dynamics, can be estimated concomitantly.
As the measure of robustness we calculated the mean value < var > and standard deviation sd of obtained relative variations var i : < var > ±sd = 8.7 ± 6.4%.
One may see that for the considered model values of < var > and sd do not exceed 9%. Thus, the fitted solution turned out to be very robust with respect to noise in the parameters, also indicating that the fitted solution is in a well-defined local minimum.
Calculating the duration of "On" and "Off" states of p53 pulses For calculating the duration of "On" and "Off" states of p53 pulses we simulated p53 model (3) with parameters from Table S1 for later time points (t > 400 h), when the system reached its steady state (see Fig. S9). Then we quantified the mean value between maxima and minima of the simulation curve (see black dashed line in Fig. S9). As a result the simulation curve was split on upper and lower parts with respect to the mean value. We designated the upper part as "On" state and the lower part as "Off" state. Figure S9 Defining "On" and "Off" states for the p53 response. Black dashed line designates the mean value between maxima and minima of the simulation curve.
Then we checked how different model parameters control the duration of "On" and "Off" states (see Fig. S10). For this we varied parameter values I, τ , K m , n, α, δ, β, κ and ν one at a time in the range, where p53 model (3) produces sustained oscillations. Then for each varied parameter value we calculated the duration of "On" and "Off" states of the p53 response as described above.